Mercurial > octave-nkf
annotate liboctave/EIG.cc @ 7961:a5d1e27ee1f4 ss-3-1-51
3.1.51 snapshot
author | John W. Eaton <jwe@octave.org> |
---|---|
date | Tue, 22 Jul 2008 11:40:48 -0400 |
parents | 29980c6b8604 |
children | 18c4ded8612a |
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); | |
462 | 68 } |
69 | |
5275 | 70 octave_idx_type |
4725 | 71 EIG::init (const Matrix& a, bool calc_ev) |
462 | 72 { |
5822 | 73 if (a.any_element_is_inf_or_nan ()) |
74 { | |
75 (*current_liboctave_error_handler) | |
76 ("EIG: matrix contains Inf or NaN values"); | |
77 return -1; | |
78 } | |
79 | |
2815 | 80 if (a.is_symmetric ()) |
4725 | 81 return symmetric_init (a, calc_ev); |
2815 | 82 |
5275 | 83 octave_idx_type n = a.rows (); |
1934 | 84 |
85 if (n != a.cols ()) | |
462 | 86 { |
87 (*current_liboctave_error_handler) ("EIG requires square matrix"); | |
88 return -1; | |
89 } | |
90 | |
5275 | 91 octave_idx_type info = 0; |
462 | 92 |
1934 | 93 Matrix atmp = a; |
94 double *tmp_data = atmp.fortran_vec (); | |
95 | |
96 Array<double> wr (n); | |
97 double *pwr = wr.fortran_vec (); | |
98 | |
99 Array<double> wi (n); | |
100 double *pwi = wi.fortran_vec (); | |
101 | |
5275 | 102 volatile octave_idx_type nvr = calc_ev ? n : 0; |
4725 | 103 Matrix vr (nvr, nvr); |
462 | 104 double *pvr = vr.fortran_vec (); |
1934 | 105 |
5275 | 106 octave_idx_type lwork = -1; |
4800 | 107 double dummy_work; |
462 | 108 |
1365 | 109 double *dummy = 0; |
5275 | 110 octave_idx_type idummy = 1; |
462 | 111 |
4552 | 112 F77_XFCN (dgeev, DGEEV, (F77_CONST_CHAR_ARG2 ("N", 1), |
4725 | 113 F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), |
4552 | 114 n, tmp_data, n, pwr, pwi, dummy, |
4800 | 115 idummy, pvr, n, &dummy_work, lwork, info |
4552 | 116 F77_CHAR_ARG_LEN (1) |
117 F77_CHAR_ARG_LEN (1))); | |
462 | 118 |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
119 if (info == 0) |
462 | 120 { |
5275 | 121 lwork = static_cast<octave_idx_type> (dummy_work); |
4800 | 122 Array<double> work (lwork); |
123 double *pwork = work.fortran_vec (); | |
124 | |
125 F77_XFCN (dgeev, DGEEV, (F77_CONST_CHAR_ARG2 ("N", 1), | |
126 F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), | |
127 n, tmp_data, n, pwr, pwi, dummy, | |
128 idummy, pvr, n, pwork, lwork, info | |
129 F77_CHAR_ARG_LEN (1) | |
130 F77_CHAR_ARG_LEN (1))); | |
131 | |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
132 if (info < 0) |
4800 | 133 { |
134 (*current_liboctave_error_handler) ("unrecoverable error in dgeev"); | |
135 return info; | |
136 } | |
137 | |
2815 | 138 if (info > 0) |
139 { | |
4800 | 140 (*current_liboctave_error_handler) ("dgeev failed to converge"); |
141 return info; | |
142 } | |
1934 | 143 |
4800 | 144 lambda.resize (n); |
145 v.resize (nvr, nvr); | |
146 | |
5275 | 147 for (octave_idx_type j = 0; j < n; j++) |
4800 | 148 { |
149 if (wi.elem (j) == 0.0) | |
462 | 150 { |
4800 | 151 lambda.elem (j) = Complex (wr.elem (j)); |
5275 | 152 for (octave_idx_type i = 0; i < nvr; i++) |
4800 | 153 v.elem (i, j) = vr.elem (i, j); |
154 } | |
155 else | |
156 { | |
157 if (j+1 >= n) | |
2815 | 158 { |
4800 | 159 (*current_liboctave_error_handler) ("EIG: internal error"); |
160 return -1; | |
2815 | 161 } |
4800 | 162 |
163 lambda.elem(j) = Complex (wr.elem(j), wi.elem(j)); | |
164 lambda.elem(j+1) = Complex (wr.elem(j+1), wi.elem(j+1)); | |
2815 | 165 |
5275 | 166 for (octave_idx_type i = 0; i < nvr; i++) |
4800 | 167 { |
168 double real_part = vr.elem (i, j); | |
169 double imag_part = vr.elem (i, j+1); | |
170 v.elem (i, j) = Complex (real_part, imag_part); | |
171 v.elem (i, j+1) = Complex (real_part, -imag_part); | |
1934 | 172 } |
4800 | 173 j++; |
2815 | 174 } |
175 } | |
176 } | |
4800 | 177 else |
178 (*current_liboctave_error_handler) ("dgeev workspace query failed"); | |
2815 | 179 |
180 return info; | |
181 } | |
182 | |
5275 | 183 octave_idx_type |
4725 | 184 EIG::symmetric_init (const Matrix& a, bool calc_ev) |
2815 | 185 { |
5275 | 186 octave_idx_type n = a.rows (); |
1934 | 187 |
2815 | 188 if (n != a.cols ()) |
189 { | |
190 (*current_liboctave_error_handler) ("EIG requires square matrix"); | |
191 return -1; | |
192 } | |
193 | |
5275 | 194 octave_idx_type info = 0; |
2815 | 195 |
196 Matrix atmp = a; | |
197 double *tmp_data = atmp.fortran_vec (); | |
198 | |
3585 | 199 ColumnVector wr (n); |
2815 | 200 double *pwr = wr.fortran_vec (); |
201 | |
5275 | 202 octave_idx_type lwork = -1; |
4800 | 203 double dummy_work; |
2815 | 204 |
4725 | 205 F77_XFCN (dsyev, DSYEV, (F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), |
4552 | 206 F77_CONST_CHAR_ARG2 ("U", 1), |
4800 | 207 n, tmp_data, n, pwr, &dummy_work, lwork, info |
4552 | 208 F77_CHAR_ARG_LEN (1) |
209 F77_CHAR_ARG_LEN (1))); | |
2815 | 210 |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
211 if (info == 0) |
2815 | 212 { |
5275 | 213 lwork = static_cast<octave_idx_type> (dummy_work); |
4800 | 214 Array<double> work (lwork); |
215 double *pwork = work.fortran_vec (); | |
216 | |
217 F77_XFCN (dsyev, DSYEV, (F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), | |
218 F77_CONST_CHAR_ARG2 ("U", 1), | |
219 n, tmp_data, n, pwr, pwork, lwork, info | |
220 F77_CHAR_ARG_LEN (1) | |
221 F77_CHAR_ARG_LEN (1))); | |
222 | |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
223 if (info < 0) |
4800 | 224 { |
225 (*current_liboctave_error_handler) ("unrecoverable error in dsyev"); | |
226 return info; | |
227 } | |
228 | |
229 if (info > 0) | |
230 { | |
231 (*current_liboctave_error_handler) ("dsyev failed to converge"); | |
232 return info; | |
233 } | |
234 | |
3585 | 235 lambda = ComplexColumnVector (wr); |
4725 | 236 v = calc_ev ? ComplexMatrix (atmp) : ComplexMatrix (); |
462 | 237 } |
4800 | 238 else |
239 (*current_liboctave_error_handler) ("dsyev workspace query failed"); | |
462 | 240 |
241 return info; | |
242 } | |
243 | |
5275 | 244 octave_idx_type |
4725 | 245 EIG::init (const ComplexMatrix& a, bool calc_ev) |
462 | 246 { |
5822 | 247 if (a.any_element_is_inf_or_nan ()) |
248 { | |
249 (*current_liboctave_error_handler) | |
250 ("EIG: matrix contains Inf or NaN values"); | |
251 return -1; | |
252 } | |
253 | |
2815 | 254 if (a.is_hermitian ()) |
4725 | 255 return hermitian_init (a, calc_ev); |
2815 | 256 |
5275 | 257 octave_idx_type n = a.rows (); |
1934 | 258 |
259 if (n != a.cols ()) | |
462 | 260 { |
261 (*current_liboctave_error_handler) ("EIG requires square matrix"); | |
262 return -1; | |
263 } | |
264 | |
5275 | 265 octave_idx_type info = 0; |
462 | 266 |
1934 | 267 ComplexMatrix atmp = a; |
268 Complex *tmp_data = atmp.fortran_vec (); | |
462 | 269 |
2815 | 270 ComplexColumnVector w (n); |
271 Complex *pw = w.fortran_vec (); | |
272 | |
5275 | 273 octave_idx_type nvr = calc_ev ? n : 0; |
4725 | 274 ComplexMatrix vtmp (nvr, nvr); |
2815 | 275 Complex *pv = vtmp.fortran_vec (); |
276 | |
5275 | 277 octave_idx_type lwork = -1; |
4800 | 278 Complex dummy_work; |
1934 | 279 |
5275 | 280 octave_idx_type lrwork = 2*n; |
2815 | 281 Array<double> rwork (lrwork); |
1934 | 282 double *prwork = rwork.fortran_vec (); |
462 | 283 |
1365 | 284 Complex *dummy = 0; |
5275 | 285 octave_idx_type idummy = 1; |
462 | 286 |
4552 | 287 F77_XFCN (zgeev, ZGEEV, (F77_CONST_CHAR_ARG2 ("N", 1), |
4725 | 288 F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), |
4552 | 289 n, tmp_data, n, pw, dummy, idummy, |
4800 | 290 pv, n, &dummy_work, lwork, prwork, info |
4552 | 291 F77_CHAR_ARG_LEN (1) |
292 F77_CHAR_ARG_LEN (1))); | |
2815 | 293 |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
294 if (info == 0) |
2815 | 295 { |
5275 | 296 lwork = static_cast<octave_idx_type> (dummy_work.real ()); |
4800 | 297 Array<Complex> work (lwork); |
298 Complex *pwork = work.fortran_vec (); | |
299 | |
300 F77_XFCN (zgeev, ZGEEV, (F77_CONST_CHAR_ARG2 ("N", 1), | |
301 F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), | |
302 n, tmp_data, n, pw, dummy, idummy, | |
303 pv, n, pwork, lwork, prwork, info | |
304 F77_CHAR_ARG_LEN (1) | |
305 F77_CHAR_ARG_LEN (1))); | |
306 | |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
307 if (info < 0) |
4800 | 308 { |
309 (*current_liboctave_error_handler) ("unrecoverable error in zgeev"); | |
310 return info; | |
311 } | |
312 | |
313 if (info > 0) | |
314 { | |
315 (*current_liboctave_error_handler) ("zgeev failed to converge"); | |
316 return info; | |
317 } | |
318 | |
2815 | 319 lambda = w; |
320 v = vtmp; | |
321 } | |
4800 | 322 else |
323 (*current_liboctave_error_handler) ("zgeev workspace query failed"); | |
2815 | 324 |
325 return info; | |
326 } | |
327 | |
5275 | 328 octave_idx_type |
4725 | 329 EIG::hermitian_init (const ComplexMatrix& a, bool calc_ev) |
2815 | 330 { |
5275 | 331 octave_idx_type n = a.rows (); |
2815 | 332 |
333 if (n != a.cols ()) | |
334 { | |
335 (*current_liboctave_error_handler) ("EIG requires square matrix"); | |
336 return -1; | |
337 } | |
338 | |
5275 | 339 octave_idx_type info = 0; |
462 | 340 |
2815 | 341 ComplexMatrix atmp = a; |
342 Complex *tmp_data = atmp.fortran_vec (); | |
343 | |
3585 | 344 ColumnVector wr (n); |
345 double *pwr = wr.fortran_vec (); | |
2815 | 346 |
5275 | 347 octave_idx_type lwork = -1; |
4800 | 348 Complex dummy_work; |
2815 | 349 |
5275 | 350 octave_idx_type lrwork = 3*n; |
2815 | 351 Array<double> rwork (lrwork); |
352 double *prwork = rwork.fortran_vec (); | |
353 | |
4725 | 354 F77_XFCN (zheev, ZHEEV, (F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), |
4552 | 355 F77_CONST_CHAR_ARG2 ("U", 1), |
4800 | 356 n, tmp_data, n, pwr, &dummy_work, lwork, |
357 prwork, info | |
4552 | 358 F77_CHAR_ARG_LEN (1) |
359 F77_CHAR_ARG_LEN (1))); | |
2815 | 360 |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
361 if (info == 0) |
2815 | 362 { |
5275 | 363 lwork = static_cast<octave_idx_type> (dummy_work.real ()); |
4800 | 364 Array<Complex> work (lwork); |
365 Complex *pwork = work.fortran_vec (); | |
366 | |
367 F77_XFCN (zheev, ZHEEV, (F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), | |
368 F77_CONST_CHAR_ARG2 ("U", 1), | |
369 n, tmp_data, n, pwr, pwork, lwork, prwork, info | |
370 F77_CHAR_ARG_LEN (1) | |
371 F77_CHAR_ARG_LEN (1))); | |
372 | |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
373 if (info < 0) |
4800 | 374 { |
375 (*current_liboctave_error_handler) ("unrecoverable error in zheev"); | |
376 return info; | |
377 } | |
378 | |
379 if (info > 0) | |
380 { | |
381 (*current_liboctave_error_handler) ("zheev failed to converge"); | |
382 return info; | |
383 } | |
384 | |
3585 | 385 lambda = ComplexColumnVector (wr); |
4725 | 386 v = calc_ev ? ComplexMatrix (atmp) : ComplexMatrix (); |
2815 | 387 } |
4800 | 388 else |
389 (*current_liboctave_error_handler) ("zheev workspace query failed"); | |
462 | 390 |
391 return info; | |
392 } | |
393 | |
394 /* | |
395 ;;; Local Variables: *** | |
396 ;;; mode: C++ *** | |
397 ;;; End: *** | |
398 */ |