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