Bug Summary

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