File: | liboctave/numeric/fEIG.cc |
Location: | line 279, 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 "fEIG.h" | |||
28 | #include "fColVector.h" | |||
29 | #include "f77-fcn.h" | |||
30 | #include "lo-error.h" | |||
31 | ||||
32 | extern "C" | |||
33 | { | |||
34 | F77_RET_Tint | |||
35 | F77_FUNC (sgeev, SGEEV)sgeev_ (F77_CONST_CHAR_ARG_DECLconst char *, | |||
36 | F77_CONST_CHAR_ARG_DECLconst char *, | |||
37 | const octave_idx_type&, float*, | |||
38 | const octave_idx_type&, float*, float*, float*, | |||
39 | const octave_idx_type&, float*, | |||
40 | const octave_idx_type&, float*, | |||
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 (cgeev, CGEEV)cgeev_ (F77_CONST_CHAR_ARG_DECLconst char *, | |||
47 | F77_CONST_CHAR_ARG_DECLconst char *, | |||
48 | const octave_idx_type&, FloatComplex*, | |||
49 | const octave_idx_type&, FloatComplex*, FloatComplex*, | |||
50 | const octave_idx_type&, FloatComplex*, | |||
51 | const octave_idx_type&, FloatComplex*, | |||
52 | const octave_idx_type&, float*, 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 (ssyev, SSYEV)ssyev_ (F77_CONST_CHAR_ARG_DECLconst char *, | |||
58 | F77_CONST_CHAR_ARG_DECLconst char *, | |||
59 | const octave_idx_type&, float*, | |||
60 | const octave_idx_type&, float*, float*, | |||
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 (cheev, CHEEV)cheev_ (F77_CONST_CHAR_ARG_DECLconst char *, | |||
67 | F77_CONST_CHAR_ARG_DECLconst char *, | |||
68 | const octave_idx_type&, FloatComplex*, | |||
69 | const octave_idx_type&, float*, FloatComplex*, | |||
70 | const octave_idx_type&, float*, octave_idx_type& | |||
71 | F77_CHAR_ARG_LEN_DECL, long | |||
72 | F77_CHAR_ARG_LEN_DECL, long); | |||
73 | ||||
74 | F77_RET_Tint | |||
75 | F77_FUNC (spotrf, SPOTRF)spotrf_ (F77_CONST_CHAR_ARG_DECLconst char *, | |||
76 | const octave_idx_type&, float*, | |||
77 | const octave_idx_type&, octave_idx_type& | |||
78 | F77_CHAR_ARG_LEN_DECL, long | |||
79 | F77_CHAR_ARG_LEN_DECL, long); | |||
80 | ||||
81 | F77_RET_Tint | |||
82 | F77_FUNC (cpotrf, CPOTRF)cpotrf_ (F77_CONST_CHAR_ARG_DECLconst char *, | |||
83 | const octave_idx_type&, FloatComplex*, | |||
84 | const octave_idx_type&, octave_idx_type& | |||
85 | F77_CHAR_ARG_LEN_DECL, long | |||
86 | F77_CHAR_ARG_LEN_DECL, long); | |||
87 | ||||
88 | F77_RET_Tint | |||
89 | F77_FUNC (sggev, SGGEV)sggev_ (F77_CONST_CHAR_ARG_DECLconst char *, | |||
90 | F77_CONST_CHAR_ARG_DECLconst char *, | |||
91 | const octave_idx_type&, float*, | |||
92 | const octave_idx_type&, float*, | |||
93 | const octave_idx_type&, float*, float*, float*, | |||
94 | float*, const octave_idx_type&, float*, | |||
95 | const octave_idx_type&, float*, | |||
96 | const octave_idx_type&, octave_idx_type& | |||
97 | F77_CHAR_ARG_LEN_DECL, long | |||
98 | F77_CHAR_ARG_LEN_DECL, long); | |||
99 | ||||
100 | F77_RET_Tint | |||
101 | F77_FUNC (ssygv, SSYGV)ssygv_ (const octave_idx_type&, | |||
102 | F77_CONST_CHAR_ARG_DECLconst char *, | |||
103 | F77_CONST_CHAR_ARG_DECLconst char *, | |||
104 | const octave_idx_type&, float*, | |||
105 | const octave_idx_type&, float*, | |||
106 | const octave_idx_type&, float*, float*, | |||
107 | const octave_idx_type&, octave_idx_type& | |||
108 | F77_CHAR_ARG_LEN_DECL, long | |||
109 | F77_CHAR_ARG_LEN_DECL, long); | |||
110 | ||||
111 | F77_RET_Tint | |||
112 | F77_FUNC (cggev, CGGEV)cggev_ (F77_CONST_CHAR_ARG_DECLconst char *, | |||
113 | F77_CONST_CHAR_ARG_DECLconst char *, | |||
114 | const octave_idx_type&, FloatComplex*, | |||
115 | const octave_idx_type&, FloatComplex*, | |||
116 | const octave_idx_type&, FloatComplex*, | |||
117 | FloatComplex*, FloatComplex*, | |||
118 | const octave_idx_type&, FloatComplex*, | |||
119 | const octave_idx_type&, FloatComplex*, | |||
120 | const octave_idx_type&, float*, octave_idx_type& | |||
121 | F77_CHAR_ARG_LEN_DECL, long | |||
122 | F77_CHAR_ARG_LEN_DECL, long); | |||
123 | ||||
124 | F77_RET_Tint | |||
125 | F77_FUNC (chegv, CHEGV)chegv_ (const octave_idx_type&, | |||
126 | F77_CONST_CHAR_ARG_DECLconst char *, | |||
127 | F77_CONST_CHAR_ARG_DECLconst char *, | |||
128 | const octave_idx_type&, FloatComplex*, | |||
129 | const octave_idx_type&, FloatComplex*, | |||
130 | const octave_idx_type&, float*, FloatComplex*, | |||
131 | const octave_idx_type&, float*, octave_idx_type& | |||
132 | F77_CHAR_ARG_LEN_DECL, long | |||
133 | F77_CHAR_ARG_LEN_DECL, long); | |||
134 | } | |||
135 | ||||
136 | octave_idx_type | |||
137 | FloatEIG::init (const FloatMatrix& a, bool calc_ev) | |||
138 | { | |||
139 | if (a.any_element_is_inf_or_nan ()) | |||
140 | { | |||
141 | (*current_liboctave_error_handler) | |||
142 | ("EIG: matrix contains Inf or NaN values"); | |||
143 | return -1; | |||
144 | } | |||
145 | ||||
146 | if (a.is_symmetric ()) | |||
147 | return symmetric_init (a, calc_ev); | |||
148 | ||||
149 | octave_idx_type n = a.rows (); | |||
150 | ||||
151 | if (n != a.cols ()) | |||
152 | { | |||
153 | (*current_liboctave_error_handler) ("EIG requires square matrix"); | |||
154 | return -1; | |||
155 | } | |||
156 | ||||
157 | octave_idx_type info = 0; | |||
158 | ||||
159 | FloatMatrix atmp = a; | |||
160 | float *tmp_data = atmp.fortran_vec (); | |||
161 | ||||
162 | Array<float> wr (dim_vector (n, 1)); | |||
163 | float *pwr = wr.fortran_vec (); | |||
164 | ||||
165 | Array<float> wi (dim_vector (n, 1)); | |||
166 | float *pwi = wi.fortran_vec (); | |||
167 | ||||
168 | volatile octave_idx_type nvr = calc_ev ? n : 0; | |||
169 | FloatMatrix vr (nvr, nvr); | |||
170 | float *pvr = vr.fortran_vec (); | |||
171 | ||||
172 | octave_idx_type lwork = -1; | |||
173 | float dummy_work; | |||
174 | ||||
175 | float *dummy = 0; | |||
176 | octave_idx_type idummy = 1; | |||
177 | ||||
178 | F77_XFCN (sgeev, SGEEV, (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" , "sgeev_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; sgeev_ ("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) | |||
179 | 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" , "sgeev_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; sgeev_ ("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) | |||
180 | 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" , "sgeev_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; sgeev_ ("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) | |||
181 | 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" , "sgeev_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; sgeev_ ("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_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" , "sgeev_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; sgeev_ ("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 | 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" , "sgeev_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; sgeev_ ("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 | ||||
185 | if (info == 0) | |||
186 | { | |||
187 | lwork = static_cast<octave_idx_type> (dummy_work); | |||
188 | Array<float> work (dim_vector (lwork, 1)); | |||
189 | float *pwork = work.fortran_vec (); | |||
190 | ||||
191 | F77_XFCN (sgeev, SGEEV, (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" , "sgeev_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; sgeev_ ("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) | |||
192 | 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" , "sgeev_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; sgeev_ ("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) | |||
193 | 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" , "sgeev_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; sgeev_ ("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) | |||
194 | 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" , "sgeev_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; sgeev_ ("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_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" , "sgeev_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; sgeev_ ("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 | 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" , "sgeev_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; sgeev_ ("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 | ||||
198 | if (info < 0) | |||
199 | { | |||
200 | (*current_liboctave_error_handler) ("unrecoverable error in sgeev"); | |||
201 | return info; | |||
202 | } | |||
203 | ||||
204 | if (info > 0) | |||
205 | { | |||
206 | (*current_liboctave_error_handler) ("sgeev failed to converge"); | |||
207 | return info; | |||
208 | } | |||
209 | ||||
210 | lambda.resize (n); | |||
211 | v.resize (nvr, nvr); | |||
212 | ||||
213 | for (octave_idx_type j = 0; j < n; j++) | |||
214 | { | |||
215 | if (wi.elem (j) == 0.0) | |||
216 | { | |||
217 | lambda.elem (j) = FloatComplex (wr.elem (j)); | |||
218 | for (octave_idx_type i = 0; i < nvr; i++) | |||
219 | v.elem (i, j) = vr.elem (i, j); | |||
220 | } | |||
221 | else | |||
222 | { | |||
223 | if (j+1 >= n) | |||
224 | { | |||
225 | (*current_liboctave_error_handler) ("EIG: internal error"); | |||
226 | return -1; | |||
227 | } | |||
228 | ||||
229 | lambda.elem (j) = FloatComplex (wr.elem (j), wi.elem (j)); | |||
230 | lambda.elem (j+1) = FloatComplex (wr.elem (j+1), wi.elem (j+1)); | |||
231 | ||||
232 | for (octave_idx_type i = 0; i < nvr; i++) | |||
233 | { | |||
234 | float real_part = vr.elem (i, j); | |||
235 | float imag_part = vr.elem (i, j+1); | |||
236 | v.elem (i, j) = FloatComplex (real_part, imag_part); | |||
237 | v.elem (i, j+1) = FloatComplex (real_part, -imag_part); | |||
238 | } | |||
239 | j++; | |||
240 | } | |||
241 | } | |||
242 | } | |||
243 | else | |||
244 | (*current_liboctave_error_handler) ("sgeev workspace query failed"); | |||
245 | ||||
246 | return info; | |||
247 | } | |||
248 | ||||
249 | octave_idx_type | |||
250 | FloatEIG::symmetric_init (const FloatMatrix& a, bool calc_ev) | |||
251 | { | |||
252 | octave_idx_type n = a.rows (); | |||
253 | ||||
254 | if (n != a.cols ()) | |||
| ||||
255 | { | |||
256 | (*current_liboctave_error_handler) ("EIG requires square matrix"); | |||
257 | return -1; | |||
258 | } | |||
259 | ||||
260 | octave_idx_type info = 0; | |||
261 | ||||
262 | FloatMatrix atmp = a; | |||
263 | float *tmp_data = atmp.fortran_vec (); | |||
264 | ||||
265 | FloatColumnVector wr (n); | |||
266 | float *pwr = wr.fortran_vec (); | |||
267 | ||||
268 | octave_idx_type lwork = -1; | |||
269 | float dummy_work; | |||
270 | ||||
271 | F77_XFCN (ssyev, SSYEV, (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" , "ssyev_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; ssyev_ (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) | |||
272 | 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" , "ssyev_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; ssyev_ (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) | |||
273 | 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" , "ssyev_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; ssyev_ (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) | |||
274 | 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" , "ssyev_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; ssyev_ (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) | |||
275 | 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" , "ssyev_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; ssyev_ (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 | ||||
277 | if (info == 0) | |||
278 | { | |||
279 | lwork = static_cast<octave_idx_type> (dummy_work); | |||
| ||||
280 | Array<float> work (dim_vector (lwork, 1)); | |||
281 | float *pwork = work.fortran_vec (); | |||
282 | ||||
283 | F77_XFCN (ssyev, SSYEV, (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" , "ssyev_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; ssyev_ (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) | |||
284 | 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" , "ssyev_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; ssyev_ (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) | |||
285 | 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" , "ssyev_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; ssyev_ (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) | |||
286 | 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" , "ssyev_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; ssyev_ (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) | |||
287 | 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" , "ssyev_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; ssyev_ (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 | ||||
289 | if (info < 0) | |||
290 | { | |||
291 | (*current_liboctave_error_handler) ("unrecoverable error in ssyev"); | |||
292 | return info; | |||
293 | } | |||
294 | ||||
295 | if (info > 0) | |||
296 | { | |||
297 | (*current_liboctave_error_handler) ("ssyev failed to converge"); | |||
298 | return info; | |||
299 | } | |||
300 | ||||
301 | lambda = FloatComplexColumnVector (wr); | |||
302 | v = calc_ev ? FloatComplexMatrix (atmp) : FloatComplexMatrix (); | |||
303 | } | |||
304 | else | |||
305 | (*current_liboctave_error_handler) ("ssyev workspace query failed"); | |||
306 | ||||
307 | return info; | |||
308 | } | |||
309 | ||||
310 | octave_idx_type | |||
311 | FloatEIG::init (const FloatComplexMatrix& a, bool calc_ev) | |||
312 | { | |||
313 | if (a.any_element_is_inf_or_nan ()) | |||
314 | { | |||
315 | (*current_liboctave_error_handler) | |||
316 | ("EIG: matrix contains Inf or NaN values"); | |||
317 | return -1; | |||
318 | } | |||
319 | ||||
320 | if (a.is_hermitian ()) | |||
321 | return hermitian_init (a, calc_ev); | |||
322 | ||||
323 | octave_idx_type n = a.rows (); | |||
324 | ||||
325 | if (n != a.cols ()) | |||
326 | { | |||
327 | (*current_liboctave_error_handler) ("EIG requires square matrix"); | |||
328 | return -1; | |||
329 | } | |||
330 | ||||
331 | octave_idx_type info = 0; | |||
332 | ||||
333 | FloatComplexMatrix atmp = a; | |||
334 | FloatComplex *tmp_data = atmp.fortran_vec (); | |||
335 | ||||
336 | FloatComplexColumnVector w (n); | |||
337 | FloatComplex *pw = w.fortran_vec (); | |||
338 | ||||
339 | octave_idx_type nvr = calc_ev ? n : 0; | |||
340 | FloatComplexMatrix vtmp (nvr, nvr); | |||
341 | FloatComplex *pv = vtmp.fortran_vec (); | |||
342 | ||||
343 | octave_idx_type lwork = -1; | |||
344 | FloatComplex dummy_work; | |||
345 | ||||
346 | octave_idx_type lrwork = 2*n; | |||
347 | Array<float> rwork (dim_vector (lrwork, 1)); | |||
348 | float *prwork = rwork.fortran_vec (); | |||
349 | ||||
350 | FloatComplex *dummy = 0; | |||
351 | octave_idx_type idummy = 1; | |||
352 | ||||
353 | F77_XFCN (cgeev, CGEEV, (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" , "cgeev_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; cgeev_ ("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) | |||
354 | 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" , "cgeev_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; cgeev_ ("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) | |||
355 | 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" , "cgeev_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; cgeev_ ("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) | |||
356 | 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" , "cgeev_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; cgeev_ ("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) | |||
357 | 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" , "cgeev_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; cgeev_ ("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_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" , "cgeev_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; cgeev_ ("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 | ||||
360 | if (info == 0) | |||
361 | { | |||
362 | lwork = static_cast<octave_idx_type> (dummy_work.real ()); | |||
363 | Array<FloatComplex> work (dim_vector (lwork, 1)); | |||
364 | FloatComplex *pwork = work.fortran_vec (); | |||
365 | ||||
366 | F77_XFCN (cgeev, CGEEV, (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" , "cgeev_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; cgeev_ ("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) | |||
367 | 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" , "cgeev_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; cgeev_ ("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) | |||
368 | 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" , "cgeev_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; cgeev_ ("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) | |||
369 | 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" , "cgeev_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; cgeev_ ("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) | |||
370 | 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" , "cgeev_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; cgeev_ ("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_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" , "cgeev_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; cgeev_ ("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 | ||||
373 | if (info < 0) | |||
374 | { | |||
375 | (*current_liboctave_error_handler) ("unrecoverable error in cgeev"); | |||
376 | return info; | |||
377 | } | |||
378 | ||||
379 | if (info > 0) | |||
380 | { | |||
381 | (*current_liboctave_error_handler) ("cgeev failed to converge"); | |||
382 | return info; | |||
383 | } | |||
384 | ||||
385 | lambda = w; | |||
386 | v = vtmp; | |||
387 | } | |||
388 | else | |||
389 | (*current_liboctave_error_handler) ("cgeev workspace query failed"); | |||
390 | ||||
391 | return info; | |||
392 | } | |||
393 | ||||
394 | octave_idx_type | |||
395 | FloatEIG::hermitian_init (const FloatComplexMatrix& a, bool calc_ev) | |||
396 | { | |||
397 | octave_idx_type n = a.rows (); | |||
398 | ||||
399 | if (n != a.cols ()) | |||
400 | { | |||
401 | (*current_liboctave_error_handler) ("EIG requires square matrix"); | |||
402 | return -1; | |||
403 | } | |||
404 | ||||
405 | octave_idx_type info = 0; | |||
406 | ||||
407 | FloatComplexMatrix atmp = a; | |||
408 | FloatComplex *tmp_data = atmp.fortran_vec (); | |||
409 | ||||
410 | FloatColumnVector wr (n); | |||
411 | float *pwr = wr.fortran_vec (); | |||
412 | ||||
413 | octave_idx_type lwork = -1; | |||
414 | FloatComplex dummy_work; | |||
415 | ||||
416 | octave_idx_type lrwork = 3*n; | |||
417 | Array<float> rwork (dim_vector (lrwork, 1)); | |||
418 | float *prwork = rwork.fortran_vec (); | |||
419 | ||||
420 | F77_XFCN (cheev, CHEEV, (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" , "cheev_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; cheev_ (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) | |||
421 | 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" , "cheev_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; cheev_ (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) | |||
422 | 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" , "cheev_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; cheev_ (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) | |||
423 | 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" , "cheev_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; cheev_ (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) | |||
424 | 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" , "cheev_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; cheev_ (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_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" , "cheev_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; cheev_ (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 | ||||
427 | if (info == 0) | |||
428 | { | |||
429 | lwork = static_cast<octave_idx_type> (dummy_work.real ()); | |||
430 | Array<FloatComplex> work (dim_vector (lwork, 1)); | |||
431 | FloatComplex *pwork = work.fortran_vec (); | |||
432 | ||||
433 | F77_XFCN (cheev, CHEEV, (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" , "cheev_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; cheev_ (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) | |||
434 | 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" , "cheev_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; cheev_ (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) | |||
435 | 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" , "cheev_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; cheev_ (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) | |||
436 | 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" , "cheev_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; cheev_ (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) | |||
437 | 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" , "cheev_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; cheev_ (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 | ||||
439 | if (info < 0) | |||
440 | { | |||
441 | (*current_liboctave_error_handler) ("unrecoverable error in cheev"); | |||
442 | return info; | |||
443 | } | |||
444 | ||||
445 | if (info > 0) | |||
446 | { | |||
447 | (*current_liboctave_error_handler) ("cheev failed to converge"); | |||
448 | return info; | |||
449 | } | |||
450 | ||||
451 | lambda = FloatComplexColumnVector (wr); | |||
452 | v = calc_ev ? FloatComplexMatrix (atmp) : FloatComplexMatrix (); | |||
453 | } | |||
454 | else | |||
455 | (*current_liboctave_error_handler) ("cheev workspace query failed"); | |||
456 | ||||
457 | return info; | |||
458 | } | |||
459 | ||||
460 | octave_idx_type | |||
461 | FloatEIG::init (const FloatMatrix& a, const FloatMatrix& b, bool calc_ev) | |||
462 | { | |||
463 | if (a.any_element_is_inf_or_nan () || b.any_element_is_inf_or_nan ()) | |||
464 | { | |||
465 | (*current_liboctave_error_handler) | |||
466 | ("EIG: matrix contains Inf or NaN values"); | |||
467 | return -1; | |||
468 | } | |||
469 | ||||
470 | octave_idx_type n = a.rows (); | |||
471 | octave_idx_type nb = b.rows (); | |||
472 | ||||
473 | if (n != a.cols () || nb != b.cols ()) | |||
474 | { | |||
475 | (*current_liboctave_error_handler) ("EIG requires square matrix"); | |||
476 | return -1; | |||
477 | } | |||
478 | ||||
479 | if (n != nb) | |||
480 | { | |||
481 | (*current_liboctave_error_handler) ("EIG requires same size matrices"); | |||
482 | return -1; | |||
483 | } | |||
484 | ||||
485 | octave_idx_type info = 0; | |||
486 | ||||
487 | FloatMatrix tmp = b; | |||
488 | float *tmp_data = tmp.fortran_vec (); | |||
489 | ||||
490 | F77_XFCN (spotrf, SPOTRF, (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" , "spotrf_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; spotrf_ ("L", n, tmp_data, n, info , 1 , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
491 | 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" , "spotrf_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; spotrf_ ("L", n, tmp_data, n, info , 1 , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
492 | 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" , "spotrf_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; spotrf_ ("L", n, tmp_data, n, info , 1 , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
493 | 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" , "spotrf_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; spotrf_ ("L", n, tmp_data, n, info , 1 , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
494 | 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" , "spotrf_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; spotrf_ ("L", n, tmp_data, n, info , 1 , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0); | |||
495 | ||||
496 | if (a.is_symmetric () && b.is_symmetric () && info == 0) | |||
497 | return symmetric_init (a, b, calc_ev); | |||
498 | ||||
499 | FloatMatrix atmp = a; | |||
500 | float *atmp_data = atmp.fortran_vec (); | |||
501 | ||||
502 | FloatMatrix btmp = b; | |||
503 | float *btmp_data = btmp.fortran_vec (); | |||
504 | ||||
505 | Array<float> ar (dim_vector (n, 1)); | |||
506 | float *par = ar.fortran_vec (); | |||
507 | ||||
508 | Array<float> ai (dim_vector (n, 1)); | |||
509 | float *pai = ai.fortran_vec (); | |||
510 | ||||
511 | Array<float> beta (dim_vector (n, 1)); | |||
512 | float *pbeta = beta.fortran_vec (); | |||
513 | ||||
514 | volatile octave_idx_type nvr = calc_ev ? n : 0; | |||
515 | FloatMatrix vr (nvr, nvr); | |||
516 | float *pvr = vr.fortran_vec (); | |||
517 | ||||
518 | octave_idx_type lwork = -1; | |||
519 | float dummy_work; | |||
520 | ||||
521 | float *dummy = 0; | |||
522 | octave_idx_type idummy = 1; | |||
523 | ||||
524 | F77_XFCN (sggev, SGGEV, (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" , "sggev_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; sggev_ ("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) | |||
525 | 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" , "sggev_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; sggev_ ("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) | |||
526 | 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" , "sggev_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; sggev_ ("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) | |||
527 | 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" , "sggev_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; sggev_ ("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) | |||
528 | 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" , "sggev_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; sggev_ ("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 | &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" , "sggev_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; sggev_ ("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 | 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" , "sggev_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; sggev_ ("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 | 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" , "sggev_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; sggev_ ("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 | ||||
533 | if (info == 0) | |||
534 | { | |||
535 | lwork = static_cast<octave_idx_type> (dummy_work); | |||
536 | Array<float> work (dim_vector (lwork, 1)); | |||
537 | float *pwork = work.fortran_vec (); | |||
538 | ||||
539 | F77_XFCN (sggev, SGGEV, (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" , "sggev_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; sggev_ ("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) | |||
540 | 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" , "sggev_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; sggev_ ("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) | |||
541 | 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" , "sggev_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; sggev_ ("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) | |||
542 | 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" , "sggev_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; sggev_ ("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) | |||
543 | 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" , "sggev_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; sggev_ ("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 | 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" , "sggev_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; sggev_ ("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 | 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" , "sggev_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; sggev_ ("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 | 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" , "sggev_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; sggev_ ("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 | ||||
548 | if (info < 0) | |||
549 | { | |||
550 | (*current_liboctave_error_handler) ("unrecoverable error in sggev"); | |||
551 | return info; | |||
552 | } | |||
553 | ||||
554 | if (info > 0) | |||
555 | { | |||
556 | (*current_liboctave_error_handler) ("sggev failed to converge"); | |||
557 | return info; | |||
558 | } | |||
559 | ||||
560 | lambda.resize (n); | |||
561 | v.resize (nvr, nvr); | |||
562 | ||||
563 | for (octave_idx_type j = 0; j < n; j++) | |||
564 | { | |||
565 | if (ai.elem (j) == 0.0) | |||
566 | { | |||
567 | lambda.elem (j) = FloatComplex (ar.elem (j) / beta.elem (j)); | |||
568 | for (octave_idx_type i = 0; i < nvr; i++) | |||
569 | v.elem (i, j) = vr.elem (i, j); | |||
570 | } | |||
571 | else | |||
572 | { | |||
573 | if (j+1 >= n) | |||
574 | { | |||
575 | (*current_liboctave_error_handler) ("EIG: internal error"); | |||
576 | return -1; | |||
577 | } | |||
578 | ||||
579 | lambda.elem (j) = FloatComplex (ar.elem (j) / beta.elem (j), | |||
580 | ai.elem (j) / beta.elem (j)); | |||
581 | lambda.elem (j+1) = FloatComplex (ar.elem (j+1) / beta.elem (j+1), | |||
582 | ai.elem (j+1) / beta.elem (j+1)); | |||
583 | ||||
584 | for (octave_idx_type i = 0; i < nvr; i++) | |||
585 | { | |||
586 | float real_part = vr.elem (i, j); | |||
587 | float imag_part = vr.elem (i, j+1); | |||
588 | v.elem (i, j) = FloatComplex (real_part, imag_part); | |||
589 | v.elem (i, j+1) = FloatComplex (real_part, -imag_part); | |||
590 | } | |||
591 | j++; | |||
592 | } | |||
593 | } | |||
594 | } | |||
595 | else | |||
596 | (*current_liboctave_error_handler) ("sggev workspace query failed"); | |||
597 | ||||
598 | return info; | |||
599 | } | |||
600 | ||||
601 | octave_idx_type | |||
602 | FloatEIG::symmetric_init (const FloatMatrix& a, const FloatMatrix& b, | |||
603 | bool calc_ev) | |||
604 | { | |||
605 | octave_idx_type n = a.rows (); | |||
606 | octave_idx_type nb = b.rows (); | |||
607 | ||||
608 | if (n != a.cols () || nb != b.cols ()) | |||
609 | { | |||
610 | (*current_liboctave_error_handler) ("EIG requires square matrix"); | |||
611 | return -1; | |||
612 | } | |||
613 | ||||
614 | if (n != nb) | |||
615 | { | |||
616 | (*current_liboctave_error_handler) ("EIG requires same size matrices"); | |||
617 | return -1; | |||
618 | } | |||
619 | ||||
620 | octave_idx_type info = 0; | |||
621 | ||||
622 | FloatMatrix atmp = a; | |||
623 | float *atmp_data = atmp.fortran_vec (); | |||
624 | ||||
625 | FloatMatrix btmp = b; | |||
626 | float *btmp_data = btmp.fortran_vec (); | |||
627 | ||||
628 | FloatColumnVector wr (n); | |||
629 | float *pwr = wr.fortran_vec (); | |||
630 | ||||
631 | octave_idx_type lwork = -1; | |||
632 | float dummy_work; | |||
633 | ||||
634 | F77_XFCN (ssygv, SSYGV, (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" , "ssygv_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; ssygv_ (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) | |||
635 | 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" , "ssygv_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; ssygv_ (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) | |||
636 | 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" , "ssygv_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; ssygv_ (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) | |||
637 | 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" , "ssygv_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; ssygv_ (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) | |||
638 | 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" , "ssygv_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; ssygv_ (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_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" , "ssygv_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; ssygv_ (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 | 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" , "ssygv_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; ssygv_ (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 | ||||
642 | if (info == 0) | |||
643 | { | |||
644 | lwork = static_cast<octave_idx_type> (dummy_work); | |||
645 | Array<float> work (dim_vector (lwork, 1)); | |||
646 | float *pwork = work.fortran_vec (); | |||
647 | ||||
648 | F77_XFCN (ssygv, SSYGV, (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" , "ssygv_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; ssygv_ (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) | |||
649 | 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" , "ssygv_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; ssygv_ (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) | |||
650 | 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" , "ssygv_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; ssygv_ (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) | |||
651 | 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" , "ssygv_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; ssygv_ (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) | |||
652 | 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" , "ssygv_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; ssygv_ (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_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" , "ssygv_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; ssygv_ (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 | 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" , "ssygv_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; ssygv_ (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 | ||||
656 | if (info < 0) | |||
657 | { | |||
658 | (*current_liboctave_error_handler) ("unrecoverable error in ssygv"); | |||
659 | return info; | |||
660 | } | |||
661 | ||||
662 | if (info > 0) | |||
663 | { | |||
664 | (*current_liboctave_error_handler) ("ssygv failed to converge"); | |||
665 | return info; | |||
666 | } | |||
667 | ||||
668 | lambda = FloatComplexColumnVector (wr); | |||
669 | v = calc_ev ? FloatComplexMatrix (atmp) : FloatComplexMatrix (); | |||
670 | } | |||
671 | else | |||
672 | (*current_liboctave_error_handler) ("ssygv workspace query failed"); | |||
673 | ||||
674 | return info; | |||
675 | } | |||
676 | ||||
677 | octave_idx_type | |||
678 | FloatEIG::init (const FloatComplexMatrix& a, const FloatComplexMatrix& b, | |||
679 | bool calc_ev) | |||
680 | { | |||
681 | if (a.any_element_is_inf_or_nan () || b.any_element_is_inf_or_nan ()) | |||
682 | { | |||
683 | (*current_liboctave_error_handler) | |||
684 | ("EIG: matrix contains Inf or NaN values"); | |||
685 | return -1; | |||
686 | } | |||
687 | ||||
688 | octave_idx_type n = a.rows (); | |||
689 | octave_idx_type nb = b.rows (); | |||
690 | ||||
691 | if (n != a.cols () || nb != b.cols ()) | |||
692 | { | |||
693 | (*current_liboctave_error_handler) ("EIG requires square matrix"); | |||
694 | return -1; | |||
695 | } | |||
696 | ||||
697 | if (n != nb) | |||
698 | { | |||
699 | (*current_liboctave_error_handler) ("EIG requires same size matrices"); | |||
700 | return -1; | |||
701 | } | |||
702 | ||||
703 | octave_idx_type info = 0; | |||
704 | ||||
705 | FloatComplexMatrix tmp = b; | |||
706 | FloatComplex *tmp_data = tmp.fortran_vec (); | |||
707 | ||||
708 | F77_XFCN (cpotrf, CPOTRF, (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" , "cpotrf_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; cpotrf_ ("L", n, tmp_data, n, info , 1 , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
709 | 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" , "cpotrf_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; cpotrf_ ("L", n, tmp_data, n, info , 1 , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
710 | 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" , "cpotrf_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; cpotrf_ ("L", n, tmp_data, n, info , 1 , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
711 | 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" , "cpotrf_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; cpotrf_ ("L", n, tmp_data, n, info , 1 , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0) | |||
712 | 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" , "cpotrf_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; cpotrf_ ("L", n, tmp_data, n, info , 1 , 1); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0); | |||
713 | ||||
714 | if (a.is_hermitian () && b.is_hermitian () && info == 0) | |||
715 | return hermitian_init (a, calc_ev); | |||
716 | ||||
717 | FloatComplexMatrix atmp = a; | |||
718 | FloatComplex *atmp_data = atmp.fortran_vec (); | |||
719 | ||||
720 | FloatComplexMatrix btmp = b; | |||
721 | FloatComplex *btmp_data = btmp.fortran_vec (); | |||
722 | ||||
723 | FloatComplexColumnVector alpha (n); | |||
724 | FloatComplex *palpha = alpha.fortran_vec (); | |||
725 | ||||
726 | FloatComplexColumnVector beta (n); | |||
727 | FloatComplex *pbeta = beta.fortran_vec (); | |||
728 | ||||
729 | octave_idx_type nvr = calc_ev ? n : 0; | |||
730 | FloatComplexMatrix vtmp (nvr, nvr); | |||
731 | FloatComplex *pv = vtmp.fortran_vec (); | |||
732 | ||||
733 | octave_idx_type lwork = -1; | |||
734 | FloatComplex dummy_work; | |||
735 | ||||
736 | octave_idx_type lrwork = 8*n; | |||
737 | Array<float> rwork (dim_vector (lrwork, 1)); | |||
738 | float *prwork = rwork.fortran_vec (); | |||
739 | ||||
740 | FloatComplex *dummy = 0; | |||
741 | octave_idx_type idummy = 1; | |||
742 | ||||
743 | F77_XFCN (cggev, CGGEV, (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" , "cggev_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; cggev_ ("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) | |||
744 | 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" , "cggev_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; cggev_ ("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) | |||
745 | 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" , "cggev_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; cggev_ ("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) | |||
746 | 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" , "cggev_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; cggev_ ("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 | 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" , "cggev_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; cggev_ ("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 | 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" , "cggev_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; cggev_ ("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 | 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" , "cggev_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; cggev_ ("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 | ||||
751 | if (info == 0) | |||
752 | { | |||
753 | lwork = static_cast<octave_idx_type> (dummy_work.real ()); | |||
754 | Array<FloatComplex> work (dim_vector (lwork, 1)); | |||
755 | FloatComplex *pwork = work.fortran_vec (); | |||
756 | ||||
757 | F77_XFCN (cggev, CGGEV, (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" , "cggev_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; cggev_ ("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) | |||
758 | 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" , "cggev_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; cggev_ ("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) | |||
759 | 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" , "cggev_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; cggev_ ("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) | |||
760 | 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" , "cggev_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; cggev_ ("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 | 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" , "cggev_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; cggev_ ("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 | 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" , "cggev_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; cggev_ ("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 | 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" , "cggev_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; cggev_ ("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 | ||||
765 | if (info < 0) | |||
766 | { | |||
767 | (*current_liboctave_error_handler) ("unrecoverable error in cggev"); | |||
768 | return info; | |||
769 | } | |||
770 | ||||
771 | if (info > 0) | |||
772 | { | |||
773 | (*current_liboctave_error_handler) ("cggev failed to converge"); | |||
774 | return info; | |||
775 | } | |||
776 | ||||
777 | lambda.resize (n); | |||
778 | ||||
779 | for (octave_idx_type j = 0; j < n; j++) | |||
780 | lambda.elem (j) = alpha.elem (j) / beta.elem (j); | |||
781 | ||||
782 | v = vtmp; | |||
783 | } | |||
784 | else | |||
785 | (*current_liboctave_error_handler) ("cggev workspace query failed"); | |||
786 | ||||
787 | return info; | |||
788 | } | |||
789 | ||||
790 | octave_idx_type | |||
791 | FloatEIG::hermitian_init (const FloatComplexMatrix& a, | |||
792 | const FloatComplexMatrix& b, bool calc_ev) | |||
793 | { | |||
794 | octave_idx_type n = a.rows (); | |||
795 | octave_idx_type nb = b.rows (); | |||
796 | ||||
797 | if (n != a.cols () || nb != b.cols ()) | |||
798 | { | |||
799 | (*current_liboctave_error_handler) ("EIG requires square matrix"); | |||
800 | return -1; | |||
801 | } | |||
802 | ||||
803 | if (n != nb) | |||
804 | { | |||
805 | (*current_liboctave_error_handler) ("EIG requires same size matrices"); | |||
806 | return -1; | |||
807 | } | |||
808 | ||||
809 | octave_idx_type info = 0; | |||
810 | ||||
811 | FloatComplexMatrix atmp = a; | |||
812 | FloatComplex *atmp_data = atmp.fortran_vec (); | |||
813 | ||||
814 | FloatComplexMatrix btmp = b; | |||
815 | FloatComplex *btmp_data = btmp.fortran_vec (); | |||
816 | ||||
817 | FloatColumnVector wr (n); | |||
818 | float *pwr = wr.fortran_vec (); | |||
819 | ||||
820 | octave_idx_type lwork = -1; | |||
821 | FloatComplex dummy_work; | |||
822 | ||||
823 | octave_idx_type lrwork = 3*n; | |||
824 | Array<float> rwork (dim_vector (lrwork, 1)); | |||
825 | float *prwork = rwork.fortran_vec (); | |||
826 | ||||
827 | F77_XFCN (chegv, CHEGV, (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" , "chegv_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; chegv_ (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) | |||
828 | 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" , "chegv_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; chegv_ (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) | |||
829 | 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" , "chegv_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; chegv_ (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) | |||
830 | 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" , "chegv_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; chegv_ (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 | 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" , "chegv_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; chegv_ (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 | 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" , "chegv_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; chegv_ (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 | 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" , "chegv_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; chegv_ (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 | 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" , "chegv_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; chegv_ (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 | ||||
836 | if (info == 0) | |||
837 | { | |||
838 | lwork = static_cast<octave_idx_type> (dummy_work.real ()); | |||
839 | Array<FloatComplex> work (dim_vector (lwork, 1)); | |||
840 | FloatComplex *pwork = work.fortran_vec (); | |||
841 | ||||
842 | F77_XFCN (chegv, CHEGV, (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" , "chegv_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; chegv_ (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) | |||
843 | 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" , "chegv_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; chegv_ (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) | |||
844 | 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" , "chegv_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; chegv_ (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) | |||
845 | 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" , "chegv_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; chegv_ (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 | 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" , "chegv_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; chegv_ (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 | 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" , "chegv_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; chegv_ (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 | 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" , "chegv_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; chegv_ (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 | ||||
850 | if (info < 0) | |||
851 | { | |||
852 | (*current_liboctave_error_handler) ("unrecoverable error in zhegv"); | |||
853 | return info; | |||
854 | } | |||
855 | ||||
856 | if (info > 0) | |||
857 | { | |||
858 | (*current_liboctave_error_handler) ("zhegv failed to converge"); | |||
859 | return info; | |||
860 | } | |||
861 | ||||
862 | lambda = FloatComplexColumnVector (wr); | |||
863 | v = calc_ev ? FloatComplexMatrix (atmp) : FloatComplexMatrix (); | |||
864 | } | |||
865 | else | |||
866 | (*current_liboctave_error_handler) ("zhegv workspace query failed"); | |||
867 | ||||
868 | return info; | |||
869 | } |