Mercurial > octave-nkf
annotate liboctave/EIG.cc @ 11518:141b3fb5cef7
style fixes
author | John W. Eaton <jwe@octave.org> |
---|---|
date | Thu, 13 Jan 2011 16:52:30 -0500 |
parents | 8a5e980da6aa |
children | fd0a3ac60b0e |
rev | line source |
---|---|
462 | 1 /* |
2 | |
7017 | 3 Copyright (C) 1994, 1995, 1996, 1997, 1998, 2000, 2002, 2003, 2004, |
8920 | 4 2005, 2006, 2007, 2008 John W. Eaton |
462 | 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 | |
7016 | 10 Free Software Foundation; either version 3 of the License, or (at your |
11 option) any later version. | |
462 | 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 | |
7016 | 19 along with Octave; see the file COPYING. If not, see |
20 <http://www.gnu.org/licenses/>. | |
462 | 21 |
22 */ | |
23 | |
24 #ifdef HAVE_CONFIG_H | |
1192 | 25 #include <config.h> |
462 | 26 #endif |
27 | |
28 #include "EIG.h" | |
2815 | 29 #include "dColVector.h" |
1847 | 30 #include "f77-fcn.h" |
462 | 31 #include "lo-error.h" |
32 | |
33 extern "C" | |
34 { | |
4552 | 35 F77_RET_T |
36 F77_FUNC (dgeev, DGEEV) (F77_CONST_CHAR_ARG_DECL, | |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
37 F77_CONST_CHAR_ARG_DECL, |
11495 | 38 const octave_idx_type&, double*, |
39 const octave_idx_type&, double*, double*, | |
40 double*, const octave_idx_type&, double*, | |
41 const octave_idx_type&, double*, | |
42 const octave_idx_type&, octave_idx_type& | |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
43 F77_CHAR_ARG_LEN_DECL |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
44 F77_CHAR_ARG_LEN_DECL); |
462 | 45 |
4552 | 46 F77_RET_T |
47 F77_FUNC (zgeev, ZGEEV) (F77_CONST_CHAR_ARG_DECL, | |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
48 F77_CONST_CHAR_ARG_DECL, |
11495 | 49 const octave_idx_type&, Complex*, |
50 const octave_idx_type&, Complex*, | |
51 Complex*, const octave_idx_type&, Complex*, | |
52 const octave_idx_type&, Complex*, | |
53 const octave_idx_type&, double*, octave_idx_type& | |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
54 F77_CHAR_ARG_LEN_DECL |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
55 F77_CHAR_ARG_LEN_DECL); |
2815 | 56 |
4552 | 57 F77_RET_T |
58 F77_FUNC (dsyev, DSYEV) (F77_CONST_CHAR_ARG_DECL, | |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
59 F77_CONST_CHAR_ARG_DECL, |
11495 | 60 const octave_idx_type&, double*, |
61 const octave_idx_type&, double*, double*, | |
62 const octave_idx_type&, octave_idx_type& | |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
63 F77_CHAR_ARG_LEN_DECL |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
64 F77_CHAR_ARG_LEN_DECL); |
2815 | 65 |
4552 | 66 F77_RET_T |
67 F77_FUNC (zheev, ZHEEV) (F77_CONST_CHAR_ARG_DECL, | |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
68 F77_CONST_CHAR_ARG_DECL, |
11495 | 69 const octave_idx_type&, Complex*, |
70 const octave_idx_type&, double*, | |
71 Complex*, const octave_idx_type&, double*, | |
72 octave_idx_type& | |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
73 F77_CHAR_ARG_LEN_DECL |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
74 F77_CHAR_ARG_LEN_DECL); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
75 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
76 F77_RET_T |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
77 F77_FUNC (dpotrf, DPOTRF) (F77_CONST_CHAR_ARG_DECL, |
11495 | 78 const octave_idx_type&, double*, |
79 const octave_idx_type&, octave_idx_type& | |
80 F77_CHAR_ARG_LEN_DECL | |
81 F77_CHAR_ARG_LEN_DECL); | |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
82 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
83 F77_RET_T |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
84 F77_FUNC (zpotrf, ZPOTRF) (F77_CONST_CHAR_ARG_DECL, |
11495 | 85 const octave_idx_type&, |
86 Complex*, const octave_idx_type&, | |
87 octave_idx_type& | |
88 F77_CHAR_ARG_LEN_DECL | |
89 F77_CHAR_ARG_LEN_DECL); | |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
90 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
91 F77_RET_T |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
92 F77_FUNC (dggev, DGGEV) (F77_CONST_CHAR_ARG_DECL, |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
93 F77_CONST_CHAR_ARG_DECL, |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
94 const octave_idx_type&, |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
95 double*, const octave_idx_type&, |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
96 double*, const octave_idx_type&, |
11495 | 97 double*, double*, double *, double*, |
98 const octave_idx_type&, double*, | |
99 const octave_idx_type&, double*, | |
100 const octave_idx_type&, octave_idx_type& | |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
101 F77_CHAR_ARG_LEN_DECL |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
102 F77_CHAR_ARG_LEN_DECL); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
103 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
104 F77_RET_T |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
105 F77_FUNC (dsygv, DSYGV) (const octave_idx_type&, |
11518 | 106 F77_CONST_CHAR_ARG_DECL, |
107 F77_CONST_CHAR_ARG_DECL, | |
11495 | 108 const octave_idx_type&, double*, |
109 const octave_idx_type&, double*, | |
110 const octave_idx_type&, double*, double*, | |
111 const octave_idx_type&, octave_idx_type& | |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
112 F77_CHAR_ARG_LEN_DECL |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
113 F77_CHAR_ARG_LEN_DECL); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
114 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
115 F77_RET_T |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
116 F77_FUNC (zggev, ZGGEV) (F77_CONST_CHAR_ARG_DECL, |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
117 F77_CONST_CHAR_ARG_DECL, |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
118 const octave_idx_type&, |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
119 Complex*, const octave_idx_type&, |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
120 Complex*, const octave_idx_type&, |
11495 | 121 Complex*, Complex*, Complex*, |
122 const octave_idx_type&, Complex*, | |
123 const octave_idx_type&, Complex*, | |
124 const octave_idx_type&, double*, octave_idx_type& | |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
125 F77_CHAR_ARG_LEN_DECL |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
126 F77_CHAR_ARG_LEN_DECL); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
127 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
128 F77_RET_T |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
129 F77_FUNC (zhegv, ZHEGV) (const octave_idx_type&, |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
130 F77_CONST_CHAR_ARG_DECL, |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
131 F77_CONST_CHAR_ARG_DECL, |
11495 | 132 const octave_idx_type&, Complex*, |
133 const octave_idx_type&, Complex*, | |
134 const octave_idx_type&, double*, Complex*, | |
135 const octave_idx_type&, double*, octave_idx_type& | |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
136 F77_CHAR_ARG_LEN_DECL |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
137 F77_CHAR_ARG_LEN_DECL); |
462 | 138 } |
139 | |
5275 | 140 octave_idx_type |
4725 | 141 EIG::init (const Matrix& a, bool calc_ev) |
462 | 142 { |
5822 | 143 if (a.any_element_is_inf_or_nan ()) |
144 { | |
145 (*current_liboctave_error_handler) | |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
146 ("EIG: matrix contains Inf or NaN values"); |
5822 | 147 return -1; |
148 } | |
149 | |
2815 | 150 if (a.is_symmetric ()) |
4725 | 151 return symmetric_init (a, calc_ev); |
2815 | 152 |
5275 | 153 octave_idx_type n = a.rows (); |
1934 | 154 |
155 if (n != a.cols ()) | |
462 | 156 { |
157 (*current_liboctave_error_handler) ("EIG requires square matrix"); | |
158 return -1; | |
159 } | |
160 | |
5275 | 161 octave_idx_type info = 0; |
462 | 162 |
1934 | 163 Matrix atmp = a; |
164 double *tmp_data = atmp.fortran_vec (); | |
165 | |
10350
12884915a8e4
merge MArray classes & improve Array interface
Jaroslav Hajek <highegg@gmail.com>
parents:
10314
diff
changeset
|
166 Array<double> wr (n, 1); |
1934 | 167 double *pwr = wr.fortran_vec (); |
168 | |
10350
12884915a8e4
merge MArray classes & improve Array interface
Jaroslav Hajek <highegg@gmail.com>
parents:
10314
diff
changeset
|
169 Array<double> wi (n, 1); |
1934 | 170 double *pwi = wi.fortran_vec (); |
171 | |
9020
728e7943752d
EIG.cc: avoid volatile decl for tmp variable
John W. Eaton <jwe@octave.org>
parents:
8920
diff
changeset
|
172 octave_idx_type tnvr = calc_ev ? n : 0; |
728e7943752d
EIG.cc: avoid volatile decl for tmp variable
John W. Eaton <jwe@octave.org>
parents:
8920
diff
changeset
|
173 Matrix vr (tnvr, tnvr); |
462 | 174 double *pvr = vr.fortran_vec (); |
1934 | 175 |
5275 | 176 octave_idx_type lwork = -1; |
4800 | 177 double dummy_work; |
462 | 178 |
1365 | 179 double *dummy = 0; |
5275 | 180 octave_idx_type idummy = 1; |
462 | 181 |
4552 | 182 F77_XFCN (dgeev, DGEEV, (F77_CONST_CHAR_ARG2 ("N", 1), |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
183 F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
184 n, tmp_data, n, pwr, pwi, dummy, |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
185 idummy, pvr, n, &dummy_work, lwork, info |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
186 F77_CHAR_ARG_LEN (1) |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
187 F77_CHAR_ARG_LEN (1))); |
462 | 188 |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
189 if (info == 0) |
462 | 190 { |
5275 | 191 lwork = static_cast<octave_idx_type> (dummy_work); |
10350
12884915a8e4
merge MArray classes & improve Array interface
Jaroslav Hajek <highegg@gmail.com>
parents:
10314
diff
changeset
|
192 Array<double> work (lwork, 1); |
4800 | 193 double *pwork = work.fortran_vec (); |
194 | |
195 F77_XFCN (dgeev, DGEEV, (F77_CONST_CHAR_ARG2 ("N", 1), | |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
196 F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
197 n, tmp_data, n, pwr, pwi, dummy, |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
198 idummy, pvr, n, pwork, lwork, info |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
199 F77_CHAR_ARG_LEN (1) |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
200 F77_CHAR_ARG_LEN (1))); |
4800 | 201 |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
202 if (info < 0) |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
203 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
204 (*current_liboctave_error_handler) ("unrecoverable error in dgeev"); |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
205 return info; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
206 } |
4800 | 207 |
2815 | 208 if (info > 0) |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
209 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
210 (*current_liboctave_error_handler) ("dgeev failed to converge"); |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
211 return info; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
212 } |
1934 | 213 |
4800 | 214 lambda.resize (n); |
9020
728e7943752d
EIG.cc: avoid volatile decl for tmp variable
John W. Eaton <jwe@octave.org>
parents:
8920
diff
changeset
|
215 octave_idx_type nvr = calc_ev ? n : 0; |
4800 | 216 v.resize (nvr, nvr); |
217 | |
5275 | 218 for (octave_idx_type j = 0; j < n; j++) |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
219 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
220 if (wi.elem (j) == 0.0) |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
221 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
222 lambda.elem (j) = Complex (wr.elem (j)); |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
223 for (octave_idx_type i = 0; i < nvr; i++) |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
224 v.elem (i, j) = vr.elem (i, j); |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
225 } |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
226 else |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
227 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
228 if (j+1 >= n) |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
229 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
230 (*current_liboctave_error_handler) ("EIG: internal error"); |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
231 return -1; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
232 } |
4800 | 233 |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
234 lambda.elem(j) = Complex (wr.elem(j), wi.elem(j)); |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
235 lambda.elem(j+1) = Complex (wr.elem(j+1), wi.elem(j+1)); |
2815 | 236 |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
237 for (octave_idx_type i = 0; i < nvr; i++) |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
238 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
239 double real_part = vr.elem (i, j); |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
240 double imag_part = vr.elem (i, j+1); |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
241 v.elem (i, j) = Complex (real_part, imag_part); |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
242 v.elem (i, j+1) = Complex (real_part, -imag_part); |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
243 } |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
244 j++; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
245 } |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
246 } |
2815 | 247 } |
4800 | 248 else |
249 (*current_liboctave_error_handler) ("dgeev workspace query failed"); | |
2815 | 250 |
251 return info; | |
252 } | |
253 | |
5275 | 254 octave_idx_type |
4725 | 255 EIG::symmetric_init (const Matrix& a, bool calc_ev) |
2815 | 256 { |
5275 | 257 octave_idx_type n = a.rows (); |
1934 | 258 |
2815 | 259 if (n != a.cols ()) |
260 { | |
261 (*current_liboctave_error_handler) ("EIG requires square matrix"); | |
262 return -1; | |
263 } | |
264 | |
5275 | 265 octave_idx_type info = 0; |
2815 | 266 |
267 Matrix atmp = a; | |
268 double *tmp_data = atmp.fortran_vec (); | |
269 | |
3585 | 270 ColumnVector wr (n); |
2815 | 271 double *pwr = wr.fortran_vec (); |
272 | |
5275 | 273 octave_idx_type lwork = -1; |
4800 | 274 double dummy_work; |
2815 | 275 |
4725 | 276 F77_XFCN (dsyev, DSYEV, (F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
277 F77_CONST_CHAR_ARG2 ("U", 1), |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
278 n, tmp_data, n, pwr, &dummy_work, lwork, info |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
279 F77_CHAR_ARG_LEN (1) |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
280 F77_CHAR_ARG_LEN (1))); |
2815 | 281 |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
282 if (info == 0) |
2815 | 283 { |
5275 | 284 lwork = static_cast<octave_idx_type> (dummy_work); |
10350
12884915a8e4
merge MArray classes & improve Array interface
Jaroslav Hajek <highegg@gmail.com>
parents:
10314
diff
changeset
|
285 Array<double> work (lwork, 1); |
4800 | 286 double *pwork = work.fortran_vec (); |
287 | |
288 F77_XFCN (dsyev, DSYEV, (F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), | |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
289 F77_CONST_CHAR_ARG2 ("U", 1), |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
290 n, tmp_data, n, pwr, pwork, lwork, info |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
291 F77_CHAR_ARG_LEN (1) |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
292 F77_CHAR_ARG_LEN (1))); |
4800 | 293 |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
294 if (info < 0) |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
295 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
296 (*current_liboctave_error_handler) ("unrecoverable error in dsyev"); |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
297 return info; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
298 } |
4800 | 299 |
300 if (info > 0) | |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
301 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
302 (*current_liboctave_error_handler) ("dsyev failed to converge"); |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
303 return info; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
304 } |
4800 | 305 |
3585 | 306 lambda = ComplexColumnVector (wr); |
4725 | 307 v = calc_ev ? ComplexMatrix (atmp) : ComplexMatrix (); |
462 | 308 } |
4800 | 309 else |
310 (*current_liboctave_error_handler) ("dsyev workspace query failed"); | |
462 | 311 |
312 return info; | |
313 } | |
314 | |
5275 | 315 octave_idx_type |
4725 | 316 EIG::init (const ComplexMatrix& a, bool calc_ev) |
462 | 317 { |
5822 | 318 if (a.any_element_is_inf_or_nan ()) |
319 { | |
320 (*current_liboctave_error_handler) | |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
321 ("EIG: matrix contains Inf or NaN values"); |
5822 | 322 return -1; |
323 } | |
324 | |
2815 | 325 if (a.is_hermitian ()) |
4725 | 326 return hermitian_init (a, calc_ev); |
2815 | 327 |
5275 | 328 octave_idx_type n = a.rows (); |
1934 | 329 |
330 if (n != a.cols ()) | |
462 | 331 { |
332 (*current_liboctave_error_handler) ("EIG requires square matrix"); | |
333 return -1; | |
334 } | |
335 | |
5275 | 336 octave_idx_type info = 0; |
462 | 337 |
1934 | 338 ComplexMatrix atmp = a; |
339 Complex *tmp_data = atmp.fortran_vec (); | |
462 | 340 |
2815 | 341 ComplexColumnVector w (n); |
342 Complex *pw = w.fortran_vec (); | |
343 | |
5275 | 344 octave_idx_type nvr = calc_ev ? n : 0; |
4725 | 345 ComplexMatrix vtmp (nvr, nvr); |
2815 | 346 Complex *pv = vtmp.fortran_vec (); |
347 | |
5275 | 348 octave_idx_type lwork = -1; |
4800 | 349 Complex dummy_work; |
1934 | 350 |
5275 | 351 octave_idx_type lrwork = 2*n; |
10350
12884915a8e4
merge MArray classes & improve Array interface
Jaroslav Hajek <highegg@gmail.com>
parents:
10314
diff
changeset
|
352 Array<double> rwork (lrwork, 1); |
1934 | 353 double *prwork = rwork.fortran_vec (); |
462 | 354 |
1365 | 355 Complex *dummy = 0; |
5275 | 356 octave_idx_type idummy = 1; |
462 | 357 |
4552 | 358 F77_XFCN (zgeev, ZGEEV, (F77_CONST_CHAR_ARG2 ("N", 1), |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
359 F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
360 n, tmp_data, n, pw, dummy, idummy, |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
361 pv, n, &dummy_work, lwork, prwork, info |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
362 F77_CHAR_ARG_LEN (1) |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
363 F77_CHAR_ARG_LEN (1))); |
2815 | 364 |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
365 if (info == 0) |
2815 | 366 { |
5275 | 367 lwork = static_cast<octave_idx_type> (dummy_work.real ()); |
10350
12884915a8e4
merge MArray classes & improve Array interface
Jaroslav Hajek <highegg@gmail.com>
parents:
10314
diff
changeset
|
368 Array<Complex> work (lwork, 1); |
4800 | 369 Complex *pwork = work.fortran_vec (); |
370 | |
371 F77_XFCN (zgeev, ZGEEV, (F77_CONST_CHAR_ARG2 ("N", 1), | |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
372 F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
373 n, tmp_data, n, pw, dummy, idummy, |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
374 pv, n, pwork, lwork, prwork, info |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
375 F77_CHAR_ARG_LEN (1) |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
376 F77_CHAR_ARG_LEN (1))); |
4800 | 377 |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
378 if (info < 0) |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
379 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
380 (*current_liboctave_error_handler) ("unrecoverable error in zgeev"); |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
381 return info; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
382 } |
4800 | 383 |
384 if (info > 0) | |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
385 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
386 (*current_liboctave_error_handler) ("zgeev failed to converge"); |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
387 return info; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
388 } |
4800 | 389 |
2815 | 390 lambda = w; |
391 v = vtmp; | |
392 } | |
4800 | 393 else |
394 (*current_liboctave_error_handler) ("zgeev workspace query failed"); | |
2815 | 395 |
396 return info; | |
397 } | |
398 | |
5275 | 399 octave_idx_type |
4725 | 400 EIG::hermitian_init (const ComplexMatrix& a, bool calc_ev) |
2815 | 401 { |
5275 | 402 octave_idx_type n = a.rows (); |
2815 | 403 |
404 if (n != a.cols ()) | |
405 { | |
406 (*current_liboctave_error_handler) ("EIG requires square matrix"); | |
407 return -1; | |
408 } | |
409 | |
5275 | 410 octave_idx_type info = 0; |
462 | 411 |
2815 | 412 ComplexMatrix atmp = a; |
413 Complex *tmp_data = atmp.fortran_vec (); | |
414 | |
3585 | 415 ColumnVector wr (n); |
416 double *pwr = wr.fortran_vec (); | |
2815 | 417 |
5275 | 418 octave_idx_type lwork = -1; |
4800 | 419 Complex dummy_work; |
2815 | 420 |
5275 | 421 octave_idx_type lrwork = 3*n; |
10350
12884915a8e4
merge MArray classes & improve Array interface
Jaroslav Hajek <highegg@gmail.com>
parents:
10314
diff
changeset
|
422 Array<double> rwork (lrwork, 1); |
2815 | 423 double *prwork = rwork.fortran_vec (); |
424 | |
4725 | 425 F77_XFCN (zheev, ZHEEV, (F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
426 F77_CONST_CHAR_ARG2 ("U", 1), |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
427 n, tmp_data, n, pwr, &dummy_work, lwork, |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
428 prwork, info |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
429 F77_CHAR_ARG_LEN (1) |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
430 F77_CHAR_ARG_LEN (1))); |
2815 | 431 |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
432 if (info == 0) |
2815 | 433 { |
5275 | 434 lwork = static_cast<octave_idx_type> (dummy_work.real ()); |
10350
12884915a8e4
merge MArray classes & improve Array interface
Jaroslav Hajek <highegg@gmail.com>
parents:
10314
diff
changeset
|
435 Array<Complex> work (lwork, 1); |
4800 | 436 Complex *pwork = work.fortran_vec (); |
437 | |
438 F77_XFCN (zheev, ZHEEV, (F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), | |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
439 F77_CONST_CHAR_ARG2 ("U", 1), |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
440 n, tmp_data, n, pwr, pwork, lwork, prwork, info |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
441 F77_CHAR_ARG_LEN (1) |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
442 F77_CHAR_ARG_LEN (1))); |
4800 | 443 |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
444 if (info < 0) |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
445 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
446 (*current_liboctave_error_handler) ("unrecoverable error in zheev"); |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
447 return info; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
448 } |
4800 | 449 |
450 if (info > 0) | |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
451 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
452 (*current_liboctave_error_handler) ("zheev failed to converge"); |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
453 return info; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
454 } |
4800 | 455 |
3585 | 456 lambda = ComplexColumnVector (wr); |
4725 | 457 v = calc_ev ? ComplexMatrix (atmp) : ComplexMatrix (); |
2815 | 458 } |
4800 | 459 else |
460 (*current_liboctave_error_handler) ("zheev workspace query failed"); | |
462 | 461 |
462 return info; | |
463 } | |
464 | |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
465 octave_idx_type |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
466 EIG::init (const Matrix& a, const Matrix& b, bool calc_ev) |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
467 { |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
468 if (a.any_element_is_inf_or_nan () || b.any_element_is_inf_or_nan ()) |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
469 { |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
470 (*current_liboctave_error_handler) |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
471 ("EIG: matrix contains Inf or NaN values"); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
472 return -1; |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
473 } |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
474 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
475 octave_idx_type n = a.rows (); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
476 octave_idx_type nb = b.rows (); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
477 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
478 if (n != a.cols () || nb != b.cols ()) |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
479 { |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
480 (*current_liboctave_error_handler) ("EIG requires square matrix"); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
481 return -1; |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
482 } |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
483 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
484 if (n != nb) |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
485 { |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
486 (*current_liboctave_error_handler) ("EIG requires same size matrices"); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
487 return -1; |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
488 } |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
489 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
490 octave_idx_type info = 0; |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
491 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
492 Matrix tmp = b; |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
493 double *tmp_data = tmp.fortran_vec (); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
494 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
495 F77_XFCN (dpotrf, DPOTRF, (F77_CONST_CHAR_ARG2 ("L", 1), |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
496 n, tmp_data, n, |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
497 info |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
498 F77_CHAR_ARG_LEN (1) |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
499 F77_CHAR_ARG_LEN (1))); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
500 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
501 if (a.is_symmetric () && b.is_symmetric () && info == 0) |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
502 return symmetric_init (a, b, calc_ev); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
503 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
504 Matrix atmp = a; |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
505 double *atmp_data = atmp.fortran_vec (); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
506 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
507 Matrix btmp = b; |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
508 double *btmp_data = btmp.fortran_vec (); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
509 |
10350
12884915a8e4
merge MArray classes & improve Array interface
Jaroslav Hajek <highegg@gmail.com>
parents:
10314
diff
changeset
|
510 Array<double> ar (n, 1); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
511 double *par = ar.fortran_vec (); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
512 |
10350
12884915a8e4
merge MArray classes & improve Array interface
Jaroslav Hajek <highegg@gmail.com>
parents:
10314
diff
changeset
|
513 Array<double> ai (n, 1); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
514 double *pai = ai.fortran_vec (); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
515 |
10350
12884915a8e4
merge MArray classes & improve Array interface
Jaroslav Hajek <highegg@gmail.com>
parents:
10314
diff
changeset
|
516 Array<double> beta (n, 1); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
517 double *pbeta = beta.fortran_vec (); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
518 |
9020
728e7943752d
EIG.cc: avoid volatile decl for tmp variable
John W. Eaton <jwe@octave.org>
parents:
8920
diff
changeset
|
519 octave_idx_type tnvr = calc_ev ? n : 0; |
728e7943752d
EIG.cc: avoid volatile decl for tmp variable
John W. Eaton <jwe@octave.org>
parents:
8920
diff
changeset
|
520 Matrix vr (tnvr, tnvr); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
521 double *pvr = vr.fortran_vec (); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
522 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
523 octave_idx_type lwork = -1; |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
524 double dummy_work; |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
525 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
526 double *dummy = 0; |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
527 octave_idx_type idummy = 1; |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
528 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
529 F77_XFCN (dggev, DGGEV, (F77_CONST_CHAR_ARG2 ("N", 1), |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
530 F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
531 n, atmp_data, n, btmp_data, n, |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
532 par, pai, pbeta, |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
533 dummy, idummy, pvr, n, |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
534 &dummy_work, lwork, info |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
535 F77_CHAR_ARG_LEN (1) |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
536 F77_CHAR_ARG_LEN (1))); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
537 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
538 if (info == 0) |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
539 { |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
540 lwork = static_cast<octave_idx_type> (dummy_work); |
10350
12884915a8e4
merge MArray classes & improve Array interface
Jaroslav Hajek <highegg@gmail.com>
parents:
10314
diff
changeset
|
541 Array<double> work (lwork, 1); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
542 double *pwork = work.fortran_vec (); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
543 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
544 F77_XFCN (dggev, DGGEV, (F77_CONST_CHAR_ARG2 ("N", 1), |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
545 F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
546 n, atmp_data, n, btmp_data, n, |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
547 par, pai, pbeta, |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
548 dummy, idummy, pvr, n, |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
549 pwork, lwork, info |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
550 F77_CHAR_ARG_LEN (1) |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
551 F77_CHAR_ARG_LEN (1))); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
552 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
553 if (info < 0) |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
554 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
555 (*current_liboctave_error_handler) ("unrecoverable error in dggev"); |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
556 return info; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
557 } |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
558 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
559 if (info > 0) |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
560 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
561 (*current_liboctave_error_handler) ("dggev failed to converge"); |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
562 return info; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
563 } |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
564 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
565 lambda.resize (n); |
9020
728e7943752d
EIG.cc: avoid volatile decl for tmp variable
John W. Eaton <jwe@octave.org>
parents:
8920
diff
changeset
|
566 octave_idx_type nvr = calc_ev ? n : 0; |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
567 v.resize (nvr, nvr); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
568 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
569 for (octave_idx_type j = 0; j < n; j++) |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
570 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
571 if (ai.elem (j) == 0.0) |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
572 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
573 lambda.elem (j) = Complex (ar.elem (j) / beta.elem (j)); |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
574 for (octave_idx_type i = 0; i < nvr; i++) |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
575 v.elem (i, j) = vr.elem (i, j); |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
576 } |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
577 else |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
578 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
579 if (j+1 >= n) |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
580 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
581 (*current_liboctave_error_handler) ("EIG: internal error"); |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
582 return -1; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
583 } |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
584 |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
585 lambda.elem(j) = Complex (ar.elem(j) / beta.elem (j), |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
586 ai.elem(j) / beta.elem (j)); |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
587 lambda.elem(j+1) = Complex (ar.elem(j+1) / beta.elem (j+1), |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
588 ai.elem(j+1) / beta.elem (j+1)); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
589 |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
590 for (octave_idx_type i = 0; i < nvr; i++) |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
591 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
592 double real_part = vr.elem (i, j); |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
593 double imag_part = vr.elem (i, j+1); |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
594 v.elem (i, j) = Complex (real_part, imag_part); |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
595 v.elem (i, j+1) = Complex (real_part, -imag_part); |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
596 } |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
597 j++; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
598 } |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
599 } |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
600 } |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
601 else |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
602 (*current_liboctave_error_handler) ("dggev workspace query failed"); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
603 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
604 return info; |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
605 } |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
606 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
607 octave_idx_type |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
608 EIG::symmetric_init (const Matrix& a, const Matrix& b, bool calc_ev) |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
609 { |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
610 octave_idx_type n = a.rows (); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
611 octave_idx_type nb = b.rows (); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
612 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
613 if (n != a.cols () || nb != b.cols ()) |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
614 { |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
615 (*current_liboctave_error_handler) ("EIG requires square matrix"); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
616 return -1; |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
617 } |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
618 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
619 if (n != nb) |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
620 { |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
621 (*current_liboctave_error_handler) ("EIG requires same size matrices"); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
622 return -1; |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
623 } |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
624 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
625 octave_idx_type info = 0; |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
626 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
627 Matrix atmp = a; |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
628 double *atmp_data = atmp.fortran_vec (); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
629 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
630 Matrix btmp = b; |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
631 double *btmp_data = btmp.fortran_vec (); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
632 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
633 ColumnVector wr (n); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
634 double *pwr = wr.fortran_vec (); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
635 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
636 octave_idx_type lwork = -1; |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
637 double dummy_work; |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
638 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
639 F77_XFCN (dsygv, DSYGV, (1, F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
640 F77_CONST_CHAR_ARG2 ("U", 1), |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
641 n, atmp_data, n, |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
642 btmp_data, n, |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
643 pwr, &dummy_work, lwork, info |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
644 F77_CHAR_ARG_LEN (1) |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
645 F77_CHAR_ARG_LEN (1))); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
646 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
647 if (info == 0) |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
648 { |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
649 lwork = static_cast<octave_idx_type> (dummy_work); |
10350
12884915a8e4
merge MArray classes & improve Array interface
Jaroslav Hajek <highegg@gmail.com>
parents:
10314
diff
changeset
|
650 Array<double> work (lwork, 1); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
651 double *pwork = work.fortran_vec (); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
652 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
653 F77_XFCN (dsygv, DSYGV, (1, F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
654 F77_CONST_CHAR_ARG2 ("U", 1), |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
655 n, atmp_data, n, |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
656 btmp_data, n, |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
657 pwr, pwork, lwork, info |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
658 F77_CHAR_ARG_LEN (1) |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
659 F77_CHAR_ARG_LEN (1))); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
660 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
661 if (info < 0) |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
662 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
663 (*current_liboctave_error_handler) ("unrecoverable error in dsygv"); |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
664 return info; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
665 } |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
666 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
667 if (info > 0) |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
668 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
669 (*current_liboctave_error_handler) ("dsygv failed to converge"); |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
670 return info; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
671 } |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
672 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
673 lambda = ComplexColumnVector (wr); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
674 v = calc_ev ? ComplexMatrix (atmp) : ComplexMatrix (); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
675 } |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
676 else |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
677 (*current_liboctave_error_handler) ("dsygv workspace query failed"); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
678 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
679 return info; |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
680 } |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
681 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
682 octave_idx_type |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
683 EIG::init (const ComplexMatrix& a, const ComplexMatrix& b, bool calc_ev) |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
684 { |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
685 if (a.any_element_is_inf_or_nan () || b.any_element_is_inf_or_nan ()) |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
686 { |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
687 (*current_liboctave_error_handler) |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
688 ("EIG: matrix contains Inf or NaN values"); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
689 return -1; |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
690 } |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
691 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
692 octave_idx_type n = a.rows (); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
693 octave_idx_type nb = b.rows (); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
694 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
695 if (n != a.cols () || nb != b.cols()) |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
696 { |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
697 (*current_liboctave_error_handler) ("EIG requires square matrix"); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
698 return -1; |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
699 } |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
700 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
701 if (n != nb) |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
702 { |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
703 (*current_liboctave_error_handler) ("EIG requires same size matrices"); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
704 return -1; |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
705 } |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
706 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
707 octave_idx_type info = 0; |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
708 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
709 ComplexMatrix tmp = b; |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
710 Complex*tmp_data = tmp.fortran_vec (); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
711 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
712 F77_XFCN (zpotrf, ZPOTRF, (F77_CONST_CHAR_ARG2 ("L", 1), |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
713 n, tmp_data, n, |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
714 info |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
715 F77_CHAR_ARG_LEN (1) |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
716 F77_CHAR_ARG_LEN (1))); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
717 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
718 if (a.is_hermitian () && b.is_hermitian () && info == 0) |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
719 return hermitian_init (a, calc_ev); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
720 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
721 ComplexMatrix atmp = a; |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
722 Complex *atmp_data = atmp.fortran_vec (); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
723 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
724 ComplexMatrix btmp = b; |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
725 Complex *btmp_data = btmp.fortran_vec (); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
726 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
727 ComplexColumnVector alpha (n); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
728 Complex *palpha = alpha.fortran_vec (); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
729 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
730 ComplexColumnVector beta (n); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
731 Complex *pbeta = beta.fortran_vec (); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
732 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
733 octave_idx_type nvr = calc_ev ? n : 0; |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
734 ComplexMatrix vtmp (nvr, nvr); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
735 Complex *pv = vtmp.fortran_vec (); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
736 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
737 octave_idx_type lwork = -1; |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
738 Complex dummy_work; |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
739 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
740 octave_idx_type lrwork = 8*n; |
10350
12884915a8e4
merge MArray classes & improve Array interface
Jaroslav Hajek <highegg@gmail.com>
parents:
10314
diff
changeset
|
741 Array<double> rwork (lrwork, 1); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
742 double *prwork = rwork.fortran_vec (); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
743 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
744 Complex *dummy = 0; |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
745 octave_idx_type idummy = 1; |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
746 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
747 F77_XFCN (zggev, ZGGEV, (F77_CONST_CHAR_ARG2 ("N", 1), |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
748 F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
749 n, atmp_data, n, btmp_data, n, |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
750 palpha, pbeta, dummy, idummy, |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
751 pv, n, &dummy_work, lwork, prwork, info |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
752 F77_CHAR_ARG_LEN (1) |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
753 F77_CHAR_ARG_LEN (1))); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
754 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
755 if (info == 0) |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
756 { |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
757 lwork = static_cast<octave_idx_type> (dummy_work.real ()); |
10350
12884915a8e4
merge MArray classes & improve Array interface
Jaroslav Hajek <highegg@gmail.com>
parents:
10314
diff
changeset
|
758 Array<Complex> work (lwork, 1); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
759 Complex *pwork = work.fortran_vec (); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
760 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
761 F77_XFCN (zggev, ZGGEV, (F77_CONST_CHAR_ARG2 ("N", 1), |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
762 F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
763 n, atmp_data, n, btmp_data, n, |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
764 palpha, pbeta, dummy, idummy, |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
765 pv, n, pwork, lwork, prwork, info |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
766 F77_CHAR_ARG_LEN (1) |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
767 F77_CHAR_ARG_LEN (1))); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
768 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
769 if (info < 0) |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
770 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
771 (*current_liboctave_error_handler) ("unrecoverable error in zggev"); |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
772 return info; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
773 } |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
774 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
775 if (info > 0) |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
776 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
777 (*current_liboctave_error_handler) ("zggev failed to converge"); |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
778 return info; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
779 } |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
780 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
781 lambda.resize (n); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
782 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
783 for (octave_idx_type j = 0; j < n; j++) |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
784 lambda.elem (j) = alpha.elem (j) / beta.elem(j); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
785 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
786 v = vtmp; |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
787 } |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
788 else |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
789 (*current_liboctave_error_handler) ("zggev workspace query failed"); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
790 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
791 return info; |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
792 } |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
793 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
794 octave_idx_type |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
795 EIG::hermitian_init (const ComplexMatrix& a, const ComplexMatrix& b, bool calc_ev) |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
796 { |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
797 octave_idx_type n = a.rows (); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
798 octave_idx_type nb = b.rows (); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
799 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
800 if (n != a.cols () || nb != b.cols ()) |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
801 { |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
802 (*current_liboctave_error_handler) ("EIG requires square matrix"); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
803 return -1; |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
804 } |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
805 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
806 if (n != nb) |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
807 { |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
808 (*current_liboctave_error_handler) ("EIG requires same size matrices"); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
809 return -1; |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
810 } |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
811 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
812 octave_idx_type info = 0; |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
813 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
814 ComplexMatrix atmp = a; |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
815 Complex *atmp_data = atmp.fortran_vec (); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
816 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
817 ComplexMatrix btmp = b; |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
818 Complex *btmp_data = btmp.fortran_vec (); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
819 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
820 ColumnVector wr (n); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
821 double *pwr = wr.fortran_vec (); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
822 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
823 octave_idx_type lwork = -1; |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
824 Complex dummy_work; |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
825 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
826 octave_idx_type lrwork = 3*n; |
10350
12884915a8e4
merge MArray classes & improve Array interface
Jaroslav Hajek <highegg@gmail.com>
parents:
10314
diff
changeset
|
827 Array<double> rwork (lrwork, 1); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
828 double *prwork = rwork.fortran_vec (); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
829 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
830 F77_XFCN (zhegv, ZHEGV, (1, F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
831 F77_CONST_CHAR_ARG2 ("U", 1), |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
832 n, atmp_data, n, |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
833 btmp_data, n, |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
834 pwr, &dummy_work, lwork, |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
835 prwork, info |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
836 F77_CHAR_ARG_LEN (1) |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
837 F77_CHAR_ARG_LEN (1))); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
838 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
839 if (info == 0) |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
840 { |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
841 lwork = static_cast<octave_idx_type> (dummy_work.real ()); |
10350
12884915a8e4
merge MArray classes & improve Array interface
Jaroslav Hajek <highegg@gmail.com>
parents:
10314
diff
changeset
|
842 Array<Complex> work (lwork, 1); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
843 Complex *pwork = work.fortran_vec (); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
844 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
845 F77_XFCN (zhegv, ZHEGV, (1, F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
846 F77_CONST_CHAR_ARG2 ("U", 1), |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
847 n, atmp_data, n, |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
848 btmp_data, n, |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
849 pwr, pwork, lwork, prwork, info |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
850 F77_CHAR_ARG_LEN (1) |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
851 F77_CHAR_ARG_LEN (1))); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
852 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
853 if (info < 0) |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
854 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
855 (*current_liboctave_error_handler) ("unrecoverable error in zhegv"); |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
856 return info; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
857 } |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
858 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
859 if (info > 0) |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
860 { |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
861 (*current_liboctave_error_handler) ("zhegv failed to converge"); |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
862 return info; |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
863 } |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
864 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
865 lambda = ComplexColumnVector (wr); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
866 v = calc_ev ? ComplexMatrix (atmp) : ComplexMatrix (); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
867 } |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
868 else |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
869 (*current_liboctave_error_handler) ("zhegv workspace query failed"); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
870 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
871 return info; |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
872 } |