Bug Summary

File:liboctave/numeric/fEIG.cc
Location:line 279, column 13
Description:Assigned value is garbage or undefined

Annotated Source Code

1/*
2
3Copyright (C) 1994-2013 John W. Eaton
4
5This file is part of Octave.
6
7Octave is free software; you can redistribute it and/or modify it
8under the terms of the GNU General Public License as published by the
9Free Software Foundation; either version 3 of the License, or (at your
10option) any later version.
11
12Octave is distributed in the hope that it will be useful, but WITHOUT
13ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
14FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
15for more details.
16
17You should have received a copy of the GNU General Public License
18along 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
32extern "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
136octave_idx_type
137FloatEIG::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
249octave_idx_type
250FloatEIG::symmetric_init (const FloatMatrix& a, bool calc_ev)
251{
252 octave_idx_type n = a.rows ();
253
254 if (n != a.cols ())
1
Taking false branch
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;
2
'dummy_work' declared without an initial value
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)
3
Taking true branch
278 {
279 lwork = static_cast<octave_idx_type> (dummy_work);
4
Assigned value is garbage or undefined
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
310octave_idx_type
311FloatEIG::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
394octave_idx_type
395FloatEIG::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
460octave_idx_type
461FloatEIG::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
601octave_idx_type
602FloatEIG::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
677octave_idx_type
678FloatEIG::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
790octave_idx_type
791FloatEIG::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}