Mercurial > octave
annotate liboctave/numeric/EIG.cc @ 31249:de6fc38c78c6
Make Jacobian types offered by dlsode.f accessible by lsode (bug #31626).
* liboctave/numeric/LSODE-opts.in: Add options "jacobian type", "lower jacobian
subdiagonals", and "upper jacobian subdiagonals".
* liboctave/numeric/LSODE.cc (file scope, lsode_j,
LSODE::do_integrate (double)): Handle new configurable Jacobian types.
* build-aux/mk-opts.pl: Don't implicitly convert to integer in condition.
author | Olaf Till <olaf.till@uni-jena.de> |
---|---|
date | Fri, 12 Nov 2010 08:53:05 +0100 |
parents | 796f54d4ddbf |
children | 597f3ee61a48 |
rev | line source |
---|---|
27923
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
1 //////////////////////////////////////////////////////////////////////// |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
2 // |
30564
796f54d4ddbf
update Octave Project Developers copyright for the new year
John W. Eaton <jwe@octave.org>
parents:
30394
diff
changeset
|
3 // Copyright (C) 1994-2022 The Octave Project Developers |
27923
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
4 // |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
5 // See the file COPYRIGHT.md in the top-level directory of this |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
6 // distribution or <https://octave.org/copyright/>. |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
7 // |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
8 // This file is part of Octave. |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
9 // |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
10 // Octave is free software: you can redistribute it and/or modify it |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
11 // under the terms of the GNU General Public License as published by |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
12 // the Free Software Foundation, either version 3 of the License, or |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
13 // (at your option) any later version. |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
14 // |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
15 // Octave is distributed in the hope that it will be useful, but |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
16 // WITHOUT ANY WARRANTY; without even the implied warranty of |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
18 // GNU General Public License for more details. |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
19 // |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
20 // You should have received a copy of the GNU General Public License |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
21 // along with Octave; see the file COPYING. If not, see |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
22 // <https://www.gnu.org/licenses/>. |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
23 // |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
24 //////////////////////////////////////////////////////////////////////// |
462 | 25 |
21724
aba2e6293dd8
use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents:
21301
diff
changeset
|
26 #if defined (HAVE_CONFIG_H) |
21301
40de9f8f23a6
Use '#include "config.h"' rather than <config.h>.
Rik <rik@octave.org>
parents:
21202
diff
changeset
|
27 # include "config.h" |
462 | 28 #endif |
29 | |
23475
d691ed308237
maint: Clean up #includes in liboctave/numeric directory.
Rik <rik@octave.org>
parents:
23450
diff
changeset
|
30 #include "Array.h" |
462 | 31 #include "EIG.h" |
2815 | 32 #include "dColVector.h" |
23475
d691ed308237
maint: Clean up #includes in liboctave/numeric directory.
Rik <rik@octave.org>
parents:
23450
diff
changeset
|
33 #include "dMatrix.h" |
462 | 34 #include "lo-error.h" |
22322
93b3cdd36854
move most f77 function decls to separate header files
John W. Eaton <jwe@octave.org>
parents:
22305
diff
changeset
|
35 #include "lo-lapack-proto.h" |
462 | 36 |
5275 | 37 octave_idx_type |
22305
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
38 EIG::init (const Matrix& a, bool calc_rev, bool calc_lev, bool balance) |
462 | 39 { |
5822 | 40 if (a.any_element_is_inf_or_nan ()) |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
41 (*current_liboctave_error_handler) |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
42 ("EIG: matrix contains Inf or NaN values"); |
5822 | 43 |
23596
b63c3a09aee7
maint: Deprecate is_symmetric and replace with issymmetric.
Rik <rik@octave.org>
parents:
23594
diff
changeset
|
44 if (a.issymmetric ()) |
22305
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
45 return symmetric_init (a, calc_rev, calc_lev); |
2815 | 46 |
22988
cd33c785e80e
put to_f77_int inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
22953
diff
changeset
|
47 F77_INT n = octave::to_f77_int (a.rows ()); |
cd33c785e80e
put to_f77_int inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
22953
diff
changeset
|
48 F77_INT a_nc = octave::to_f77_int (a.cols ()); |
1934 | 49 |
22953
fd649fd3db75
use F77_INT instead of octave_idx_type for liboctave EIG classes
John W. Eaton <jwe@octave.org>
parents:
22755
diff
changeset
|
50 if (n != a_nc) |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
51 (*current_liboctave_error_handler) ("EIG requires square matrix"); |
462 | 52 |
22953
fd649fd3db75
use F77_INT instead of octave_idx_type for liboctave EIG classes
John W. Eaton <jwe@octave.org>
parents:
22755
diff
changeset
|
53 F77_INT info = 0; |
462 | 54 |
1934 | 55 Matrix atmp = a; |
56 double *tmp_data = atmp.fortran_vec (); | |
57 | |
11570
57632dea2446
attempt better backward compatibility for Array constructors
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
58 Array<double> wr (dim_vector (n, 1)); |
1934 | 59 double *pwr = wr.fortran_vec (); |
60 | |
11570
57632dea2446
attempt better backward compatibility for Array constructors
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
61 Array<double> wi (dim_vector (n, 1)); |
1934 | 62 double *pwi = wi.fortran_vec (); |
63 | |
23450
855122b993da
maint: Wrap tertiary operator in parentheses "(COND ? x : y)".
Rik <rik@octave.org>
parents:
23220
diff
changeset
|
64 F77_INT tnvr = (calc_rev ? n : 0); |
9020
728e7943752d
EIG.cc: avoid volatile decl for tmp variable
John W. Eaton <jwe@octave.org>
parents:
8920
diff
changeset
|
65 Matrix vr (tnvr, tnvr); |
462 | 66 double *pvr = vr.fortran_vec (); |
1934 | 67 |
23450
855122b993da
maint: Wrap tertiary operator in parentheses "(COND ? x : y)".
Rik <rik@octave.org>
parents:
23220
diff
changeset
|
68 F77_INT tnvl = (calc_lev ? n : 0); |
22305
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
69 Matrix vl (tnvl, tnvl); |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
70 double *pvl = vl.fortran_vec (); |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
71 |
22953
fd649fd3db75
use F77_INT instead of octave_idx_type for liboctave EIG classes
John W. Eaton <jwe@octave.org>
parents:
22755
diff
changeset
|
72 F77_INT lwork = -1; |
4800 | 73 double dummy_work; |
462 | 74 |
22953
fd649fd3db75
use F77_INT instead of octave_idx_type for liboctave EIG classes
John W. Eaton <jwe@octave.org>
parents:
22755
diff
changeset
|
75 F77_INT ilo; |
fd649fd3db75
use F77_INT instead of octave_idx_type for liboctave EIG classes
John W. Eaton <jwe@octave.org>
parents:
22755
diff
changeset
|
76 F77_INT ihi; |
22305
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
77 |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
78 Array<double> scale (dim_vector (n, 1)); |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
79 double *pscale = scale.fortran_vec (); |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
80 |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
81 double abnrm; |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
82 |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
83 Array<double> rconde (dim_vector (n, 1)); |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
84 double *prconde = rconde.fortran_vec (); |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
85 |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
86 Array<double> rcondv (dim_vector (n, 1)); |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
87 double *prcondv = rcondv.fortran_vec (); |
462 | 88 |
22953
fd649fd3db75
use F77_INT instead of octave_idx_type for liboctave EIG classes
John W. Eaton <jwe@octave.org>
parents:
22755
diff
changeset
|
89 F77_INT dummy_iwork; |
22305
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
90 |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
91 F77_XFCN (dgeevx, DGEEVX, (F77_CONST_CHAR_ARG2 (balance ? "B" : "N", 1), |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
92 F77_CONST_CHAR_ARG2 (calc_lev ? "V" : "N", 1), |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
93 F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1), |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
94 F77_CONST_CHAR_ARG2 ("N", 1), |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
95 n, tmp_data, n, pwr, pwi, pvl, |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
96 n, pvr, n, ilo, ihi, pscale, |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
97 abnrm, prconde, prcondv, &dummy_work, |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
98 lwork, &dummy_iwork, info |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
99 F77_CHAR_ARG_LEN (1) |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
100 F77_CHAR_ARG_LEN (1) |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
101 F77_CHAR_ARG_LEN (1) |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
102 F77_CHAR_ARG_LEN (1))); |
462 | 103 |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
104 if (info != 0) |
22305
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
105 (*current_liboctave_error_handler) ("dgeevx workspace query failed"); |
4800 | 106 |
22953
fd649fd3db75
use F77_INT instead of octave_idx_type for liboctave EIG classes
John W. Eaton <jwe@octave.org>
parents:
22755
diff
changeset
|
107 lwork = static_cast<F77_INT> (dummy_work); |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
108 Array<double> work (dim_vector (lwork, 1)); |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
109 double *pwork = work.fortran_vec (); |
4800 | 110 |
22305
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
111 F77_XFCN (dgeevx, DGEEVX, (F77_CONST_CHAR_ARG2 (balance ? "B" : "N", 1), |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
112 F77_CONST_CHAR_ARG2 (calc_lev ? "V" : "N", 1), |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
113 F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1), |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
114 F77_CONST_CHAR_ARG2 ("N", 1), |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
115 n, tmp_data, n, pwr, pwi, pvl, |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
116 n, pvr, n, ilo, ihi, pscale, |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
117 abnrm, prconde, prcondv, pwork, |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
118 lwork, &dummy_iwork, info |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
119 F77_CHAR_ARG_LEN (1) |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
120 F77_CHAR_ARG_LEN (1) |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
121 F77_CHAR_ARG_LEN (1) |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
122 F77_CHAR_ARG_LEN (1))); |
4800 | 123 |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
124 if (info < 0) |
22305
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
125 (*current_liboctave_error_handler) ("unrecoverable error in dgeevx"); |
1934 | 126 |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
127 if (info > 0) |
22305
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
128 (*current_liboctave_error_handler) ("dgeevx failed to converge"); |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
129 |
30044
1f34286a2637
maint: use "m_" prefix for member variables of classes associated with EIG.
Rik <rik@octave.org>
parents:
29358
diff
changeset
|
130 m_lambda.resize (n); |
23450
855122b993da
maint: Wrap tertiary operator in parentheses "(COND ? x : y)".
Rik <rik@octave.org>
parents:
23220
diff
changeset
|
131 F77_INT nvr = (calc_rev ? n : 0); |
30044
1f34286a2637
maint: use "m_" prefix for member variables of classes associated with EIG.
Rik <rik@octave.org>
parents:
29358
diff
changeset
|
132 m_v.resize (nvr, nvr); |
23450
855122b993da
maint: Wrap tertiary operator in parentheses "(COND ? x : y)".
Rik <rik@octave.org>
parents:
23220
diff
changeset
|
133 F77_INT nvl = (calc_lev ? n : 0); |
30044
1f34286a2637
maint: use "m_" prefix for member variables of classes associated with EIG.
Rik <rik@octave.org>
parents:
29358
diff
changeset
|
134 m_w.resize (nvl, nvl); |
4800 | 135 |
22953
fd649fd3db75
use F77_INT instead of octave_idx_type for liboctave EIG classes
John W. Eaton <jwe@octave.org>
parents:
22755
diff
changeset
|
136 for (F77_INT j = 0; j < n; j++) |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
137 { |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
138 if (wi.elem (j) == 0.0) |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
139 { |
30044
1f34286a2637
maint: use "m_" prefix for member variables of classes associated with EIG.
Rik <rik@octave.org>
parents:
29358
diff
changeset
|
140 m_lambda.elem (j) = Complex (wr.elem (j)); |
22953
fd649fd3db75
use F77_INT instead of octave_idx_type for liboctave EIG classes
John W. Eaton <jwe@octave.org>
parents:
22755
diff
changeset
|
141 for (F77_INT i = 0; i < nvr; i++) |
30044
1f34286a2637
maint: use "m_" prefix for member variables of classes associated with EIG.
Rik <rik@octave.org>
parents:
29358
diff
changeset
|
142 m_v.elem (i, j) = vr.elem (i, j); |
22305
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
143 |
22953
fd649fd3db75
use F77_INT instead of octave_idx_type for liboctave EIG classes
John W. Eaton <jwe@octave.org>
parents:
22755
diff
changeset
|
144 for (F77_INT i = 0; i < nvl; i++) |
30044
1f34286a2637
maint: use "m_" prefix for member variables of classes associated with EIG.
Rik <rik@octave.org>
parents:
29358
diff
changeset
|
145 m_w.elem (i, j) = vl.elem (i, j); |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
146 } |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
147 else |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
148 { |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
149 if (j+1 >= n) |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
150 (*current_liboctave_error_handler) ("EIG: internal error"); |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
151 |
30044
1f34286a2637
maint: use "m_" prefix for member variables of classes associated with EIG.
Rik <rik@octave.org>
parents:
29358
diff
changeset
|
152 m_lambda.elem (j) = Complex (wr.elem (j), wi.elem (j)); |
1f34286a2637
maint: use "m_" prefix for member variables of classes associated with EIG.
Rik <rik@octave.org>
parents:
29358
diff
changeset
|
153 m_lambda.elem (j+1) = Complex (wr.elem (j+1), wi.elem (j+1)); |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
154 |
22953
fd649fd3db75
use F77_INT instead of octave_idx_type for liboctave EIG classes
John W. Eaton <jwe@octave.org>
parents:
22755
diff
changeset
|
155 for (F77_INT i = 0; i < nvr; i++) |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
156 { |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
157 double real_part = vr.elem (i, j); |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
158 double imag_part = vr.elem (i, j+1); |
30044
1f34286a2637
maint: use "m_" prefix for member variables of classes associated with EIG.
Rik <rik@octave.org>
parents:
29358
diff
changeset
|
159 m_v.elem (i, j) = Complex (real_part, imag_part); |
1f34286a2637
maint: use "m_" prefix for member variables of classes associated with EIG.
Rik <rik@octave.org>
parents:
29358
diff
changeset
|
160 m_v.elem (i, j+1) = Complex (real_part, -imag_part); |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
161 } |
22305
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
162 |
22953
fd649fd3db75
use F77_INT instead of octave_idx_type for liboctave EIG classes
John W. Eaton <jwe@octave.org>
parents:
22755
diff
changeset
|
163 for (F77_INT i = 0; i < nvl; i++) |
22305
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
164 { |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
165 double real_part = vl.elem (i, j); |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
166 double imag_part = vl.elem (i, j+1); |
30044
1f34286a2637
maint: use "m_" prefix for member variables of classes associated with EIG.
Rik <rik@octave.org>
parents:
29358
diff
changeset
|
167 m_w.elem (i, j) = Complex (real_part, imag_part); |
1f34286a2637
maint: use "m_" prefix for member variables of classes associated with EIG.
Rik <rik@octave.org>
parents:
29358
diff
changeset
|
168 m_w.elem (i, j+1) = Complex (real_part, -imag_part); |
22305
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
169 } |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
170 j++; |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
171 } |
2815 | 172 } |
173 | |
174 return info; | |
175 } | |
176 | |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11570
diff
changeset
|
177 octave_idx_type |
22305
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
178 EIG::symmetric_init (const Matrix& a, bool calc_rev, bool calc_lev) |
2815 | 179 { |
22988
cd33c785e80e
put to_f77_int inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
22953
diff
changeset
|
180 F77_INT n = octave::to_f77_int (a.rows ()); |
cd33c785e80e
put to_f77_int inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
22953
diff
changeset
|
181 F77_INT a_nc = octave::to_f77_int (a.cols ()); |
1934 | 182 |
22953
fd649fd3db75
use F77_INT instead of octave_idx_type for liboctave EIG classes
John W. Eaton <jwe@octave.org>
parents:
22755
diff
changeset
|
183 if (n != a_nc) |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
184 (*current_liboctave_error_handler) ("EIG requires square matrix"); |
2815 | 185 |
22953
fd649fd3db75
use F77_INT instead of octave_idx_type for liboctave EIG classes
John W. Eaton <jwe@octave.org>
parents:
22755
diff
changeset
|
186 F77_INT info = 0; |
2815 | 187 |
188 Matrix atmp = a; | |
189 double *tmp_data = atmp.fortran_vec (); | |
190 | |
3585 | 191 ColumnVector wr (n); |
2815 | 192 double *pwr = wr.fortran_vec (); |
193 | |
22953
fd649fd3db75
use F77_INT instead of octave_idx_type for liboctave EIG classes
John W. Eaton <jwe@octave.org>
parents:
22755
diff
changeset
|
194 F77_INT lwork = -1; |
4800 | 195 double dummy_work; |
2815 | 196 |
22305
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
197 F77_XFCN (dsyev, DSYEV, (F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1), |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
198 F77_CONST_CHAR_ARG2 ("U", 1), |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
199 n, tmp_data, n, pwr, &dummy_work, lwork, info |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
200 F77_CHAR_ARG_LEN (1) |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
201 F77_CHAR_ARG_LEN (1))); |
2815 | 202 |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
203 if (info != 0) |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
204 (*current_liboctave_error_handler) ("dsyev workspace query failed"); |
4800 | 205 |
22953
fd649fd3db75
use F77_INT instead of octave_idx_type for liboctave EIG classes
John W. Eaton <jwe@octave.org>
parents:
22755
diff
changeset
|
206 lwork = static_cast<F77_INT> (dummy_work); |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
207 Array<double> work (dim_vector (lwork, 1)); |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
208 double *pwork = work.fortran_vec (); |
4800 | 209 |
22305
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
210 F77_XFCN (dsyev, DSYEV, (F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1), |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
211 F77_CONST_CHAR_ARG2 ("U", 1), |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
212 n, tmp_data, n, pwr, pwork, lwork, info |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
213 F77_CHAR_ARG_LEN (1) |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
214 F77_CHAR_ARG_LEN (1))); |
4800 | 215 |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
216 if (info < 0) |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
217 (*current_liboctave_error_handler) ("unrecoverable error in dsyev"); |
4800 | 218 |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
219 if (info > 0) |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
220 (*current_liboctave_error_handler) ("dsyev failed to converge"); |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
221 |
30044
1f34286a2637
maint: use "m_" prefix for member variables of classes associated with EIG.
Rik <rik@octave.org>
parents:
29358
diff
changeset
|
222 m_lambda = ComplexColumnVector (wr); |
1f34286a2637
maint: use "m_" prefix for member variables of classes associated with EIG.
Rik <rik@octave.org>
parents:
29358
diff
changeset
|
223 m_v = (calc_rev ? ComplexMatrix (atmp) : ComplexMatrix ()); |
1f34286a2637
maint: use "m_" prefix for member variables of classes associated with EIG.
Rik <rik@octave.org>
parents:
29358
diff
changeset
|
224 m_w = (calc_lev ? ComplexMatrix (atmp) : ComplexMatrix ()); |
462 | 225 |
226 return info; | |
227 } | |
228 | |
5275 | 229 octave_idx_type |
22305
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
230 EIG::init (const ComplexMatrix& a, bool calc_rev, bool calc_lev, bool balance) |
462 | 231 { |
5822 | 232 if (a.any_element_is_inf_or_nan ()) |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
233 (*current_liboctave_error_handler) |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
234 ("EIG: matrix contains Inf or NaN values"); |
5822 | 235 |
23594
af5b813503cb
maint: Deprecate is_hermitian and replace with ishermitian.
Rik <rik@octave.org>
parents:
23475
diff
changeset
|
236 if (a.ishermitian ()) |
22305
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
237 return hermitian_init (a, calc_rev, calc_lev); |
2815 | 238 |
22988
cd33c785e80e
put to_f77_int inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
22953
diff
changeset
|
239 F77_INT n = octave::to_f77_int (a.rows ()); |
cd33c785e80e
put to_f77_int inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
22953
diff
changeset
|
240 F77_INT a_nc = octave::to_f77_int (a.cols ()); |
1934 | 241 |
22953
fd649fd3db75
use F77_INT instead of octave_idx_type for liboctave EIG classes
John W. Eaton <jwe@octave.org>
parents:
22755
diff
changeset
|
242 if (n != a_nc) |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
243 (*current_liboctave_error_handler) ("EIG requires square matrix"); |
462 | 244 |
22953
fd649fd3db75
use F77_INT instead of octave_idx_type for liboctave EIG classes
John W. Eaton <jwe@octave.org>
parents:
22755
diff
changeset
|
245 F77_INT info = 0; |
462 | 246 |
1934 | 247 ComplexMatrix atmp = a; |
248 Complex *tmp_data = atmp.fortran_vec (); | |
462 | 249 |
22305
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
250 ComplexColumnVector wr (n); |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
251 Complex *pw = wr.fortran_vec (); |
2815 | 252 |
23450
855122b993da
maint: Wrap tertiary operator in parentheses "(COND ? x : y)".
Rik <rik@octave.org>
parents:
23220
diff
changeset
|
253 F77_INT nvr = (calc_rev ? n : 0); |
22305
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
254 ComplexMatrix vrtmp (nvr, nvr); |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
255 Complex *pvr = vrtmp.fortran_vec (); |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
256 |
23450
855122b993da
maint: Wrap tertiary operator in parentheses "(COND ? x : y)".
Rik <rik@octave.org>
parents:
23220
diff
changeset
|
257 F77_INT nvl = (calc_lev ? n : 0); |
22305
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
258 ComplexMatrix vltmp (nvl, nvl); |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
259 Complex *pvl = vltmp.fortran_vec (); |
2815 | 260 |
22953
fd649fd3db75
use F77_INT instead of octave_idx_type for liboctave EIG classes
John W. Eaton <jwe@octave.org>
parents:
22755
diff
changeset
|
261 F77_INT lwork = -1; |
4800 | 262 Complex dummy_work; |
1934 | 263 |
22953
fd649fd3db75
use F77_INT instead of octave_idx_type for liboctave EIG classes
John W. Eaton <jwe@octave.org>
parents:
22755
diff
changeset
|
264 F77_INT lrwork = 2*n; |
11570
57632dea2446
attempt better backward compatibility for Array constructors
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
265 Array<double> rwork (dim_vector (lrwork, 1)); |
1934 | 266 double *prwork = rwork.fortran_vec (); |
462 | 267 |
22953
fd649fd3db75
use F77_INT instead of octave_idx_type for liboctave EIG classes
John W. Eaton <jwe@octave.org>
parents:
22755
diff
changeset
|
268 F77_INT ilo; |
fd649fd3db75
use F77_INT instead of octave_idx_type for liboctave EIG classes
John W. Eaton <jwe@octave.org>
parents:
22755
diff
changeset
|
269 F77_INT ihi; |
22305
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
270 |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
271 Array<double> scale (dim_vector (n, 1)); |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
272 double *pscale = scale.fortran_vec (); |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
273 |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
274 double abnrm; |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
275 |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
276 Array<double> rconde (dim_vector (n, 1)); |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
277 double *prconde = rconde.fortran_vec (); |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
278 |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
279 Array<double> rcondv (dim_vector (n, 1)); |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
280 double *prcondv = rcondv.fortran_vec (); |
462 | 281 |
22305
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
282 F77_XFCN (zgeevx, ZGEEVX, (F77_CONST_CHAR_ARG2 (balance ? "B" : "N", 1), |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
283 F77_CONST_CHAR_ARG2 (calc_lev ? "V" : "N", 1), |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
284 F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1), |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
285 F77_CONST_CHAR_ARG2 ("N", 1), |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
286 n, F77_DBLE_CMPLX_ARG (tmp_data), n, |
23475
d691ed308237
maint: Clean up #includes in liboctave/numeric directory.
Rik <rik@octave.org>
parents:
23450
diff
changeset
|
287 F77_DBLE_CMPLX_ARG (pw), F77_DBLE_CMPLX_ARG (pvl), |
d691ed308237
maint: Clean up #includes in liboctave/numeric directory.
Rik <rik@octave.org>
parents:
23450
diff
changeset
|
288 n, F77_DBLE_CMPLX_ARG (pvr), n, ilo, ihi, |
d691ed308237
maint: Clean up #includes in liboctave/numeric directory.
Rik <rik@octave.org>
parents:
23450
diff
changeset
|
289 pscale, abnrm, prconde, prcondv, |
d691ed308237
maint: Clean up #includes in liboctave/numeric directory.
Rik <rik@octave.org>
parents:
23450
diff
changeset
|
290 F77_DBLE_CMPLX_ARG (&dummy_work), lwork, prwork, |
d691ed308237
maint: Clean up #includes in liboctave/numeric directory.
Rik <rik@octave.org>
parents:
23450
diff
changeset
|
291 info |
22305
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
292 F77_CHAR_ARG_LEN (1) |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
293 F77_CHAR_ARG_LEN (1) |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
294 F77_CHAR_ARG_LEN (1) |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
295 F77_CHAR_ARG_LEN (1))); |
2815 | 296 |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
297 if (info != 0) |
22305
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
298 (*current_liboctave_error_handler) ("zgeevx workspace query failed"); |
4800 | 299 |
22953
fd649fd3db75
use F77_INT instead of octave_idx_type for liboctave EIG classes
John W. Eaton <jwe@octave.org>
parents:
22755
diff
changeset
|
300 lwork = static_cast<F77_INT> (dummy_work.real ()); |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
301 Array<Complex> work (dim_vector (lwork, 1)); |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
302 Complex *pwork = work.fortran_vec (); |
4800 | 303 |
22305
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
304 F77_XFCN (zgeevx, ZGEEVX, (F77_CONST_CHAR_ARG2 (balance ? "B" : "N", 1), |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
305 F77_CONST_CHAR_ARG2 (calc_lev ? "V" : "N", 1), |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
306 F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1), |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
307 F77_CONST_CHAR_ARG2 ("N", 1), |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
308 n, F77_DBLE_CMPLX_ARG (tmp_data), n, |
23475
d691ed308237
maint: Clean up #includes in liboctave/numeric directory.
Rik <rik@octave.org>
parents:
23450
diff
changeset
|
309 F77_DBLE_CMPLX_ARG (pw), F77_DBLE_CMPLX_ARG (pvl), |
d691ed308237
maint: Clean up #includes in liboctave/numeric directory.
Rik <rik@octave.org>
parents:
23450
diff
changeset
|
310 n, F77_DBLE_CMPLX_ARG (pvr), n, ilo, ihi, |
d691ed308237
maint: Clean up #includes in liboctave/numeric directory.
Rik <rik@octave.org>
parents:
23450
diff
changeset
|
311 pscale, abnrm, prconde, prcondv, |
22305
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
312 F77_DBLE_CMPLX_ARG (pwork), lwork, prwork, info |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
313 F77_CHAR_ARG_LEN (1) |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
314 F77_CHAR_ARG_LEN (1) |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
315 F77_CHAR_ARG_LEN (1) |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
316 F77_CHAR_ARG_LEN (1))); |
4800 | 317 |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
318 if (info < 0) |
22305
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
319 (*current_liboctave_error_handler) ("unrecoverable error in zgeevx"); |
4800 | 320 |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
321 if (info > 0) |
22305
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
322 (*current_liboctave_error_handler) ("zgeevx failed to converge"); |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
323 |
30044
1f34286a2637
maint: use "m_" prefix for member variables of classes associated with EIG.
Rik <rik@octave.org>
parents:
29358
diff
changeset
|
324 m_lambda = wr; |
1f34286a2637
maint: use "m_" prefix for member variables of classes associated with EIG.
Rik <rik@octave.org>
parents:
29358
diff
changeset
|
325 m_v = vrtmp; |
1f34286a2637
maint: use "m_" prefix for member variables of classes associated with EIG.
Rik <rik@octave.org>
parents:
29358
diff
changeset
|
326 m_w = vltmp; |
2815 | 327 |
328 return info; | |
329 } | |
330 | |
5275 | 331 octave_idx_type |
22305
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
332 EIG::hermitian_init (const ComplexMatrix& a, bool calc_rev, bool calc_lev) |
2815 | 333 { |
22988
cd33c785e80e
put to_f77_int inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
22953
diff
changeset
|
334 F77_INT n = octave::to_f77_int (a.rows ()); |
cd33c785e80e
put to_f77_int inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
22953
diff
changeset
|
335 F77_INT a_nc = octave::to_f77_int (a.cols ()); |
2815 | 336 |
22953
fd649fd3db75
use F77_INT instead of octave_idx_type for liboctave EIG classes
John W. Eaton <jwe@octave.org>
parents:
22755
diff
changeset
|
337 if (n != a_nc) |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
338 (*current_liboctave_error_handler) ("EIG requires square matrix"); |
2815 | 339 |
22953
fd649fd3db75
use F77_INT instead of octave_idx_type for liboctave EIG classes
John W. Eaton <jwe@octave.org>
parents:
22755
diff
changeset
|
340 F77_INT info = 0; |
462 | 341 |
2815 | 342 ComplexMatrix atmp = a; |
343 Complex *tmp_data = atmp.fortran_vec (); | |
344 | |
3585 | 345 ColumnVector wr (n); |
346 double *pwr = wr.fortran_vec (); | |
2815 | 347 |
22953
fd649fd3db75
use F77_INT instead of octave_idx_type for liboctave EIG classes
John W. Eaton <jwe@octave.org>
parents:
22755
diff
changeset
|
348 F77_INT lwork = -1; |
4800 | 349 Complex dummy_work; |
2815 | 350 |
22953
fd649fd3db75
use F77_INT instead of octave_idx_type for liboctave EIG classes
John W. Eaton <jwe@octave.org>
parents:
22755
diff
changeset
|
351 F77_INT lrwork = 3*n; |
11570
57632dea2446
attempt better backward compatibility for Array constructors
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
352 Array<double> rwork (dim_vector (lrwork, 1)); |
2815 | 353 double *prwork = rwork.fortran_vec (); |
354 | |
22305
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
355 F77_XFCN (zheev, ZHEEV, (F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1), |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
356 F77_CONST_CHAR_ARG2 ("U", 1), |
22305
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
357 n, F77_DBLE_CMPLX_ARG (tmp_data), n, pwr, |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
358 F77_DBLE_CMPLX_ARG (&dummy_work), lwork, |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
359 prwork, info |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
360 F77_CHAR_ARG_LEN (1) |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
361 F77_CHAR_ARG_LEN (1))); |
2815 | 362 |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
363 if (info != 0) |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
364 (*current_liboctave_error_handler) ("zheev workspace query failed"); |
4800 | 365 |
22953
fd649fd3db75
use F77_INT instead of octave_idx_type for liboctave EIG classes
John W. Eaton <jwe@octave.org>
parents:
22755
diff
changeset
|
366 lwork = static_cast<F77_INT> (dummy_work.real ()); |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
367 Array<Complex> work (dim_vector (lwork, 1)); |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
368 Complex *pwork = work.fortran_vec (); |
4800 | 369 |
22305
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
370 F77_XFCN (zheev, ZHEEV, (F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1), |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
371 F77_CONST_CHAR_ARG2 ("U", 1), |
22305
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
372 n, F77_DBLE_CMPLX_ARG (tmp_data), n, pwr, |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
373 F77_DBLE_CMPLX_ARG (pwork), lwork, prwork, info |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
374 F77_CHAR_ARG_LEN (1) |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
375 F77_CHAR_ARG_LEN (1))); |
4800 | 376 |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
377 if (info < 0) |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
378 (*current_liboctave_error_handler) ("unrecoverable error in zheev"); |
4800 | 379 |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
380 if (info > 0) |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
381 (*current_liboctave_error_handler) ("zheev failed to converge"); |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
382 |
30044
1f34286a2637
maint: use "m_" prefix for member variables of classes associated with EIG.
Rik <rik@octave.org>
parents:
29358
diff
changeset
|
383 m_lambda = ComplexColumnVector (wr); |
1f34286a2637
maint: use "m_" prefix for member variables of classes associated with EIG.
Rik <rik@octave.org>
parents:
29358
diff
changeset
|
384 m_v = (calc_rev ? ComplexMatrix (atmp) : ComplexMatrix ()); |
1f34286a2637
maint: use "m_" prefix for member variables of classes associated with EIG.
Rik <rik@octave.org>
parents:
29358
diff
changeset
|
385 m_w = (calc_lev ? ComplexMatrix (atmp) : ComplexMatrix ()); |
462 | 386 |
387 return info; | |
388 } | |
389 | |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
390 octave_idx_type |
22305
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
391 EIG::init (const Matrix& a, const Matrix& b, bool calc_rev, bool calc_lev, |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
392 bool force_qz) |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
393 { |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
394 if (a.any_element_is_inf_or_nan () || b.any_element_is_inf_or_nan ()) |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
395 (*current_liboctave_error_handler) |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
396 ("EIG: matrix contains Inf or NaN values"); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
397 |
22988
cd33c785e80e
put to_f77_int inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
22953
diff
changeset
|
398 F77_INT n = octave::to_f77_int (a.rows ()); |
cd33c785e80e
put to_f77_int inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
22953
diff
changeset
|
399 F77_INT nb = octave::to_f77_int (b.rows ()); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
400 |
22988
cd33c785e80e
put to_f77_int inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
22953
diff
changeset
|
401 F77_INT a_nc = octave::to_f77_int (a.cols ()); |
cd33c785e80e
put to_f77_int inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
22953
diff
changeset
|
402 F77_INT b_nc = octave::to_f77_int (b.cols ()); |
23085
dffa2b8c9482
maint: strip trailing whitespace from source files.
John W. Eaton <jwe@octave.org>
parents:
23084
diff
changeset
|
403 |
22953
fd649fd3db75
use F77_INT instead of octave_idx_type for liboctave EIG classes
John W. Eaton <jwe@octave.org>
parents:
22755
diff
changeset
|
404 if (n != a_nc || nb != b_nc) |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
405 (*current_liboctave_error_handler) ("EIG requires square matrix"); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
406 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
407 if (n != nb) |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
408 (*current_liboctave_error_handler) ("EIG requires same size matrices"); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
409 |
22953
fd649fd3db75
use F77_INT instead of octave_idx_type for liboctave EIG classes
John W. Eaton <jwe@octave.org>
parents:
22755
diff
changeset
|
410 F77_INT info = 0; |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
411 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
412 Matrix tmp = b; |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
413 double *tmp_data = tmp.fortran_vec (); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
414 |
22305
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
415 if (! force_qz) |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
416 { |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
417 F77_XFCN (dpotrf, DPOTRF, (F77_CONST_CHAR_ARG2 ("L", 1), |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
418 n, tmp_data, n, |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
419 info |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
420 F77_CHAR_ARG_LEN (1))); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
421 |
23596
b63c3a09aee7
maint: Deprecate is_symmetric and replace with issymmetric.
Rik <rik@octave.org>
parents:
23594
diff
changeset
|
422 if (a.issymmetric () && b.issymmetric () && info == 0) |
22305
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
423 return symmetric_init (a, b, calc_rev, calc_lev); |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
424 } |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
425 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
426 Matrix atmp = a; |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
427 double *atmp_data = atmp.fortran_vec (); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
428 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
429 Matrix btmp = b; |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
430 double *btmp_data = btmp.fortran_vec (); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
431 |
11570
57632dea2446
attempt better backward compatibility for Array constructors
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
432 Array<double> ar (dim_vector (n, 1)); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
433 double *par = ar.fortran_vec (); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
434 |
11570
57632dea2446
attempt better backward compatibility for Array constructors
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
435 Array<double> ai (dim_vector (n, 1)); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
436 double *pai = ai.fortran_vec (); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
437 |
11570
57632dea2446
attempt better backward compatibility for Array constructors
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
438 Array<double> beta (dim_vector (n, 1)); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
439 double *pbeta = beta.fortran_vec (); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
440 |
23450
855122b993da
maint: Wrap tertiary operator in parentheses "(COND ? x : y)".
Rik <rik@octave.org>
parents:
23220
diff
changeset
|
441 F77_INT tnvr = (calc_rev ? n : 0); |
9020
728e7943752d
EIG.cc: avoid volatile decl for tmp variable
John W. Eaton <jwe@octave.org>
parents:
8920
diff
changeset
|
442 Matrix vr (tnvr, tnvr); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
443 double *pvr = vr.fortran_vec (); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
444 |
23450
855122b993da
maint: Wrap tertiary operator in parentheses "(COND ? x : y)".
Rik <rik@octave.org>
parents:
23220
diff
changeset
|
445 F77_INT tnvl = (calc_lev ? n : 0); |
22305
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
446 Matrix vl (tnvl, tnvl); |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
447 double *pvl = vl.fortran_vec (); |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
448 |
22953
fd649fd3db75
use F77_INT instead of octave_idx_type for liboctave EIG classes
John W. Eaton <jwe@octave.org>
parents:
22755
diff
changeset
|
449 F77_INT lwork = -1; |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
450 double dummy_work; |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
451 |
22305
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
452 F77_XFCN (dggev, DGGEV, (F77_CONST_CHAR_ARG2 (calc_lev ? "V" : "N", 1), |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
453 F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1), |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11570
diff
changeset
|
454 n, atmp_data, n, btmp_data, n, |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
455 par, pai, pbeta, |
22305
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
456 pvl, n, pvr, n, |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
457 &dummy_work, lwork, info |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
458 F77_CHAR_ARG_LEN (1) |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
459 F77_CHAR_ARG_LEN (1))); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
460 |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
461 if (info != 0) |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
462 (*current_liboctave_error_handler) ("dggev workspace query failed"); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
463 |
22953
fd649fd3db75
use F77_INT instead of octave_idx_type for liboctave EIG classes
John W. Eaton <jwe@octave.org>
parents:
22755
diff
changeset
|
464 lwork = static_cast<F77_INT> (dummy_work); |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
465 Array<double> work (dim_vector (lwork, 1)); |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
466 double *pwork = work.fortran_vec (); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
467 |
22305
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
468 F77_XFCN (dggev, DGGEV, (F77_CONST_CHAR_ARG2 (calc_lev ? "V" : "N", 1), |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
469 F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1), |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
470 n, atmp_data, n, btmp_data, n, |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
471 par, pai, pbeta, |
22305
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
472 pvl, n, pvr, n, |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
473 pwork, lwork, info |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
474 F77_CHAR_ARG_LEN (1) |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
475 F77_CHAR_ARG_LEN (1))); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
476 |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
477 if (info < 0) |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
478 (*current_liboctave_error_handler) ("unrecoverable error in dggev"); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
479 |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
480 if (info > 0) |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
481 (*current_liboctave_error_handler) ("dggev failed to converge"); |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
482 |
30044
1f34286a2637
maint: use "m_" prefix for member variables of classes associated with EIG.
Rik <rik@octave.org>
parents:
29358
diff
changeset
|
483 m_lambda.resize (n); |
23450
855122b993da
maint: Wrap tertiary operator in parentheses "(COND ? x : y)".
Rik <rik@octave.org>
parents:
23220
diff
changeset
|
484 F77_INT nvr = (calc_rev ? n : 0); |
30044
1f34286a2637
maint: use "m_" prefix for member variables of classes associated with EIG.
Rik <rik@octave.org>
parents:
29358
diff
changeset
|
485 m_v.resize (nvr, nvr); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
486 |
23450
855122b993da
maint: Wrap tertiary operator in parentheses "(COND ? x : y)".
Rik <rik@octave.org>
parents:
23220
diff
changeset
|
487 F77_INT nvl = (calc_lev ? n : 0); |
30044
1f34286a2637
maint: use "m_" prefix for member variables of classes associated with EIG.
Rik <rik@octave.org>
parents:
29358
diff
changeset
|
488 m_w.resize (nvl, nvl); |
22305
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
489 |
22953
fd649fd3db75
use F77_INT instead of octave_idx_type for liboctave EIG classes
John W. Eaton <jwe@octave.org>
parents:
22755
diff
changeset
|
490 for (F77_INT j = 0; j < n; j++) |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
491 { |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
492 if (ai.elem (j) == 0.0) |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
493 { |
30044
1f34286a2637
maint: use "m_" prefix for member variables of classes associated with EIG.
Rik <rik@octave.org>
parents:
29358
diff
changeset
|
494 m_lambda.elem (j) = Complex (ar.elem (j) / beta.elem (j)); |
22953
fd649fd3db75
use F77_INT instead of octave_idx_type for liboctave EIG classes
John W. Eaton <jwe@octave.org>
parents:
22755
diff
changeset
|
495 for (F77_INT i = 0; i < nvr; i++) |
30044
1f34286a2637
maint: use "m_" prefix for member variables of classes associated with EIG.
Rik <rik@octave.org>
parents:
29358
diff
changeset
|
496 m_v.elem (i, j) = vr.elem (i, j); |
22953
fd649fd3db75
use F77_INT instead of octave_idx_type for liboctave EIG classes
John W. Eaton <jwe@octave.org>
parents:
22755
diff
changeset
|
497 for (F77_INT i = 0; i < nvl; i++) |
30044
1f34286a2637
maint: use "m_" prefix for member variables of classes associated with EIG.
Rik <rik@octave.org>
parents:
29358
diff
changeset
|
498 m_w.elem (i, j) = vl.elem (i, j); |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
499 } |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
500 else |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
501 { |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
502 if (j+1 >= n) |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
503 (*current_liboctave_error_handler) ("EIG: internal error"); |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
504 |
30044
1f34286a2637
maint: use "m_" prefix for member variables of classes associated with EIG.
Rik <rik@octave.org>
parents:
29358
diff
changeset
|
505 m_lambda.elem (j) = Complex (ar.elem (j) / beta.elem (j), |
1f34286a2637
maint: use "m_" prefix for member variables of classes associated with EIG.
Rik <rik@octave.org>
parents:
29358
diff
changeset
|
506 ai.elem (j) / beta.elem (j)); |
1f34286a2637
maint: use "m_" prefix for member variables of classes associated with EIG.
Rik <rik@octave.org>
parents:
29358
diff
changeset
|
507 m_lambda.elem (j+1) = Complex (ar.elem (j+1) / beta.elem (j+1), |
1f34286a2637
maint: use "m_" prefix for member variables of classes associated with EIG.
Rik <rik@octave.org>
parents:
29358
diff
changeset
|
508 ai.elem (j+1) / beta.elem (j+1)); |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
509 |
22953
fd649fd3db75
use F77_INT instead of octave_idx_type for liboctave EIG classes
John W. Eaton <jwe@octave.org>
parents:
22755
diff
changeset
|
510 for (F77_INT i = 0; i < nvr; i++) |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
511 { |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
512 double real_part = vr.elem (i, j); |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
513 double imag_part = vr.elem (i, j+1); |
30044
1f34286a2637
maint: use "m_" prefix for member variables of classes associated with EIG.
Rik <rik@octave.org>
parents:
29358
diff
changeset
|
514 m_v.elem (i, j) = Complex (real_part, imag_part); |
1f34286a2637
maint: use "m_" prefix for member variables of classes associated with EIG.
Rik <rik@octave.org>
parents:
29358
diff
changeset
|
515 m_v.elem (i, j+1) = Complex (real_part, -imag_part); |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
516 } |
22953
fd649fd3db75
use F77_INT instead of octave_idx_type for liboctave EIG classes
John W. Eaton <jwe@octave.org>
parents:
22755
diff
changeset
|
517 for (F77_INT i = 0; i < nvl; i++) |
22305
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
518 { |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
519 double real_part = vl.elem (i, j); |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
520 double imag_part = vl.elem (i, j+1); |
30044
1f34286a2637
maint: use "m_" prefix for member variables of classes associated with EIG.
Rik <rik@octave.org>
parents:
29358
diff
changeset
|
521 m_w.elem (i, j) = Complex (real_part, imag_part); |
1f34286a2637
maint: use "m_" prefix for member variables of classes associated with EIG.
Rik <rik@octave.org>
parents:
29358
diff
changeset
|
522 m_w.elem (i, j+1) = Complex (real_part, -imag_part); |
22305
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
523 } |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
524 j++; |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
525 } |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
526 } |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
527 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
528 return info; |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
529 } |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
530 |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11570
diff
changeset
|
531 octave_idx_type |
22305
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
532 EIG::symmetric_init (const Matrix& a, const Matrix& b, bool calc_rev, |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
533 bool calc_lev) |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
534 { |
22988
cd33c785e80e
put to_f77_int inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
22953
diff
changeset
|
535 F77_INT n = octave::to_f77_int (a.rows ()); |
cd33c785e80e
put to_f77_int inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
22953
diff
changeset
|
536 F77_INT nb = octave::to_f77_int (b.rows ()); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
537 |
22988
cd33c785e80e
put to_f77_int inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
22953
diff
changeset
|
538 F77_INT a_nc = octave::to_f77_int (a.cols ()); |
cd33c785e80e
put to_f77_int inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
22953
diff
changeset
|
539 F77_INT b_nc = octave::to_f77_int (b.cols ()); |
22953
fd649fd3db75
use F77_INT instead of octave_idx_type for liboctave EIG classes
John W. Eaton <jwe@octave.org>
parents:
22755
diff
changeset
|
540 |
fd649fd3db75
use F77_INT instead of octave_idx_type for liboctave EIG classes
John W. Eaton <jwe@octave.org>
parents:
22755
diff
changeset
|
541 if (n != a_nc || nb != b_nc) |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
542 (*current_liboctave_error_handler) ("EIG requires square matrix"); |
8339
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 (n != nb) |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
545 (*current_liboctave_error_handler) ("EIG requires same size matrices"); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
546 |
22953
fd649fd3db75
use F77_INT instead of octave_idx_type for liboctave EIG classes
John W. Eaton <jwe@octave.org>
parents:
22755
diff
changeset
|
547 F77_INT info = 0; |
8339
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 Matrix atmp = a; |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
550 double *atmp_data = atmp.fortran_vec (); |
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 Matrix btmp = b; |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
553 double *btmp_data = btmp.fortran_vec (); |
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 ColumnVector wr (n); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
556 double *pwr = wr.fortran_vec (); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
557 |
22953
fd649fd3db75
use F77_INT instead of octave_idx_type for liboctave EIG classes
John W. Eaton <jwe@octave.org>
parents:
22755
diff
changeset
|
558 F77_INT lwork = -1; |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
559 double dummy_work; |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
560 |
22305
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
561 F77_XFCN (dsygv, DSYGV, (1, F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1), |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
562 F77_CONST_CHAR_ARG2 ("U", 1), |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11570
diff
changeset
|
563 n, atmp_data, n, |
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11570
diff
changeset
|
564 btmp_data, n, |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
565 pwr, &dummy_work, lwork, info |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
566 F77_CHAR_ARG_LEN (1) |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
567 F77_CHAR_ARG_LEN (1))); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
568 |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
569 if (info != 0) |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
570 (*current_liboctave_error_handler) ("dsygv workspace query failed"); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
571 |
22953
fd649fd3db75
use F77_INT instead of octave_idx_type for liboctave EIG classes
John W. Eaton <jwe@octave.org>
parents:
22755
diff
changeset
|
572 lwork = static_cast<F77_INT> (dummy_work); |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
573 Array<double> work (dim_vector (lwork, 1)); |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
574 double *pwork = work.fortran_vec (); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
575 |
22305
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
576 F77_XFCN (dsygv, DSYGV, (1, F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1), |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
577 F77_CONST_CHAR_ARG2 ("U", 1), |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
578 n, atmp_data, n, |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
579 btmp_data, n, |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
580 pwr, pwork, lwork, info |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
581 F77_CHAR_ARG_LEN (1) |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
582 F77_CHAR_ARG_LEN (1))); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
583 |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
584 if (info < 0) |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
585 (*current_liboctave_error_handler) ("unrecoverable error in dsygv"); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
586 |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
587 if (info > 0) |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
588 (*current_liboctave_error_handler) ("dsygv failed to converge"); |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
589 |
30044
1f34286a2637
maint: use "m_" prefix for member variables of classes associated with EIG.
Rik <rik@octave.org>
parents:
29358
diff
changeset
|
590 m_lambda = ComplexColumnVector (wr); |
1f34286a2637
maint: use "m_" prefix for member variables of classes associated with EIG.
Rik <rik@octave.org>
parents:
29358
diff
changeset
|
591 m_v = (calc_rev ? ComplexMatrix (atmp) : ComplexMatrix ()); |
1f34286a2637
maint: use "m_" prefix for member variables of classes associated with EIG.
Rik <rik@octave.org>
parents:
29358
diff
changeset
|
592 m_w = (calc_lev ? ComplexMatrix (atmp) : ComplexMatrix ()); |
8339
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 |
22305
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
598 EIG::init (const ComplexMatrix& a, const ComplexMatrix& b, bool calc_rev, |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
599 bool calc_lev, bool force_qz) |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
600 { |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
601 if (a.any_element_is_inf_or_nan () || b.any_element_is_inf_or_nan ()) |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
602 (*current_liboctave_error_handler) |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
603 ("EIG: matrix contains Inf or NaN values"); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
604 |
22988
cd33c785e80e
put to_f77_int inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
22953
diff
changeset
|
605 F77_INT n = octave::to_f77_int (a.rows ()); |
cd33c785e80e
put to_f77_int inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
22953
diff
changeset
|
606 F77_INT nb = octave::to_f77_int (b.rows ()); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
607 |
22988
cd33c785e80e
put to_f77_int inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
22953
diff
changeset
|
608 F77_INT a_nc = octave::to_f77_int (a.cols ()); |
cd33c785e80e
put to_f77_int inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
22953
diff
changeset
|
609 F77_INT b_nc = octave::to_f77_int (b.cols ()); |
22953
fd649fd3db75
use F77_INT instead of octave_idx_type for liboctave EIG classes
John W. Eaton <jwe@octave.org>
parents:
22755
diff
changeset
|
610 |
fd649fd3db75
use F77_INT instead of octave_idx_type for liboctave EIG classes
John W. Eaton <jwe@octave.org>
parents:
22755
diff
changeset
|
611 if (n != a_nc || nb != b_nc) |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
612 (*current_liboctave_error_handler) ("EIG requires square matrix"); |
8339
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 if (n != nb) |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
615 (*current_liboctave_error_handler) ("EIG requires same size matrices"); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
616 |
22953
fd649fd3db75
use F77_INT instead of octave_idx_type for liboctave EIG classes
John W. Eaton <jwe@octave.org>
parents:
22755
diff
changeset
|
617 F77_INT info = 0; |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
618 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
619 ComplexMatrix tmp = b; |
30394
f3f3e3793fb5
maint: style check C++ files in liboctave/ ahead of 7.1 release.
Rik <rik@octave.org>
parents:
30044
diff
changeset
|
620 Complex *tmp_data = tmp.fortran_vec (); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
621 |
22305
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
622 if (! force_qz) |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
623 { |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
624 F77_XFCN (zpotrf, ZPOTRF, (F77_CONST_CHAR_ARG2 ("L", 1), |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
625 n, F77_DBLE_CMPLX_ARG (tmp_data), n, |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
626 info |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
627 F77_CHAR_ARG_LEN (1))); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
628 |
23594
af5b813503cb
maint: Deprecate is_hermitian and replace with ishermitian.
Rik <rik@octave.org>
parents:
23475
diff
changeset
|
629 if (a.ishermitian () && b.ishermitian () && info == 0) |
22305
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
630 return hermitian_init (a, b, calc_rev, calc_lev); |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
631 } |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
632 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
633 ComplexMatrix atmp = a; |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
634 Complex *atmp_data = atmp.fortran_vec (); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
635 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
636 ComplexMatrix btmp = b; |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
637 Complex *btmp_data = btmp.fortran_vec (); |
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 ComplexColumnVector alpha (n); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
640 Complex *palpha = alpha.fortran_vec (); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
641 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
642 ComplexColumnVector beta (n); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
643 Complex *pbeta = beta.fortran_vec (); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
644 |
23450
855122b993da
maint: Wrap tertiary operator in parentheses "(COND ? x : y)".
Rik <rik@octave.org>
parents:
23220
diff
changeset
|
645 F77_INT nvr = (calc_rev ? n : 0); |
22305
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
646 ComplexMatrix vrtmp (nvr, nvr); |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
647 Complex *pvr = vrtmp.fortran_vec (); |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
648 |
23450
855122b993da
maint: Wrap tertiary operator in parentheses "(COND ? x : y)".
Rik <rik@octave.org>
parents:
23220
diff
changeset
|
649 F77_INT nvl = (calc_lev ? n : 0); |
22305
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
650 ComplexMatrix vltmp (nvl, nvl); |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
651 Complex *pvl = vltmp.fortran_vec (); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
652 |
22953
fd649fd3db75
use F77_INT instead of octave_idx_type for liboctave EIG classes
John W. Eaton <jwe@octave.org>
parents:
22755
diff
changeset
|
653 F77_INT lwork = -1; |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
654 Complex dummy_work; |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
655 |
22953
fd649fd3db75
use F77_INT instead of octave_idx_type for liboctave EIG classes
John W. Eaton <jwe@octave.org>
parents:
22755
diff
changeset
|
656 F77_INT lrwork = 8*n; |
11570
57632dea2446
attempt better backward compatibility for Array constructors
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
657 Array<double> rwork (dim_vector (lrwork, 1)); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
658 double *prwork = rwork.fortran_vec (); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
659 |
27009
fc73dafece57
Fix return of left-handed vectors when inputs are complex (bug #56026)
Rik <rik@octave.org>
parents:
26376
diff
changeset
|
660 F77_XFCN (zggev, ZGGEV, (F77_CONST_CHAR_ARG2 (calc_lev ? "V" : "N", 1), |
22305
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
661 F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1), |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
662 n, F77_DBLE_CMPLX_ARG (atmp_data), n, |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
663 F77_DBLE_CMPLX_ARG (btmp_data), n, |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
664 F77_DBLE_CMPLX_ARG (palpha), |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
665 F77_DBLE_CMPLX_ARG (pbeta), |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
666 F77_DBLE_CMPLX_ARG (pvl), n, |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
667 F77_DBLE_CMPLX_ARG (pvr), n, |
23475
d691ed308237
maint: Clean up #includes in liboctave/numeric directory.
Rik <rik@octave.org>
parents:
23450
diff
changeset
|
668 F77_DBLE_CMPLX_ARG (&dummy_work), lwork, prwork, |
d691ed308237
maint: Clean up #includes in liboctave/numeric directory.
Rik <rik@octave.org>
parents:
23450
diff
changeset
|
669 info |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
670 F77_CHAR_ARG_LEN (1) |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
671 F77_CHAR_ARG_LEN (1))); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
672 |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
673 if (info != 0) |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
674 (*current_liboctave_error_handler) ("zggev workspace query failed"); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
675 |
22953
fd649fd3db75
use F77_INT instead of octave_idx_type for liboctave EIG classes
John W. Eaton <jwe@octave.org>
parents:
22755
diff
changeset
|
676 lwork = static_cast<F77_INT> (dummy_work.real ()); |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
677 Array<Complex> work (dim_vector (lwork, 1)); |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
678 Complex *pwork = work.fortran_vec (); |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11570
diff
changeset
|
679 |
27009
fc73dafece57
Fix return of left-handed vectors when inputs are complex (bug #56026)
Rik <rik@octave.org>
parents:
26376
diff
changeset
|
680 F77_XFCN (zggev, ZGGEV, (F77_CONST_CHAR_ARG2 (calc_lev ? "V" : "N", 1), |
22305
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
681 F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1), |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
682 n, F77_DBLE_CMPLX_ARG (atmp_data), n, |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
683 F77_DBLE_CMPLX_ARG (btmp_data), n, |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
684 F77_DBLE_CMPLX_ARG (palpha), |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
685 F77_DBLE_CMPLX_ARG (pbeta), |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
686 F77_DBLE_CMPLX_ARG (pvl), n, |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
687 F77_DBLE_CMPLX_ARG (pvr), n, |
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
688 F77_DBLE_CMPLX_ARG (pwork), lwork, prwork, info |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
689 F77_CHAR_ARG_LEN (1) |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
690 F77_CHAR_ARG_LEN (1))); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
691 |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
692 if (info < 0) |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
693 (*current_liboctave_error_handler) ("unrecoverable error in zggev"); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
694 |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
695 if (info > 0) |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
696 (*current_liboctave_error_handler) ("zggev failed to converge"); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
697 |
30044
1f34286a2637
maint: use "m_" prefix for member variables of classes associated with EIG.
Rik <rik@octave.org>
parents:
29358
diff
changeset
|
698 m_lambda.resize (n); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
699 |
22953
fd649fd3db75
use F77_INT instead of octave_idx_type for liboctave EIG classes
John W. Eaton <jwe@octave.org>
parents:
22755
diff
changeset
|
700 for (F77_INT j = 0; j < n; j++) |
30044
1f34286a2637
maint: use "m_" prefix for member variables of classes associated with EIG.
Rik <rik@octave.org>
parents:
29358
diff
changeset
|
701 m_lambda.elem (j) = alpha.elem (j) / beta.elem (j); |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
702 |
30044
1f34286a2637
maint: use "m_" prefix for member variables of classes associated with EIG.
Rik <rik@octave.org>
parents:
29358
diff
changeset
|
703 m_v = vrtmp; |
1f34286a2637
maint: use "m_" prefix for member variables of classes associated with EIG.
Rik <rik@octave.org>
parents:
29358
diff
changeset
|
704 m_w = vltmp; |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
705 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
706 return info; |
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 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
709 octave_idx_type |
17769
49a5a4be04a1
maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
710 EIG::hermitian_init (const ComplexMatrix& a, const ComplexMatrix& b, |
22305
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
711 bool calc_rev, bool calc_lev) |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
712 { |
22988
cd33c785e80e
put to_f77_int inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
22953
diff
changeset
|
713 F77_INT n = octave::to_f77_int (a.rows ()); |
cd33c785e80e
put to_f77_int inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
22953
diff
changeset
|
714 F77_INT nb = octave::to_f77_int (b.rows ()); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
715 |
22988
cd33c785e80e
put to_f77_int inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
22953
diff
changeset
|
716 F77_INT a_nc = octave::to_f77_int (a.cols ()); |
cd33c785e80e
put to_f77_int inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
22953
diff
changeset
|
717 F77_INT b_nc = octave::to_f77_int (b.cols ()); |
22953
fd649fd3db75
use F77_INT instead of octave_idx_type for liboctave EIG classes
John W. Eaton <jwe@octave.org>
parents:
22755
diff
changeset
|
718 |
fd649fd3db75
use F77_INT instead of octave_idx_type for liboctave EIG classes
John W. Eaton <jwe@octave.org>
parents:
22755
diff
changeset
|
719 if (n != a_nc || nb != b_nc) |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
720 (*current_liboctave_error_handler) ("EIG requires square matrix"); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
721 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
722 if (n != nb) |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
723 (*current_liboctave_error_handler) ("EIG requires same size matrices"); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
724 |
22953
fd649fd3db75
use F77_INT instead of octave_idx_type for liboctave EIG classes
John W. Eaton <jwe@octave.org>
parents:
22755
diff
changeset
|
725 F77_INT info = 0; |
8339
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 ComplexMatrix atmp = a; |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
728 Complex *atmp_data = atmp.fortran_vec (); |
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 ComplexMatrix btmp = b; |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
731 Complex *btmp_data = btmp.fortran_vec (); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
732 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
733 ColumnVector wr (n); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
734 double *pwr = wr.fortran_vec (); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
735 |
22953
fd649fd3db75
use F77_INT instead of octave_idx_type for liboctave EIG classes
John W. Eaton <jwe@octave.org>
parents:
22755
diff
changeset
|
736 F77_INT lwork = -1; |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
737 Complex dummy_work; |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
738 |
22953
fd649fd3db75
use F77_INT instead of octave_idx_type for liboctave EIG classes
John W. Eaton <jwe@octave.org>
parents:
22755
diff
changeset
|
739 F77_INT lrwork = 3*n; |
11570
57632dea2446
attempt better backward compatibility for Array constructors
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
740 Array<double> rwork (dim_vector (lrwork, 1)); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
741 double *prwork = rwork.fortran_vec (); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
742 |
22305
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
743 F77_XFCN (zhegv, ZHEGV, (1, F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1), |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
744 F77_CONST_CHAR_ARG2 ("U", 1), |
22135
407c66ae1e20
reduce warnings from GCC's link-time optimization feature (bug #48531)
John W. Eaton <jwe@octave.org>
parents:
22133
diff
changeset
|
745 n, F77_DBLE_CMPLX_ARG (atmp_data), n, |
407c66ae1e20
reduce warnings from GCC's link-time optimization feature (bug #48531)
John W. Eaton <jwe@octave.org>
parents:
22133
diff
changeset
|
746 F77_DBLE_CMPLX_ARG (btmp_data), n, |
407c66ae1e20
reduce warnings from GCC's link-time optimization feature (bug #48531)
John W. Eaton <jwe@octave.org>
parents:
22133
diff
changeset
|
747 pwr, F77_DBLE_CMPLX_ARG (&dummy_work), lwork, |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
748 prwork, info |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
749 F77_CHAR_ARG_LEN (1) |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
750 F77_CHAR_ARG_LEN (1))); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
751 |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
752 if (info != 0) |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
753 (*current_liboctave_error_handler) ("zhegv workspace query failed"); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
754 |
22953
fd649fd3db75
use F77_INT instead of octave_idx_type for liboctave EIG classes
John W. Eaton <jwe@octave.org>
parents:
22755
diff
changeset
|
755 lwork = static_cast<F77_INT> (dummy_work.real ()); |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
756 Array<Complex> work (dim_vector (lwork, 1)); |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
757 Complex *pwork = work.fortran_vec (); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
758 |
22305
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
22135
diff
changeset
|
759 F77_XFCN (zhegv, ZHEGV, (1, F77_CONST_CHAR_ARG2 (calc_rev ? "V" : "N", 1), |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
760 F77_CONST_CHAR_ARG2 ("U", 1), |
22135
407c66ae1e20
reduce warnings from GCC's link-time optimization feature (bug #48531)
John W. Eaton <jwe@octave.org>
parents:
22133
diff
changeset
|
761 n, F77_DBLE_CMPLX_ARG (atmp_data), n, |
407c66ae1e20
reduce warnings from GCC's link-time optimization feature (bug #48531)
John W. Eaton <jwe@octave.org>
parents:
22133
diff
changeset
|
762 F77_DBLE_CMPLX_ARG (btmp_data), n, |
407c66ae1e20
reduce warnings from GCC's link-time optimization feature (bug #48531)
John W. Eaton <jwe@octave.org>
parents:
22133
diff
changeset
|
763 pwr, F77_DBLE_CMPLX_ARG (pwork), lwork, prwork, info |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
764 F77_CHAR_ARG_LEN (1) |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
765 F77_CHAR_ARG_LEN (1))); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
766 |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
767 if (info < 0) |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
768 (*current_liboctave_error_handler) ("unrecoverable error in zhegv"); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
769 |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
770 if (info > 0) |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
771 (*current_liboctave_error_handler) ("zhegv failed to converge"); |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
772 |
30044
1f34286a2637
maint: use "m_" prefix for member variables of classes associated with EIG.
Rik <rik@octave.org>
parents:
29358
diff
changeset
|
773 m_lambda = ComplexColumnVector (wr); |
1f34286a2637
maint: use "m_" prefix for member variables of classes associated with EIG.
Rik <rik@octave.org>
parents:
29358
diff
changeset
|
774 m_v = (calc_rev ? ComplexMatrix (atmp) : ComplexMatrix ()); |
1f34286a2637
maint: use "m_" prefix for member variables of classes associated with EIG.
Rik <rik@octave.org>
parents:
29358
diff
changeset
|
775 m_w = (calc_lev ? ComplexMatrix (atmp) : ComplexMatrix ()); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
776 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
777 return info; |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
778 } |