GNU Octave 10.1.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 
Loading...
Searching...
No Matches
fEIG.cc
Go to the documentation of this file.
1////////////////////////////////////////////////////////////////////////
2//
3// Copyright (C) 1994-2025 The Octave Project Developers
4//
5// See the file COPYRIGHT.md in the top-level directory of this
6// distribution or <https://octave.org/copyright/>.
7//
8// This file is part of Octave.
9//
10// Octave is free software: you can redistribute it and/or modify it
11// under the terms of the GNU General Public License as published by
12// the Free Software Foundation, either version 3 of the License, or
13// (at your option) any later version.
14//
15// Octave is distributed in the hope that it will be useful, but
16// WITHOUT ANY WARRANTY; without even the implied warranty of
17// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18// GNU General Public License for more details.
19//
20// You should have received a copy of the GNU General Public License
21// along with Octave; see the file COPYING. If not, see
22// <https://www.gnu.org/licenses/>.
23//
24////////////////////////////////////////////////////////////////////////
25
26#if defined (HAVE_CONFIG_H)
27# include "config.h"
28#endif
29
30#include "Array.h"
31#include "fEIG.h"
32#include "fColVector.h"
33#include "fMatrix.h"
34#include "lo-error.h"
35#include "lo-lapack-proto.h"
36
37octave_idx_type
38FloatEIG::init (const FloatMatrix& a, bool calc_rev, bool calc_lev,
39 bool balance)
40{
41 if (a.any_element_is_inf_or_nan ())
42 (*current_liboctave_error_handler)
43 ("EIG: matrix contains Inf or NaN values");
44
45 if (a.issymmetric ())
46 return symmetric_init (a, calc_rev, calc_lev);
47
48 F77_INT n = octave::to_f77_int (a.rows ());
49 F77_INT a_nc = octave::to_f77_int (a.cols ());
50
51 if (n != a_nc)
52 (*current_liboctave_error_handler) ("EIG requires square matrix");
53
54 F77_INT info = 0;
55
56 FloatMatrix atmp = a;
57 float *tmp_data = atmp.rwdata ();
58
59 Array<float> wr (dim_vector (n, 1));
60 float *pwr = wr.rwdata ();
61
62 Array<float> wi (dim_vector (n, 1));
63 float *pwi = wi.rwdata ();
64
65 F77_INT nvr = (calc_rev ? n : 0);
66 FloatMatrix vr (nvr, nvr);
67 float *pvr = vr.rwdata ();
68
69 F77_INT nvl = (calc_lev ? n : 0);
70 FloatMatrix vl (nvl, nvl);
71 float *pvl = vl.rwdata ();
72
73 F77_INT lwork = -1;
74 float dummy_work;
75
76 F77_INT ilo;
77 F77_INT ihi;
78
79 Array<float> scale (dim_vector (n, 1));
80 float *pscale = scale.rwdata ();
81
82 float abnrm;
83
84 Array<float> rconde (dim_vector (n, 1));
85 float *prconde = rconde.rwdata ();
86
87 Array<float> rcondv (dim_vector (n, 1));
88 float *prcondv = rcondv.rwdata ();
89
90 F77_INT dummy_iwork;
91
92 F77_XFCN (sgeevx, SGEEVX, (F77_CONST_CHAR_ARG2 (balance ? "B" : "N", 1),
93 F77_CONST_CHAR_ARG2 ("N", 1),
94 F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1),
95 F77_CONST_CHAR_ARG2 ("N", 1),
96 n, tmp_data, n, pwr, pwi,
97 pvl, n, pvr, n,
98 ilo, ihi, pscale, abnrm, prconde, prcondv,
99 &dummy_work, lwork, &dummy_iwork, info
100 F77_CHAR_ARG_LEN (1)
101 F77_CHAR_ARG_LEN (1)
102 F77_CHAR_ARG_LEN (1)
103 F77_CHAR_ARG_LEN (1)));
104
105 if (info != 0)
106 (*current_liboctave_error_handler) ("sgeevx workspace query failed");
107
108 lwork = static_cast<F77_INT> (dummy_work);
109 Array<float> work (dim_vector (lwork, 1));
110 float *pwork = work.rwdata ();
111
112 F77_XFCN (sgeevx, SGEEVX, (F77_CONST_CHAR_ARG2 (balance ? "B" : "N", 1),
113 F77_CONST_CHAR_ARG2 ("N", 1),
114 F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1),
115 F77_CONST_CHAR_ARG2 ("N", 1),
116 n, tmp_data, n, pwr, pwi,
117 pvl, n, pvr, n,
118 ilo, ihi, pscale, abnrm, prconde, prcondv,
119 pwork, lwork, &dummy_iwork, info
120 F77_CHAR_ARG_LEN (1)
121 F77_CHAR_ARG_LEN (1)
122 F77_CHAR_ARG_LEN (1)
123 F77_CHAR_ARG_LEN (1)));
124
125 if (info < 0)
126 (*current_liboctave_error_handler) ("unrecoverable error in sgeevx");
127
128 if (info > 0)
129 (*current_liboctave_error_handler) ("sgeevx failed to converge");
130
131 m_lambda.resize (n);
132 m_v.resize (nvr, nvr);
133 m_w.resize (nvl, nvl);
134
135 for (F77_INT j = 0; j < n; j++)
136 {
137 if (wi.elem (j) == 0.0)
138 {
139 m_lambda.elem (j) = FloatComplex (wr.elem (j));
140 for (octave_idx_type i = 0; i < nvr; i++)
141 m_v.elem (i, j) = vr.elem (i, j);
142
143 for (F77_INT i = 0; i < nvl; i++)
144 m_w.elem (i, j) = vl.elem (i, j);
145 }
146 else
147 {
148 if (j+1 >= n)
149 (*current_liboctave_error_handler) ("EIG: internal error");
150
151 m_lambda.elem (j) = FloatComplex (wr.elem (j), wi.elem (j));
152 m_lambda.elem (j+1) = FloatComplex (wr.elem (j+1), wi.elem (j+1));
153
154 for (F77_INT i = 0; i < nvr; i++)
155 {
156 float real_part = vr.elem (i, j);
157 float imag_part = vr.elem (i, j+1);
158 m_v.elem (i, j) = FloatComplex (real_part, imag_part);
159 m_v.elem (i, j+1) = FloatComplex (real_part, -imag_part);
160 }
161 for (F77_INT i = 0; i < nvl; i++)
162 {
163 float real_part = vl.elem (i, j);
164 float imag_part = vl.elem (i, j+1);
165 m_w.elem (i, j) = FloatComplex (real_part, imag_part);
166 m_w.elem (i, j+1) = FloatComplex (real_part, -imag_part);
167 }
168 j++;
169 }
170 }
171
172 return info;
173}
174
175octave_idx_type
176FloatEIG::symmetric_init (const FloatMatrix& a, bool calc_rev, bool calc_lev)
177{
178 F77_INT n = octave::to_f77_int (a.rows ());
179 F77_INT a_nc = octave::to_f77_int (a.cols ());
180
181 if (n != a_nc)
182 (*current_liboctave_error_handler) ("EIG requires square matrix");
183
184 F77_INT info = 0;
185
186 FloatMatrix atmp = a;
187 float *tmp_data = atmp.rwdata ();
188
189 FloatColumnVector wr (n);
190 float *pwr = wr.rwdata ();
191
192 F77_INT lwork = -1;
193 float dummy_work;
194
195 F77_XFCN (ssyev, SSYEV, (F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1),
196 F77_CONST_CHAR_ARG2 ("U", 1),
197 n, tmp_data, n, pwr, &dummy_work, lwork, info
198 F77_CHAR_ARG_LEN (1)
199 F77_CHAR_ARG_LEN (1)));
200
201 if (info != 0)
202 (*current_liboctave_error_handler) ("ssyev workspace query failed");
203
204 lwork = static_cast<F77_INT> (dummy_work);
205 Array<float> work (dim_vector (lwork, 1));
206 float *pwork = work.rwdata ();
207
208 F77_XFCN (ssyev, SSYEV, (F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1),
209 F77_CONST_CHAR_ARG2 ("U", 1),
210 n, tmp_data, n, pwr, pwork, lwork, info
211 F77_CHAR_ARG_LEN (1)
212 F77_CHAR_ARG_LEN (1)));
213
214 if (info < 0)
215 (*current_liboctave_error_handler) ("unrecoverable error in ssyev");
216
217 if (info > 0)
218 (*current_liboctave_error_handler) ("ssyev failed to converge");
219
220 m_lambda = FloatComplexColumnVector (wr);
221 m_v = (calc_rev ? FloatComplexMatrix (atmp) : FloatComplexMatrix ());
222 m_w = (calc_lev ? FloatComplexMatrix (atmp) : FloatComplexMatrix ());
223
224 return info;
225}
226
227octave_idx_type
228FloatEIG::init (const FloatComplexMatrix& a, bool calc_rev, bool calc_lev,
229 bool balance)
230{
231 if (a.any_element_is_inf_or_nan ())
232 (*current_liboctave_error_handler)
233 ("EIG: matrix contains Inf or NaN values");
234
235 if (a.ishermitian ())
236 return hermitian_init (a, calc_rev, calc_lev);
237
238 F77_INT n = octave::to_f77_int (a.rows ());
239 F77_INT a_nc = octave::to_f77_int (a.cols ());
240
241 if (n != a_nc)
242 (*current_liboctave_error_handler) ("EIG requires square matrix");
243
244 F77_INT info = 0;
245
246 FloatComplexMatrix atmp = a;
247 FloatComplex *tmp_data = atmp.rwdata ();
248
249 FloatComplexColumnVector wr (n);
250 FloatComplex *pw = wr.rwdata ();
251
252 F77_INT nvr = (calc_rev ? n : 0);
253 FloatComplexMatrix vrtmp (nvr, nvr);
254 FloatComplex *pvr = vrtmp.rwdata ();
255
256 F77_INT nvl = (calc_lev ? n : 0);
257 FloatComplexMatrix vltmp (nvl, nvl);
258 FloatComplex *pvl = vltmp.rwdata ();
259
260 F77_INT lwork = -1;
261 FloatComplex dummy_work;
262
263 F77_INT lrwork = 2*n;
264 Array<float> rwork (dim_vector (lrwork, 1));
265 float *prwork = rwork.rwdata ();
266
267 F77_INT ilo;
268 F77_INT ihi;
269
270 Array<float> scale (dim_vector (n, 1));
271 float *pscale = scale.rwdata ();
272
273 float abnrm;
274
275 Array<float> rconde (dim_vector (n, 1));
276 float *prconde = rconde.rwdata ();
277
278 Array<float> rcondv (dim_vector (n, 1));
279 float *prcondv = rcondv.rwdata ();
280
281 F77_XFCN (cgeevx, CGEEVX, (F77_CONST_CHAR_ARG2 (balance ? "B" : "N", 1),
282 F77_CONST_CHAR_ARG2 (calc_lev ? "V" : "N", 1),
283 F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1),
284 F77_CONST_CHAR_ARG2 ("N", 1),
285 n, F77_CMPLX_ARG (tmp_data), n, F77_CMPLX_ARG (pw),
286 F77_CMPLX_ARG (pvl), n, F77_CMPLX_ARG (pvr), n,
287 ilo, ihi, pscale, abnrm, prconde, prcondv,
288 F77_CMPLX_ARG (&dummy_work), lwork, prwork, info
289 F77_CHAR_ARG_LEN (1)
290 F77_CHAR_ARG_LEN (1)
291 F77_CHAR_ARG_LEN (1)
292 F77_CHAR_ARG_LEN (1)));
293
294 if (info != 0)
295 (*current_liboctave_error_handler) ("cgeevx workspace query failed");
296
297 lwork = static_cast<F77_INT> (dummy_work.real ());
298 Array<FloatComplex> work (dim_vector (lwork, 1));
299 FloatComplex *pwork = work.rwdata ();
300
301 F77_XFCN (cgeevx, CGEEVX, (F77_CONST_CHAR_ARG2 (balance ? "B" : "N", 1),
302 F77_CONST_CHAR_ARG2 (calc_lev ? "V" : "N", 1),
303 F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1),
304 F77_CONST_CHAR_ARG2 ("N", 1),
305 n, F77_CMPLX_ARG (tmp_data), n, F77_CMPLX_ARG (pw),
306 F77_CMPLX_ARG (pvl), n, F77_CMPLX_ARG (pvr), n,
307 ilo, ihi, pscale, abnrm, prconde, prcondv,
308 F77_CMPLX_ARG (pwork), lwork, prwork, info
309 F77_CHAR_ARG_LEN (1)
310 F77_CHAR_ARG_LEN (1)
311 F77_CHAR_ARG_LEN (1)
312 F77_CHAR_ARG_LEN (1)));
313
314 if (info < 0)
315 (*current_liboctave_error_handler) ("unrecoverable error in cgeevx");
316
317 if (info > 0)
318 (*current_liboctave_error_handler) ("cgeevx failed to converge");
319
320 m_lambda = wr;
321 m_v = vrtmp;
322 m_w = vltmp;
323
324 return info;
325}
326
327octave_idx_type
328FloatEIG::hermitian_init (const FloatComplexMatrix& a, bool calc_rev,
329 bool calc_lev)
330{
331 F77_INT n = octave::to_f77_int (a.rows ());
332 F77_INT a_nc = octave::to_f77_int (a.cols ());
333
334 if (n != a_nc)
335 (*current_liboctave_error_handler) ("EIG requires square matrix");
336
337 F77_INT info = 0;
338
339 FloatComplexMatrix atmp = a;
340 FloatComplex *tmp_data = atmp.rwdata ();
341
342 FloatColumnVector wr (n);
343 float *pwr = wr.rwdata ();
344
345 F77_INT lwork = -1;
346 FloatComplex dummy_work;
347
348 F77_INT lrwork = 3*n;
349 Array<float> rwork (dim_vector (lrwork, 1));
350 float *prwork = rwork.rwdata ();
351
352 F77_XFCN (cheev, CHEEV, (F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1),
353 F77_CONST_CHAR_ARG2 ("U", 1),
354 n, F77_CMPLX_ARG (tmp_data), n, pwr,
355 F77_CMPLX_ARG (&dummy_work), lwork,
356 prwork, info
357 F77_CHAR_ARG_LEN (1)
358 F77_CHAR_ARG_LEN (1)));
359
360 if (info != 0)
361 (*current_liboctave_error_handler) ("cheev workspace query failed");
362
363 lwork = static_cast<F77_INT> (dummy_work.real ());
364 Array<FloatComplex> work (dim_vector (lwork, 1));
365 FloatComplex *pwork = work.rwdata ();
366
367 F77_XFCN (cheev, CHEEV, (F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1),
368 F77_CONST_CHAR_ARG2 ("U", 1),
369 n, F77_CMPLX_ARG (tmp_data), n, pwr,
370 F77_CMPLX_ARG (pwork), lwork, prwork, info
371 F77_CHAR_ARG_LEN (1)
372 F77_CHAR_ARG_LEN (1)));
373
374 if (info < 0)
375 (*current_liboctave_error_handler) ("unrecoverable error in cheev");
376
377 if (info > 0)
378 (*current_liboctave_error_handler) ("cheev failed to converge");
379
380 m_lambda = FloatComplexColumnVector (wr);
381 m_v = (calc_rev ? FloatComplexMatrix (atmp) : FloatComplexMatrix ());
382 m_w = (calc_lev ? FloatComplexMatrix (atmp) : FloatComplexMatrix ());
383
384 return info;
385}
386
387octave_idx_type
388FloatEIG::init (const FloatMatrix& a, const FloatMatrix& b, bool calc_rev,
389 bool calc_lev, bool force_qz)
390{
391 if (a.any_element_is_inf_or_nan () || b.any_element_is_inf_or_nan ())
392 (*current_liboctave_error_handler)
393 ("EIG: matrix contains Inf or NaN values");
394
395 F77_INT n = octave::to_f77_int (a.rows ());
396 F77_INT nb = octave::to_f77_int (b.rows ());
397
398 F77_INT a_nc = octave::to_f77_int (a.cols ());
399 F77_INT b_nc = octave::to_f77_int (b.cols ());
400
401 if (n != a_nc || nb != b_nc)
402 (*current_liboctave_error_handler) ("EIG requires square matrix");
403
404 if (n != nb)
405 (*current_liboctave_error_handler) ("EIG requires same size matrices");
406
407 F77_INT info = 0;
408
409 FloatMatrix tmp = b;
410 float *tmp_data = tmp.rwdata ();
411 if (! force_qz)
412 {
413 F77_XFCN (spotrf, SPOTRF, (F77_CONST_CHAR_ARG2 ("L", 1),
414 n, tmp_data, n,
415 info
416 F77_CHAR_ARG_LEN (1)));
417
418 if (a.issymmetric () && b.issymmetric () && info == 0)
419 return symmetric_init (a, b, calc_rev, calc_lev);
420 }
421
422 FloatMatrix atmp = a;
423 float *atmp_data = atmp.rwdata ();
424
425 FloatMatrix btmp = b;
426 float *btmp_data = btmp.rwdata ();
427
428 Array<float> ar (dim_vector (n, 1));
429 float *par = ar.rwdata ();
430
431 Array<float> ai (dim_vector (n, 1));
432 float *pai = ai.rwdata ();
433
434 Array<float> beta (dim_vector (n, 1));
435 float *pbeta = beta.rwdata ();
436
437 F77_INT nvr = (calc_rev ? n : 0);
438 FloatMatrix vr (nvr, nvr);
439 float *pvr = vr.rwdata ();
440
441 F77_INT nvl = (calc_lev ? n : 0);
442 FloatMatrix vl (nvl, nvl);
443 float *pvl = vl.rwdata ();
444
445 F77_INT lwork = -1;
446 float dummy_work;
447
448 F77_XFCN (sggev, SGGEV, (F77_CONST_CHAR_ARG2 (calc_lev ? "V" : "N", 1),
449 F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1),
450 n, atmp_data, n, btmp_data, n,
451 par, pai, pbeta,
452 pvl, n, pvr, n,
453 &dummy_work, lwork, info
454 F77_CHAR_ARG_LEN (1)
455 F77_CHAR_ARG_LEN (1)));
456
457 if (info != 0)
458 (*current_liboctave_error_handler) ("sggev workspace query failed");
459
460 lwork = static_cast<F77_INT> (dummy_work);
461 Array<float> work (dim_vector (lwork, 1));
462 float *pwork = work.rwdata ();
463
464 F77_XFCN (sggev, SGGEV, (F77_CONST_CHAR_ARG2 (calc_lev ? "V" : "N", 1),
465 F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1),
466 n, atmp_data, n, btmp_data, n,
467 par, pai, pbeta,
468 pvl, n, pvr, n,
469 pwork, lwork, info
470 F77_CHAR_ARG_LEN (1)
471 F77_CHAR_ARG_LEN (1)));
472
473 if (info < 0)
474 (*current_liboctave_error_handler) ("unrecoverable error in sggev");
475
476 if (info > 0)
477 (*current_liboctave_error_handler) ("sggev failed to converge");
478
479 m_lambda.resize (n);
480 m_v.resize (nvr, nvr);
481 m_w.resize (nvl, nvl);
482
483
484 for (F77_INT j = 0; j < n; j++)
485 {
486 if (ai.elem (j) == 0.0)
487 {
488 m_lambda.elem (j) = FloatComplex (ar.elem (j) / beta.elem (j));
489 for (F77_INT i = 0; i < nvr; i++)
490 m_v.elem (i, j) = vr.elem (i, j);
491
492 for (F77_INT i = 0; i < nvl; i++)
493 m_w.elem (i, j) = vl.elem (i, j);
494 }
495 else
496 {
497 if (j+1 >= n)
498 (*current_liboctave_error_handler) ("EIG: internal error");
499
500 m_lambda.elem (j) = FloatComplex (ar.elem (j) / beta.elem (j),
501 ai.elem (j) / beta.elem (j));
502 m_lambda.elem (j+1) = FloatComplex (ar.elem (j+1) / beta.elem (j+1),
503 ai.elem (j+1) / beta.elem (j+1));
504
505 for (F77_INT i = 0; i < nvr; i++)
506 {
507 float real_part = vr.elem (i, j);
508 float imag_part = vr.elem (i, j+1);
509 m_v.elem (i, j) = FloatComplex (real_part, imag_part);
510 m_v.elem (i, j+1) = FloatComplex (real_part, -imag_part);
511 }
512 for (F77_INT i = 0; i < nvl; i++)
513 {
514 float real_part = vl.elem (i, j);
515 float imag_part = vl.elem (i, j+1);
516 m_w.elem (i, j) = FloatComplex (real_part, imag_part);
517 m_w.elem (i, j+1) = FloatComplex (real_part, -imag_part);
518 }
519 j++;
520 }
521 }
522
523 return info;
524}
525
526octave_idx_type
527FloatEIG::symmetric_init (const FloatMatrix& a, const FloatMatrix& b,
528 bool calc_rev, bool calc_lev)
529{
530 F77_INT n = octave::to_f77_int (a.rows ());
531 F77_INT nb = octave::to_f77_int (b.rows ());
532
533 F77_INT a_nc = octave::to_f77_int (a.cols ());
534 F77_INT b_nc = octave::to_f77_int (b.cols ());
535
536 if (n != a_nc || nb != b_nc)
537 (*current_liboctave_error_handler) ("EIG requires square matrix");
538
539 if (n != nb)
540 (*current_liboctave_error_handler) ("EIG requires same size matrices");
541
542 F77_INT info = 0;
543
544 FloatMatrix atmp = a;
545 float *atmp_data = atmp.rwdata ();
546
547 FloatMatrix btmp = b;
548 float *btmp_data = btmp.rwdata ();
549
550 FloatColumnVector wr (n);
551 float *pwr = wr.rwdata ();
552
553 F77_INT lwork = -1;
554 float dummy_work;
555
556 F77_XFCN (ssygv, SSYGV, (1, F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1),
557 F77_CONST_CHAR_ARG2 ("U", 1),
558 n, atmp_data, n,
559 btmp_data, n,
560 pwr, &dummy_work, lwork, info
561 F77_CHAR_ARG_LEN (1)
562 F77_CHAR_ARG_LEN (1)));
563
564 if (info != 0)
565 (*current_liboctave_error_handler) ("ssygv workspace query failed");
566
567 lwork = static_cast<F77_INT> (dummy_work);
568 Array<float> work (dim_vector (lwork, 1));
569 float *pwork = work.rwdata ();
570
571 F77_XFCN (ssygv, SSYGV, (1, F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1),
572 F77_CONST_CHAR_ARG2 ("U", 1),
573 n, atmp_data, n,
574 btmp_data, n,
575 pwr, pwork, lwork, info
576 F77_CHAR_ARG_LEN (1)
577 F77_CHAR_ARG_LEN (1)));
578
579 if (info < 0)
580 (*current_liboctave_error_handler) ("unrecoverable error in ssygv");
581
582 if (info > 0)
583 (*current_liboctave_error_handler) ("ssygv failed to converge");
584
585 m_lambda = FloatComplexColumnVector (wr);
586 m_v = (calc_rev ? FloatComplexMatrix (atmp) : FloatComplexMatrix ());
587 m_w = (calc_lev ? FloatComplexMatrix (atmp) : FloatComplexMatrix ());
588
589 return info;
590}
591
592octave_idx_type
593FloatEIG::init (const FloatComplexMatrix& a, const FloatComplexMatrix& b,
594 bool calc_rev, bool calc_lev, bool force_qz)
595{
596 if (a.any_element_is_inf_or_nan () || b.any_element_is_inf_or_nan ())
597 (*current_liboctave_error_handler)
598 ("EIG: matrix contains Inf or NaN values");
599
600 F77_INT n = octave::to_f77_int (a.rows ());
601 F77_INT nb = octave::to_f77_int (b.rows ());
602
603 F77_INT a_nc = octave::to_f77_int (a.cols ());
604 F77_INT b_nc = octave::to_f77_int (b.cols ());
605
606 if (n != a_nc || nb != b_nc)
607 (*current_liboctave_error_handler) ("EIG requires square matrix");
608
609 if (n != nb)
610 (*current_liboctave_error_handler) ("EIG requires same size matrices");
611
612 F77_INT info = 0;
613
614 FloatComplexMatrix tmp = b;
615 FloatComplex *tmp_data = tmp.rwdata ();
616
617 if (! force_qz)
618 {
619 F77_XFCN (cpotrf, CPOTRF, (F77_CONST_CHAR_ARG2 ("L", 1),
620 n, F77_CMPLX_ARG (tmp_data), n,
621 info
622 F77_CHAR_ARG_LEN (1)));
623
624 if (a.ishermitian () && b.ishermitian () && info == 0)
625 return hermitian_init (a, b, calc_rev, calc_lev);
626 }
627
628 FloatComplexMatrix atmp = a;
629 FloatComplex *atmp_data = atmp.rwdata ();
630
631 FloatComplexMatrix btmp = b;
632 FloatComplex *btmp_data = btmp.rwdata ();
633
634 FloatComplexColumnVector alpha (n);
635 FloatComplex *palpha = alpha.rwdata ();
636
637 FloatComplexColumnVector beta (n);
638 FloatComplex *pbeta = beta.rwdata ();
639
640 F77_INT nvr = (calc_rev ? n : 0);
641 FloatComplexMatrix vrtmp (nvr, nvr);
642 FloatComplex *pvr = vrtmp.rwdata ();
643
644 F77_INT nvl = (calc_lev ? n : 0);
645 FloatComplexMatrix vltmp (nvl, nvl);
646 FloatComplex *pvl = vltmp.rwdata ();
647
648 F77_INT lwork = -1;
649 FloatComplex dummy_work;
650
651 F77_INT lrwork = 8*n;
652 Array<float> rwork (dim_vector (lrwork, 1));
653 float *prwork = rwork.rwdata ();
654
655 F77_XFCN (cggev, CGGEV, (F77_CONST_CHAR_ARG2 (calc_lev ? "V" : "N", 1),
656 F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1),
657 n, F77_CMPLX_ARG (atmp_data), n,
658 F77_CMPLX_ARG (btmp_data), n,
659 F77_CMPLX_ARG (palpha), F77_CMPLX_ARG (pbeta),
660 F77_CMPLX_ARG (pvl), n, F77_CMPLX_ARG (pvr), n,
661 F77_CMPLX_ARG (&dummy_work), lwork, prwork, info
662 F77_CHAR_ARG_LEN (1)
663 F77_CHAR_ARG_LEN (1)));
664
665 if (info != 0)
666 (*current_liboctave_error_handler) ("cggev workspace query failed");
667
668 lwork = static_cast<F77_INT> (dummy_work.real ());
669 Array<FloatComplex> work (dim_vector (lwork, 1));
670 FloatComplex *pwork = work.rwdata ();
671
672 F77_XFCN (cggev, CGGEV, (F77_CONST_CHAR_ARG2 (calc_lev ? "V" : "N", 1),
673 F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1),
674 n, F77_CMPLX_ARG (atmp_data), n,
675 F77_CMPLX_ARG (btmp_data), n,
676 F77_CMPLX_ARG (palpha), F77_CMPLX_ARG (pbeta),
677 F77_CMPLX_ARG (pvl), n, F77_CMPLX_ARG (pvr), n,
678 F77_CMPLX_ARG (pwork), lwork, prwork, info
679 F77_CHAR_ARG_LEN (1)
680 F77_CHAR_ARG_LEN (1)));
681
682 if (info < 0)
683 (*current_liboctave_error_handler) ("unrecoverable error in cggev");
684
685 if (info > 0)
686 (*current_liboctave_error_handler) ("cggev failed to converge");
687
688 m_lambda.resize (n);
689
690 for (F77_INT j = 0; j < n; j++)
691 m_lambda.elem (j) = alpha.elem (j) / beta.elem (j);
692
693 m_v = vrtmp;
694 m_w = vltmp;
695
696 return info;
697}
698
699octave_idx_type
700FloatEIG::hermitian_init (const FloatComplexMatrix& a,
701 const FloatComplexMatrix& b,
702 bool calc_rev, bool calc_lev)
703{
704 F77_INT n = octave::to_f77_int (a.rows ());
705 F77_INT nb = octave::to_f77_int (b.rows ());
706
707 F77_INT a_nc = octave::to_f77_int (a.cols ());
708 F77_INT b_nc = octave::to_f77_int (b.cols ());
709
710 if (n != a_nc || nb != b_nc)
711 (*current_liboctave_error_handler) ("EIG requires square matrix");
712
713 if (n != nb)
714 (*current_liboctave_error_handler) ("EIG requires same size matrices");
715
716 F77_INT info = 0;
717
718 FloatComplexMatrix atmp = a;
719 FloatComplex *atmp_data = atmp.rwdata ();
720
721 FloatComplexMatrix btmp = b;
722 FloatComplex *btmp_data = btmp.rwdata ();
723
724 FloatColumnVector wr (n);
725 float *pwr = wr.rwdata ();
726
727 F77_INT lwork = -1;
728 FloatComplex dummy_work;
729
730 F77_INT lrwork = 3*n;
731 Array<float> rwork (dim_vector (lrwork, 1));
732 float *prwork = rwork.rwdata ();
733
734 F77_XFCN (chegv, CHEGV, (1, F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1),
735 F77_CONST_CHAR_ARG2 ("U", 1),
736 n, F77_CMPLX_ARG (atmp_data), n,
737 F77_CMPLX_ARG (btmp_data), n,
738 pwr, F77_CMPLX_ARG (&dummy_work), lwork,
739 prwork, info
740 F77_CHAR_ARG_LEN (1)
741 F77_CHAR_ARG_LEN (1)));
742
743 if (info != 0)
744 (*current_liboctave_error_handler) ("zhegv workspace query failed");
745
746 lwork = static_cast<F77_INT> (dummy_work.real ());
747 Array<FloatComplex> work (dim_vector (lwork, 1));
748 FloatComplex *pwork = work.rwdata ();
749
750 F77_XFCN (chegv, CHEGV, (1, F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1),
751 F77_CONST_CHAR_ARG2 ("U", 1),
752 n, F77_CMPLX_ARG (atmp_data), n,
753 F77_CMPLX_ARG (btmp_data), n,
754 pwr, F77_CMPLX_ARG (pwork), lwork, prwork, info
755 F77_CHAR_ARG_LEN (1)
756 F77_CHAR_ARG_LEN (1)));
757
758 if (info < 0)
759 (*current_liboctave_error_handler) ("unrecoverable error in zhegv");
760
761 if (info > 0)
762 (*current_liboctave_error_handler) ("zhegv failed to converge");
763
764 m_lambda = FloatComplexColumnVector (wr);
765 m_v = (calc_rev ? FloatComplexMatrix (atmp) : FloatComplexMatrix ());
766 m_w = (calc_lev ? FloatComplexMatrix (atmp) : FloatComplexMatrix ());
767
768 return info;
769}
Array
N Dimensional Array with copy-on-write semantics.
Definition Array.h:130
Array::elem
T & elem(octave_idx_type n)
Size of the specified dimension.
Definition Array.h:563
Array::rows
octave_idx_type rows() const
Definition Array.h:463
Array::cols
octave_idx_type cols() const
Definition Array.h:473
Array::rwdata
T * rwdata()
Size of the specified dimension.
Definition Array-base.cc:1787
FloatComplexColumnVector::resize
void resize(octave_idx_type n, const FloatComplex &rfv=FloatComplex(0))
Definition fCColVector.h:152
FloatComplexMatrix::resize
void resize(octave_idx_type nr, octave_idx_type nc, const FloatComplex &rfv=FloatComplex(0))
Definition fCMatrix.h:199
FloatComplexMatrix::ishermitian
bool ishermitian() const
Definition fCMatrix.cc:175
FloatComplexNDArray::any_element_is_inf_or_nan
bool any_element_is_inf_or_nan() const
Definition fCNDArray.cc:271
FloatEIG::FloatComplexMatrix
friend class FloatComplexMatrix
Definition fEIG.h:41
FloatMatrix::issymmetric
bool issymmetric() const
Definition fMatrix.cc:139
FloatNDArray::any_element_is_inf_or_nan
bool any_element_is_inf_or_nan() const
Definition fNDArray.cc:281
dim_vector
Vector representing the dimensions (size) of an Array.
Definition dim-vector.h:90
octave_idx_type
F77_CMPLX_ARG
#define F77_CMPLX_ARG(x)
Definition f77-fcn.h:310
F77_XFCN
#define F77_XFCN(f, F, args)
Definition f77-fcn.h:45
F77_INT
octave_f77_int_type F77_INT
Definition f77-fcn.h:306
scale
void scale(Matrix &m, double x, double y, double z)
Definition graphics.cc:5731
FloatComplex
std::complex< float > FloatComplex
Definition oct-cmplx.h:34