Mercurial > hg > octave-nkf

comparison liboctave/EIG.cc @ 462:07fabd96ac6a

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 1994-06-06 00:58:18 by jwe] Initial revision
author jwe
date Mon, 06 Jun 1994 00:58:18 +0000
parents
children 714fd17fca28
comparison
equal deleted inserted replaced
461:00f8b2242a18 462:07fabd96ac6a
1 // -*- C++ -*-
2 /*
3
4 Copyright (C) 1992, 1993, 1994 John W. Eaton
5
6 This file is part of Octave.
7
8 Octave is free software; you can redistribute it and/or modify it
9 under the terms of the GNU General Public License as published by the
10 Free Software Foundation; either version 2, or (at your option) any
11 later version.
12
13 Octave is distributed in the hope that it will be useful, but WITHOUT
14 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
15 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
16 for more details.
17
18 You should have received a copy of the GNU General Public License
19 along with Octave; see the file COPYING. If not, write to the Free
20 Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
21
22 */
23
24 #ifdef HAVE_CONFIG_H
25 #include "config.h"
26 #endif
27
28 #if defined (__GNUG__)
29 #pragma implementation
30 #endif
31
32 #include "EIG.h"
33 #include "mx-inlines.cc"
34 #include "lo-error.h"
35 #include "f77-uscore.h"
36
37 extern "C"
38 {
39 int F77_FCN (dgeev) (const char*, const char*, const int*, double*,
40 const int*, double*, double*, double*,
41 const int*, double*, const int*, double*,
42 const int*, int*, long, long);
43
44 int F77_FCN (zgeev) (const char*, const char*, const int*, Complex*,
45 const int*, Complex*, Complex*, const int*,
46 Complex*, const int*, Complex*, const int*,
47 double*, int*, long, long);
48 }
49
50 int
51 EIG::init (const Matrix& a)
52 {
53 int a_nr = a.rows ();
54 if (a_nr != a.cols ())
55 {
56 (*current_liboctave_error_handler) ("EIG requires square matrix");
57 return -1;
58 }
59
60 int n = a_nr;
61
62 int info;
63
64 char jobvl = 'N';
65 char jobvr = 'V';
66
67 double *tmp_data = dup (a.data (), a.length ());
68 double *wr = new double[n];
69 double *wi = new double[n];
70 Matrix vr (n, n);
71 double *pvr = vr.fortran_vec ();
72 int lwork = 8*n;
73 double *work = new double[lwork];
74
75 double dummy;
76 int idummy = 1;
77
78 F77_FCN (dgeev) (&jobvl, &jobvr, &n, tmp_data, &n, wr, wi, &dummy,
79 &idummy, pvr, &n, work, &lwork, &info, 1L, 1L);
80
81 lambda.resize (n);
82 v.resize (n, n);
83
84 for (int j = 0; j < n; j++)
85 {
86 if (wi[j] == 0.0)
87 {
88 lambda.elem (j) = Complex (wr[j]);
89 for (int i = 0; i < n; i++)
90 v.elem (i, j) = vr.elem (i, j);
91 }
92 else
93 {
94 if (j+1 >= n)
95 {
96 (*current_liboctave_error_handler) ("EIG: internal error");
97 return -1;
98 }
99
100 for (int i = 0; i < n; i++)
101 {
102 lambda.elem (j) = Complex (wr[j], wi[j]);
103 lambda.elem (j+1) = Complex (wr[j+1], wi[j+1]);
104 double real_part = vr.elem (i, j);
105 double imag_part = vr.elem (i, j+1);
106 v.elem (i, j) = Complex (real_part, imag_part);
107 v.elem (i, j+1) = Complex (real_part, -imag_part);
108 }
109 j++;
110 }
111 }
112
113 delete [] tmp_data;
114 delete [] wr;
115 delete [] wi;
116 delete [] work;
117
118 return info;
119 }
120
121 int
122 EIG::init (const ComplexMatrix& a)
123 {
124 int a_nr = a.rows ();
125 if (a_nr != a.cols ())
126 {
127 (*current_liboctave_error_handler) ("EIG requires square matrix");
128 return -1;
129 }
130
131 int n = a_nr;
132
133 int info;
134
135 char jobvl = 'N';
136 char jobvr = 'V';
137
138 lambda.resize (n);
139 v.resize (n, n);
140
141 Complex *pw = lambda.fortran_vec ();
142 Complex *pvr = v.fortran_vec ();
143
144 Complex *tmp_data = dup (a.data (), a.length ());
145
146 int lwork = 8*n;
147 Complex *work = new Complex[lwork];
148 double *rwork = new double[4*n];
149
150 Complex dummy;
151 int idummy = 1;
152
153 F77_FCN (zgeev) (&jobvl, &jobvr, &n, tmp_data, &n, pw, &dummy,
154 &idummy, pvr, &n, work, &lwork, rwork, &info, 1L,
155 1L);
156
157 delete [] tmp_data;
158 delete [] work;
159 delete [] rwork;
160
161 return info;
162 }
163
164 /*
165 ;;; Local Variables: ***
166 ;;; mode: C++ ***
167 ;;; page-delimiter: "^/\\*" ***
168 ;;; End: ***
169 */