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