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