Mercurial > octave
annotate liboctave/numeric/EIG.cc @ 21136:7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Remove statements after call to handler that are no longer reachable.
Place input validation first and immediately call handler if necessary.
Change if/error_handler/else to if/error_handler and re-indent code.
* Array-util.cc, Array.cc, CColVector.cc, CDiagMatrix.cc, CMatrix.cc,
CNDArray.cc, CRowVector.cc, CSparse.cc, DiagArray2.cc, MArray.cc,
PermMatrix.cc, Sparse.cc, Sparse.h, chMatrix.cc, chNDArray.cc, dColVector.cc,
dDiagMatrix.cc, dMatrix.cc, dNDArray.cc, dRowVector.cc, dSparse.cc,
fCColVector.cc, fCDiagMatrix.cc, fCMatrix.cc, fCNDArray.cc, fCRowVector.cc,
fColVector.cc, fDiagMatrix.cc, fMatrix.cc, fNDArray.cc, fRowVector.cc,
idx-vector.cc, CmplxAEPBAL.cc, CmplxCHOL.cc, CmplxGEPBAL.cc, CmplxHESS.cc,
CmplxLU.cc, CmplxQR.cc, CmplxSCHUR.cc, CmplxSVD.cc, DASPK.cc, EIG.cc, LSODE.cc,
Quad.cc, SparseCmplxCHOL.cc, SparseCmplxLU.cc, SparseCmplxQR.cc, SparseQR.cc,
SparsedbleCHOL.cc, SparsedbleLU.cc, base-lu.cc, bsxfun-defs.cc, dbleAEPBAL.cc,
dbleCHOL.cc, dbleGEPBAL.cc, dbleHESS.cc, dbleLU.cc, dbleQR.cc, dbleSCHUR.cc,
dbleSVD.cc, eigs-base.cc, fCmplxAEPBAL.cc, fCmplxCHOL.cc, fCmplxLU.cc,
fCmplxQR.cc, fCmplxSCHUR.cc, fEIG.cc, floatAEPBAL.cc, floatCHOL.cc,
floatGEPBAL.cc, floatHESS.cc, floatLU.cc, floatQR.cc, floatSCHUR.cc,
floatSVD.cc, lo-specfun.cc, oct-fftw.cc, oct-rand.cc, oct-spparms.cc,
sparse-base-chol.cc, sparse-dmsolve.cc, file-ops.cc, lo-sysdep.cc,
mach-info.cc, oct-env.cc, oct-syscalls.cc, cmd-edit.cc, cmd-hist.cc,
data-conv.cc, lo-ieee.cc, lo-regexp.cc, oct-base64.cc, oct-shlib.cc,
pathsearch.cc, singleton-cleanup.cc, sparse-util.cc, unwind-prot.cc:
Remove statements after call to handler that are no longer reachable.
Place input validation first and immediately call handler if necessary.
Change if/error_handler/else to if/error_handler and re-indent code.
author | Rik <rik@octave.org> |
---|---|
date | Sat, 23 Jan 2016 13:52:03 -0800 |
parents | ff904ae0285b |
children | f7121e111991 |
rev | line source |
---|---|
462 | 1 /* |
2 | |
19697
4197fc428c7d
maint: Update copyright notices for 2015.
John W. Eaton <jwe@octave.org>
parents:
17769
diff
changeset
|
3 Copyright (C) 1994-2015 John W. Eaton |
462 | 4 |
5 This file is part of Octave. | |
6 | |
7 Octave is free software; you can redistribute it and/or modify it | |
8 under the terms of the GNU General Public License as published by the | |
7016 | 9 Free Software Foundation; either version 3 of the License, or (at your |
10 option) any later version. | |
462 | 11 |
12 Octave is distributed in the hope that it will be useful, but WITHOUT | |
13 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or | |
14 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License | |
15 for more details. | |
16 | |
17 You should have received a copy of the GNU General Public License | |
7016 | 18 along with Octave; see the file COPYING. If not, see |
19 <http://www.gnu.org/licenses/>. | |
462 | 20 |
21 */ | |
22 | |
23 #ifdef HAVE_CONFIG_H | |
1192 | 24 #include <config.h> |
462 | 25 #endif |
26 | |
27 #include "EIG.h" | |
2815 | 28 #include "dColVector.h" |
1847 | 29 #include "f77-fcn.h" |
462 | 30 #include "lo-error.h" |
31 | |
32 extern "C" | |
33 { | |
4552 | 34 F77_RET_T |
35 F77_FUNC (dgeev, DGEEV) (F77_CONST_CHAR_ARG_DECL, | |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
36 F77_CONST_CHAR_ARG_DECL, |
11495 | 37 const octave_idx_type&, double*, |
38 const octave_idx_type&, double*, double*, | |
39 double*, const octave_idx_type&, double*, | |
40 const octave_idx_type&, double*, | |
41 const octave_idx_type&, octave_idx_type& | |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
42 F77_CHAR_ARG_LEN_DECL |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
43 F77_CHAR_ARG_LEN_DECL); |
462 | 44 |
4552 | 45 F77_RET_T |
46 F77_FUNC (zgeev, ZGEEV) (F77_CONST_CHAR_ARG_DECL, | |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
47 F77_CONST_CHAR_ARG_DECL, |
11495 | 48 const octave_idx_type&, Complex*, |
49 const octave_idx_type&, Complex*, | |
50 Complex*, const octave_idx_type&, Complex*, | |
51 const octave_idx_type&, Complex*, | |
52 const octave_idx_type&, double*, octave_idx_type& | |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
53 F77_CHAR_ARG_LEN_DECL |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
54 F77_CHAR_ARG_LEN_DECL); |
2815 | 55 |
4552 | 56 F77_RET_T |
57 F77_FUNC (dsyev, DSYEV) (F77_CONST_CHAR_ARG_DECL, | |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
58 F77_CONST_CHAR_ARG_DECL, |
11495 | 59 const octave_idx_type&, double*, |
60 const octave_idx_type&, double*, double*, | |
61 const octave_idx_type&, octave_idx_type& | |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
62 F77_CHAR_ARG_LEN_DECL |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
63 F77_CHAR_ARG_LEN_DECL); |
2815 | 64 |
4552 | 65 F77_RET_T |
66 F77_FUNC (zheev, ZHEEV) (F77_CONST_CHAR_ARG_DECL, | |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
67 F77_CONST_CHAR_ARG_DECL, |
11495 | 68 const octave_idx_type&, Complex*, |
69 const octave_idx_type&, double*, | |
70 Complex*, const octave_idx_type&, double*, | |
71 octave_idx_type& | |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
72 F77_CHAR_ARG_LEN_DECL |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
73 F77_CHAR_ARG_LEN_DECL); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
74 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
75 F77_RET_T |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11570
diff
changeset
|
76 F77_FUNC (dpotrf, DPOTRF) (F77_CONST_CHAR_ARG_DECL, |
11495 | 77 const octave_idx_type&, double*, |
78 const octave_idx_type&, octave_idx_type& | |
79 F77_CHAR_ARG_LEN_DECL | |
80 F77_CHAR_ARG_LEN_DECL); | |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
81 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
82 F77_RET_T |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11570
diff
changeset
|
83 F77_FUNC (zpotrf, ZPOTRF) (F77_CONST_CHAR_ARG_DECL, |
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11570
diff
changeset
|
84 const octave_idx_type&, |
11495 | 85 Complex*, const octave_idx_type&, |
86 octave_idx_type& | |
87 F77_CHAR_ARG_LEN_DECL | |
88 F77_CHAR_ARG_LEN_DECL); | |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
89 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
90 F77_RET_T |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11570
diff
changeset
|
91 F77_FUNC (dggev, DGGEV) (F77_CONST_CHAR_ARG_DECL, |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
92 F77_CONST_CHAR_ARG_DECL, |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11570
diff
changeset
|
93 const octave_idx_type&, |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
94 double*, const octave_idx_type&, |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
95 double*, const octave_idx_type&, |
11495 | 96 double*, double*, double *, double*, |
97 const octave_idx_type&, double*, | |
98 const octave_idx_type&, double*, | |
99 const octave_idx_type&, octave_idx_type& | |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
100 F77_CHAR_ARG_LEN_DECL |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
101 F77_CHAR_ARG_LEN_DECL); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
102 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
103 F77_RET_T |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
104 F77_FUNC (dsygv, DSYGV) (const octave_idx_type&, |
11518 | 105 F77_CONST_CHAR_ARG_DECL, |
106 F77_CONST_CHAR_ARG_DECL, | |
11495 | 107 const octave_idx_type&, double*, |
108 const octave_idx_type&, double*, | |
109 const octave_idx_type&, double*, double*, | |
110 const octave_idx_type&, octave_idx_type& | |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
111 F77_CHAR_ARG_LEN_DECL |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
112 F77_CHAR_ARG_LEN_DECL); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
113 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
114 F77_RET_T |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11570
diff
changeset
|
115 F77_FUNC (zggev, ZGGEV) (F77_CONST_CHAR_ARG_DECL, |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
116 F77_CONST_CHAR_ARG_DECL, |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11570
diff
changeset
|
117 const octave_idx_type&, |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
118 Complex*, const octave_idx_type&, |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
119 Complex*, const octave_idx_type&, |
11495 | 120 Complex*, Complex*, Complex*, |
121 const octave_idx_type&, Complex*, | |
122 const octave_idx_type&, Complex*, | |
123 const octave_idx_type&, double*, octave_idx_type& | |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
124 F77_CHAR_ARG_LEN_DECL |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
125 F77_CHAR_ARG_LEN_DECL); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
126 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
127 F77_RET_T |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
128 F77_FUNC (zhegv, ZHEGV) (const octave_idx_type&, |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11570
diff
changeset
|
129 F77_CONST_CHAR_ARG_DECL, |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
130 F77_CONST_CHAR_ARG_DECL, |
11495 | 131 const octave_idx_type&, Complex*, |
132 const octave_idx_type&, Complex*, | |
133 const octave_idx_type&, double*, Complex*, | |
134 const octave_idx_type&, double*, octave_idx_type& | |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
135 F77_CHAR_ARG_LEN_DECL |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
136 F77_CHAR_ARG_LEN_DECL); |
462 | 137 } |
138 | |
5275 | 139 octave_idx_type |
4725 | 140 EIG::init (const Matrix& a, bool calc_ev) |
462 | 141 { |
5822 | 142 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
|
143 (*current_liboctave_error_handler) |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
144 ("EIG: matrix contains Inf or NaN values"); |
5822 | 145 |
2815 | 146 if (a.is_symmetric ()) |
4725 | 147 return symmetric_init (a, calc_ev); |
2815 | 148 |
5275 | 149 octave_idx_type n = a.rows (); |
1934 | 150 |
151 if (n != a.cols ()) | |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
152 (*current_liboctave_error_handler) ("EIG requires square matrix"); |
462 | 153 |
5275 | 154 octave_idx_type info = 0; |
462 | 155 |
1934 | 156 Matrix atmp = a; |
157 double *tmp_data = atmp.fortran_vec (); | |
158 | |
11570
57632dea2446
attempt better backward compatibility for Array constructors
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
159 Array<double> wr (dim_vector (n, 1)); |
1934 | 160 double *pwr = wr.fortran_vec (); |
161 | |
11570
57632dea2446
attempt better backward compatibility for Array constructors
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
162 Array<double> wi (dim_vector (n, 1)); |
1934 | 163 double *pwi = wi.fortran_vec (); |
164 | |
9020
728e7943752d
EIG.cc: avoid volatile decl for tmp variable
John W. Eaton <jwe@octave.org>
parents:
8920
diff
changeset
|
165 octave_idx_type tnvr = calc_ev ? n : 0; |
728e7943752d
EIG.cc: avoid volatile decl for tmp variable
John W. Eaton <jwe@octave.org>
parents:
8920
diff
changeset
|
166 Matrix vr (tnvr, tnvr); |
462 | 167 double *pvr = vr.fortran_vec (); |
1934 | 168 |
5275 | 169 octave_idx_type lwork = -1; |
4800 | 170 double dummy_work; |
462 | 171 |
1365 | 172 double *dummy = 0; |
5275 | 173 octave_idx_type idummy = 1; |
462 | 174 |
4552 | 175 F77_XFCN (dgeev, DGEEV, (F77_CONST_CHAR_ARG2 ("N", 1), |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
176 F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
177 n, tmp_data, n, pwr, pwi, dummy, |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
178 idummy, pvr, n, &dummy_work, lwork, info |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
179 F77_CHAR_ARG_LEN (1) |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
180 F77_CHAR_ARG_LEN (1))); |
462 | 181 |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
182 if (info != 0) |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
183 (*current_liboctave_error_handler) ("dgeev workspace query failed"); |
4800 | 184 |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
185 lwork = static_cast<octave_idx_type> (dummy_work); |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
186 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
|
187 double *pwork = work.fortran_vec (); |
4800 | 188 |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
189 F77_XFCN (dgeev, DGEEV, (F77_CONST_CHAR_ARG2 ("N", 1), |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
190 F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
191 n, tmp_data, n, pwr, pwi, dummy, |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
192 idummy, pvr, n, pwork, lwork, info |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
193 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
|
194 F77_CHAR_ARG_LEN (1))); |
4800 | 195 |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
196 if (info < 0) |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
197 (*current_liboctave_error_handler) ("unrecoverable error in dgeev"); |
1934 | 198 |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
199 if (info > 0) |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
200 (*current_liboctave_error_handler) ("dgeev failed to converge"); |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
201 |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
202 lambda.resize (n); |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
203 octave_idx_type nvr = calc_ev ? n : 0; |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
204 v.resize (nvr, nvr); |
4800 | 205 |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
206 for (octave_idx_type j = 0; j < n; j++) |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
207 { |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
208 if (wi.elem (j) == 0.0) |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
209 { |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
210 lambda.elem (j) = Complex (wr.elem (j)); |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
211 for (octave_idx_type i = 0; i < nvr; i++) |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
212 v.elem (i, j) = vr.elem (i, j); |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
213 } |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
214 else |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
215 { |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
216 if (j+1 >= n) |
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) ("EIG: internal error"); |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
218 |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
219 lambda.elem (j) = Complex (wr.elem (j), wi.elem (j)); |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
220 lambda.elem (j+1) = Complex (wr.elem (j+1), wi.elem (j+1)); |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
221 |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
222 for (octave_idx_type i = 0; i < nvr; i++) |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
223 { |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
224 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
|
225 double imag_part = vr.elem (i, j+1); |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
226 v.elem (i, j) = Complex (real_part, imag_part); |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
227 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
|
228 } |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
229 j++; |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
230 } |
2815 | 231 } |
232 | |
233 return info; | |
234 } | |
235 | |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11570
diff
changeset
|
236 octave_idx_type |
4725 | 237 EIG::symmetric_init (const Matrix& a, bool calc_ev) |
2815 | 238 { |
5275 | 239 octave_idx_type n = a.rows (); |
1934 | 240 |
2815 | 241 if (n != a.cols ()) |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
242 (*current_liboctave_error_handler) ("EIG requires square matrix"); |
2815 | 243 |
5275 | 244 octave_idx_type info = 0; |
2815 | 245 |
246 Matrix atmp = a; | |
247 double *tmp_data = atmp.fortran_vec (); | |
248 | |
3585 | 249 ColumnVector wr (n); |
2815 | 250 double *pwr = wr.fortran_vec (); |
251 | |
5275 | 252 octave_idx_type lwork = -1; |
4800 | 253 double dummy_work; |
2815 | 254 |
4725 | 255 F77_XFCN (dsyev, DSYEV, (F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
256 F77_CONST_CHAR_ARG2 ("U", 1), |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
257 n, tmp_data, n, pwr, &dummy_work, lwork, info |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
258 F77_CHAR_ARG_LEN (1) |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
259 F77_CHAR_ARG_LEN (1))); |
2815 | 260 |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
261 if (info != 0) |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
262 (*current_liboctave_error_handler) ("dsyev workspace query failed"); |
4800 | 263 |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
264 lwork = static_cast<octave_idx_type> (dummy_work); |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
265 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
|
266 double *pwork = work.fortran_vec (); |
4800 | 267 |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
268 F77_XFCN (dsyev, DSYEV, (F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
269 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
|
270 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
|
271 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
|
272 F77_CHAR_ARG_LEN (1))); |
4800 | 273 |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
274 if (info < 0) |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
275 (*current_liboctave_error_handler) ("unrecoverable error in dsyev"); |
4800 | 276 |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
277 if (info > 0) |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
278 (*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
|
279 |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
280 lambda = ComplexColumnVector (wr); |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
281 v = calc_ev ? ComplexMatrix (atmp) : ComplexMatrix (); |
462 | 282 |
283 return info; | |
284 } | |
285 | |
5275 | 286 octave_idx_type |
4725 | 287 EIG::init (const ComplexMatrix& a, bool calc_ev) |
462 | 288 { |
5822 | 289 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
|
290 (*current_liboctave_error_handler) |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
291 ("EIG: matrix contains Inf or NaN values"); |
5822 | 292 |
2815 | 293 if (a.is_hermitian ()) |
4725 | 294 return hermitian_init (a, calc_ev); |
2815 | 295 |
5275 | 296 octave_idx_type n = a.rows (); |
1934 | 297 |
298 if (n != a.cols ()) | |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
299 (*current_liboctave_error_handler) ("EIG requires square matrix"); |
462 | 300 |
5275 | 301 octave_idx_type info = 0; |
462 | 302 |
1934 | 303 ComplexMatrix atmp = a; |
304 Complex *tmp_data = atmp.fortran_vec (); | |
462 | 305 |
2815 | 306 ComplexColumnVector w (n); |
307 Complex *pw = w.fortran_vec (); | |
308 | |
5275 | 309 octave_idx_type nvr = calc_ev ? n : 0; |
4725 | 310 ComplexMatrix vtmp (nvr, nvr); |
2815 | 311 Complex *pv = vtmp.fortran_vec (); |
312 | |
5275 | 313 octave_idx_type lwork = -1; |
4800 | 314 Complex dummy_work; |
1934 | 315 |
5275 | 316 octave_idx_type lrwork = 2*n; |
11570
57632dea2446
attempt better backward compatibility for Array constructors
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
317 Array<double> rwork (dim_vector (lrwork, 1)); |
1934 | 318 double *prwork = rwork.fortran_vec (); |
462 | 319 |
1365 | 320 Complex *dummy = 0; |
5275 | 321 octave_idx_type idummy = 1; |
462 | 322 |
4552 | 323 F77_XFCN (zgeev, ZGEEV, (F77_CONST_CHAR_ARG2 ("N", 1), |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
324 F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
325 n, tmp_data, n, pw, dummy, idummy, |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
326 pv, n, &dummy_work, lwork, prwork, info |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
327 F77_CHAR_ARG_LEN (1) |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
328 F77_CHAR_ARG_LEN (1))); |
2815 | 329 |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
330 if (info != 0) |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
331 (*current_liboctave_error_handler) ("zgeev workspace query failed"); |
4800 | 332 |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
333 lwork = static_cast<octave_idx_type> (dummy_work.real ()); |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
334 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
|
335 Complex *pwork = work.fortran_vec (); |
4800 | 336 |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
337 F77_XFCN (zgeev, ZGEEV, (F77_CONST_CHAR_ARG2 ("N", 1), |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
338 F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
339 n, tmp_data, n, pw, dummy, idummy, |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
340 pv, n, pwork, lwork, prwork, info |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
341 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
|
342 F77_CHAR_ARG_LEN (1))); |
4800 | 343 |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
344 if (info < 0) |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
345 (*current_liboctave_error_handler) ("unrecoverable error in zgeev"); |
4800 | 346 |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
347 if (info > 0) |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
348 (*current_liboctave_error_handler) ("zgeev failed to converge"); |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
349 |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
350 lambda = w; |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
351 v = vtmp; |
2815 | 352 |
353 return info; | |
354 } | |
355 | |
5275 | 356 octave_idx_type |
4725 | 357 EIG::hermitian_init (const ComplexMatrix& a, bool calc_ev) |
2815 | 358 { |
5275 | 359 octave_idx_type n = a.rows (); |
2815 | 360 |
361 if (n != a.cols ()) | |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
362 (*current_liboctave_error_handler) ("EIG requires square matrix"); |
2815 | 363 |
5275 | 364 octave_idx_type info = 0; |
462 | 365 |
2815 | 366 ComplexMatrix atmp = a; |
367 Complex *tmp_data = atmp.fortran_vec (); | |
368 | |
3585 | 369 ColumnVector wr (n); |
370 double *pwr = wr.fortran_vec (); | |
2815 | 371 |
5275 | 372 octave_idx_type lwork = -1; |
4800 | 373 Complex dummy_work; |
2815 | 374 |
5275 | 375 octave_idx_type lrwork = 3*n; |
11570
57632dea2446
attempt better backward compatibility for Array constructors
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
376 Array<double> rwork (dim_vector (lrwork, 1)); |
2815 | 377 double *prwork = rwork.fortran_vec (); |
378 | |
4725 | 379 F77_XFCN (zheev, ZHEEV, (F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
380 F77_CONST_CHAR_ARG2 ("U", 1), |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
381 n, tmp_data, n, pwr, &dummy_work, lwork, |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
382 prwork, info |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
383 F77_CHAR_ARG_LEN (1) |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
384 F77_CHAR_ARG_LEN (1))); |
2815 | 385 |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
386 if (info != 0) |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
387 (*current_liboctave_error_handler) ("zheev workspace query failed"); |
4800 | 388 |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
389 lwork = static_cast<octave_idx_type> (dummy_work.real ()); |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
390 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
|
391 Complex *pwork = work.fortran_vec (); |
4800 | 392 |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
393 F77_XFCN (zheev, ZHEEV, (F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
394 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
|
395 n, tmp_data, n, pwr, pwork, lwork, prwork, info |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
396 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
|
397 F77_CHAR_ARG_LEN (1))); |
4800 | 398 |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
399 if (info < 0) |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
400 (*current_liboctave_error_handler) ("unrecoverable error in zheev"); |
4800 | 401 |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
402 if (info > 0) |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
403 (*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
|
404 |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
405 lambda = ComplexColumnVector (wr); |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
406 v = calc_ev ? ComplexMatrix (atmp) : ComplexMatrix (); |
462 | 407 |
408 return info; | |
409 } | |
410 | |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
411 octave_idx_type |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
412 EIG::init (const Matrix& a, const Matrix& b, bool calc_ev) |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
413 { |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
414 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
|
415 (*current_liboctave_error_handler) |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
416 ("EIG: matrix contains Inf or NaN values"); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
417 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
418 octave_idx_type n = a.rows (); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
419 octave_idx_type nb = b.rows (); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
420 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
421 if (n != a.cols () || nb != b.cols ()) |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
422 (*current_liboctave_error_handler) ("EIG requires square matrix"); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
423 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
424 if (n != nb) |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
425 (*current_liboctave_error_handler) ("EIG requires same size matrices"); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
426 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
427 octave_idx_type info = 0; |
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 tmp = b; |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
430 double *tmp_data = tmp.fortran_vec (); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
431 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
432 F77_XFCN (dpotrf, DPOTRF, (F77_CONST_CHAR_ARG2 ("L", 1), |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11570
diff
changeset
|
433 n, tmp_data, n, |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
434 info |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
435 F77_CHAR_ARG_LEN (1) |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
436 F77_CHAR_ARG_LEN (1))); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
437 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
438 if (a.is_symmetric () && b.is_symmetric () && info == 0) |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
439 return symmetric_init (a, b, calc_ev); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
440 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
441 Matrix atmp = a; |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
442 double *atmp_data = atmp.fortran_vec (); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
443 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
444 Matrix btmp = b; |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
445 double *btmp_data = btmp.fortran_vec (); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
446 |
11570
57632dea2446
attempt better backward compatibility for Array constructors
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
447 Array<double> ar (dim_vector (n, 1)); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
448 double *par = ar.fortran_vec (); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
449 |
11570
57632dea2446
attempt better backward compatibility for Array constructors
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
450 Array<double> ai (dim_vector (n, 1)); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
451 double *pai = ai.fortran_vec (); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
452 |
11570
57632dea2446
attempt better backward compatibility for Array constructors
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
453 Array<double> beta (dim_vector (n, 1)); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
454 double *pbeta = beta.fortran_vec (); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
455 |
9020
728e7943752d
EIG.cc: avoid volatile decl for tmp variable
John W. Eaton <jwe@octave.org>
parents:
8920
diff
changeset
|
456 octave_idx_type tnvr = calc_ev ? n : 0; |
728e7943752d
EIG.cc: avoid volatile decl for tmp variable
John W. Eaton <jwe@octave.org>
parents:
8920
diff
changeset
|
457 Matrix vr (tnvr, tnvr); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
458 double *pvr = vr.fortran_vec (); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
459 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
460 octave_idx_type lwork = -1; |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
461 double dummy_work; |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
462 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
463 double *dummy = 0; |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
464 octave_idx_type idummy = 1; |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
465 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
466 F77_XFCN (dggev, DGGEV, (F77_CONST_CHAR_ARG2 ("N", 1), |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
467 F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11570
diff
changeset
|
468 n, atmp_data, n, btmp_data, n, |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
469 par, pai, pbeta, |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
470 dummy, idummy, pvr, n, |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
471 &dummy_work, lwork, info |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
472 F77_CHAR_ARG_LEN (1) |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
473 F77_CHAR_ARG_LEN (1))); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
474 |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
475 if (info != 0) |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
476 (*current_liboctave_error_handler) ("dggev workspace query failed"); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
477 |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
478 lwork = static_cast<octave_idx_type> (dummy_work); |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
479 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
|
480 double *pwork = work.fortran_vec (); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
481 |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
482 F77_XFCN (dggev, DGGEV, (F77_CONST_CHAR_ARG2 ("N", 1), |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
483 F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
484 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
|
485 par, pai, pbeta, |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
486 dummy, idummy, pvr, n, |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
487 pwork, lwork, info |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
488 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
|
489 F77_CHAR_ARG_LEN (1))); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
490 |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
491 if (info < 0) |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
492 (*current_liboctave_error_handler) ("unrecoverable error in dggev"); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
493 |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
494 if (info > 0) |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
495 (*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
|
496 |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
497 lambda.resize (n); |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
498 octave_idx_type nvr = calc_ev ? n : 0; |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
499 v.resize (nvr, nvr); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
500 |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
501 for (octave_idx_type j = 0; j < n; j++) |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
502 { |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
503 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
|
504 { |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
505 lambda.elem (j) = Complex (ar.elem (j) / beta.elem (j)); |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
506 for (octave_idx_type i = 0; i < nvr; i++) |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
507 v.elem (i, j) = vr.elem (i, j); |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
508 } |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
509 else |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
510 { |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
511 if (j+1 >= n) |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
512 (*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
|
513 |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
514 lambda.elem (j) = Complex (ar.elem (j) / beta.elem (j), |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
515 ai.elem (j) / beta.elem (j)); |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
516 lambda.elem (j+1) = Complex (ar.elem (j+1) / beta.elem (j+1), |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
517 ai.elem (j+1) / beta.elem (j+1)); |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
518 |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
519 for (octave_idx_type i = 0; i < nvr; i++) |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
520 { |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
521 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
|
522 double imag_part = vr.elem (i, j+1); |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
523 v.elem (i, j) = Complex (real_part, imag_part); |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
524 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
|
525 } |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
526 j++; |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
527 } |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
528 } |
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 return info; |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
531 } |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
532 |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11570
diff
changeset
|
533 octave_idx_type |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
534 EIG::symmetric_init (const Matrix& a, const Matrix& b, bool calc_ev) |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
535 { |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
536 octave_idx_type n = a.rows (); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
537 octave_idx_type nb = b.rows (); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
538 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
539 if (n != a.cols () || nb != b.cols ()) |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
540 (*current_liboctave_error_handler) ("EIG requires square matrix"); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
541 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
542 if (n != nb) |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
543 (*current_liboctave_error_handler) ("EIG requires same size matrices"); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
544 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
545 octave_idx_type info = 0; |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
546 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
547 Matrix atmp = a; |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
548 double *atmp_data = atmp.fortran_vec (); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
549 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
550 Matrix btmp = b; |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
551 double *btmp_data = btmp.fortran_vec (); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
552 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
553 ColumnVector wr (n); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
554 double *pwr = wr.fortran_vec (); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
555 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
556 octave_idx_type lwork = -1; |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
557 double dummy_work; |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
558 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
559 F77_XFCN (dsygv, DSYGV, (1, F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
560 F77_CONST_CHAR_ARG2 ("U", 1), |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11570
diff
changeset
|
561 n, atmp_data, n, |
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11570
diff
changeset
|
562 btmp_data, n, |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
563 pwr, &dummy_work, lwork, info |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
564 F77_CHAR_ARG_LEN (1) |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
565 F77_CHAR_ARG_LEN (1))); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
566 |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
567 if (info != 0) |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
568 (*current_liboctave_error_handler) ("dsygv workspace query failed"); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
569 |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
570 lwork = static_cast<octave_idx_type> (dummy_work); |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
571 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
|
572 double *pwork = work.fortran_vec (); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
573 |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
574 F77_XFCN (dsygv, DSYGV, (1, F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
575 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
|
576 n, atmp_data, n, |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
577 btmp_data, n, |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
578 pwr, pwork, lwork, info |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
579 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
|
580 F77_CHAR_ARG_LEN (1))); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
581 |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
582 if (info < 0) |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
583 (*current_liboctave_error_handler) ("unrecoverable error in dsygv"); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
584 |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
585 if (info > 0) |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
586 (*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
|
587 |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
588 lambda = ComplexColumnVector (wr); |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
589 v = calc_ev ? ComplexMatrix (atmp) : ComplexMatrix (); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
590 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
591 return info; |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
592 } |
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 octave_idx_type |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
595 EIG::init (const ComplexMatrix& a, const ComplexMatrix& b, bool calc_ev) |
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 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
|
598 (*current_liboctave_error_handler) |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
599 ("EIG: matrix contains Inf or NaN values"); |
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 octave_idx_type n = a.rows (); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
602 octave_idx_type nb = b.rows (); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
603 |
14846
460a3c6d8bf1
maint: Use Octave coding convention for cuddled parenthis in function calls with empty argument lists.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
604 if (n != a.cols () || nb != b.cols ()) |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
605 (*current_liboctave_error_handler) ("EIG requires square matrix"); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
606 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
607 if (n != nb) |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
608 (*current_liboctave_error_handler) ("EIG requires same size matrices"); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
609 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
610 octave_idx_type info = 0; |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
611 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
612 ComplexMatrix tmp = b; |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
613 Complex*tmp_data = tmp.fortran_vec (); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
614 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
615 F77_XFCN (zpotrf, ZPOTRF, (F77_CONST_CHAR_ARG2 ("L", 1), |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11570
diff
changeset
|
616 n, tmp_data, n, |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
617 info |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
618 F77_CHAR_ARG_LEN (1) |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
619 F77_CHAR_ARG_LEN (1))); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
620 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
621 if (a.is_hermitian () && b.is_hermitian () && info == 0) |
20490
ff904ae0285b
eig: Return correct solution for a pair of hermitian matrices (bug #45511)
Mike Miller <mtmiller@octave.org>
parents:
19697
diff
changeset
|
622 return hermitian_init (a, b, calc_ev); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
623 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
624 ComplexMatrix atmp = a; |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
625 Complex *atmp_data = atmp.fortran_vec (); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
626 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
627 ComplexMatrix btmp = b; |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
628 Complex *btmp_data = btmp.fortran_vec (); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
629 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
630 ComplexColumnVector alpha (n); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
631 Complex *palpha = alpha.fortran_vec (); |
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 ComplexColumnVector beta (n); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
634 Complex *pbeta = beta.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 octave_idx_type nvr = calc_ev ? n : 0; |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
637 ComplexMatrix vtmp (nvr, nvr); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
638 Complex *pv = vtmp.fortran_vec (); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
639 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
640 octave_idx_type lwork = -1; |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
641 Complex dummy_work; |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
642 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
643 octave_idx_type lrwork = 8*n; |
11570
57632dea2446
attempt better backward compatibility for Array constructors
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
644 Array<double> rwork (dim_vector (lrwork, 1)); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
645 double *prwork = rwork.fortran_vec (); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
646 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
647 Complex *dummy = 0; |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
648 octave_idx_type idummy = 1; |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
649 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
650 F77_XFCN (zggev, ZGGEV, (F77_CONST_CHAR_ARG2 ("N", 1), |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
651 F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11570
diff
changeset
|
652 n, atmp_data, n, btmp_data, n, |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
653 palpha, pbeta, dummy, idummy, |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
654 pv, n, &dummy_work, lwork, prwork, info |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
655 F77_CHAR_ARG_LEN (1) |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
656 F77_CHAR_ARG_LEN (1))); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
657 |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
658 if (info != 0) |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
659 (*current_liboctave_error_handler) ("zggev workspace query failed"); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
660 |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
661 lwork = static_cast<octave_idx_type> (dummy_work.real ()); |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
662 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
|
663 Complex *pwork = work.fortran_vec (); |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11570
diff
changeset
|
664 |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
665 F77_XFCN (zggev, ZGGEV, (F77_CONST_CHAR_ARG2 ("N", 1), |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
666 F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
667 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
|
668 palpha, pbeta, dummy, idummy, |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
669 pv, n, pwork, lwork, prwork, info |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
670 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
|
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) ("unrecoverable error in zggev"); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
675 |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
676 if (info > 0) |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
677 (*current_liboctave_error_handler) ("zggev failed to converge"); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
678 |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
679 lambda.resize (n); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
680 |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
681 for (octave_idx_type j = 0; j < n; j++) |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
682 lambda.elem (j) = alpha.elem (j) / beta.elem (j); |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
683 |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
684 v = vtmp; |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
685 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
686 return info; |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
687 } |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
688 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
689 octave_idx_type |
17769
49a5a4be04a1
maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
690 EIG::hermitian_init (const ComplexMatrix& a, const ComplexMatrix& b, |
49a5a4be04a1
maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
691 bool calc_ev) |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
692 { |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
693 octave_idx_type n = a.rows (); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
694 octave_idx_type nb = b.rows (); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
695 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
696 if (n != a.cols () || nb != b.cols ()) |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
697 (*current_liboctave_error_handler) ("EIG requires square matrix"); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
698 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
699 if (n != nb) |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
700 (*current_liboctave_error_handler) ("EIG requires same size matrices"); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
701 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
702 octave_idx_type info = 0; |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
703 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
704 ComplexMatrix atmp = a; |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
705 Complex *atmp_data = atmp.fortran_vec (); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
706 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
707 ComplexMatrix btmp = b; |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
708 Complex *btmp_data = btmp.fortran_vec (); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
709 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
710 ColumnVector wr (n); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
711 double *pwr = wr.fortran_vec (); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
712 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
713 octave_idx_type lwork = -1; |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
714 Complex dummy_work; |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
715 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
716 octave_idx_type lrwork = 3*n; |
11570
57632dea2446
attempt better backward compatibility for Array constructors
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
717 Array<double> rwork (dim_vector (lrwork, 1)); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
718 double *prwork = rwork.fortran_vec (); |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
719 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
720 F77_XFCN (zhegv, ZHEGV, (1, F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
721 F77_CONST_CHAR_ARG2 ("U", 1), |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11570
diff
changeset
|
722 n, atmp_data, n, |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
723 btmp_data, n, |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
724 pwr, &dummy_work, lwork, |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
725 prwork, info |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
726 F77_CHAR_ARG_LEN (1) |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
727 F77_CHAR_ARG_LEN (1))); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
728 |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
729 if (info != 0) |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
730 (*current_liboctave_error_handler) ("zhegv workspace query failed"); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
731 |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
732 lwork = static_cast<octave_idx_type> (dummy_work.real ()); |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
733 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
|
734 Complex *pwork = work.fortran_vec (); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
735 |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
736 F77_XFCN (zhegv, ZHEGV, (1, F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1), |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
737 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
|
738 n, atmp_data, n, |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
739 btmp_data, n, |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
740 pwr, pwork, lwork, prwork, info |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
741 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
|
742 F77_CHAR_ARG_LEN (1))); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
743 |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
744 if (info < 0) |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
745 (*current_liboctave_error_handler) ("unrecoverable error in zhegv"); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
746 |
21136
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
747 if (info > 0) |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
748 (*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
|
749 |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
750 lambda = ComplexColumnVector (wr); |
7cac4e7458f2
maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents:
20490
diff
changeset
|
751 v = calc_ev ? ComplexMatrix (atmp) : ComplexMatrix (); |
8339
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
752 |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
753 return info; |
18c4ded8612a
Add generalized eigenvalue functions
Jarkko Kaleva <d3roga@gmail.com>
parents:
7482
diff
changeset
|
754 } |