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