Mercurial > jwe > octave

comparison liboctave/numeric/fEIG.cc @ 21136:7cac4e7458f2

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
maint: clean up code around calls to current_liboctave_error_handler. Remove statements after call to handler that are no longer reachable. Place input validation first and immediately call handler if necessary. Change if/error_handler/else to if/error_handler and re-indent code. * Array-util.cc, Array.cc, CColVector.cc, CDiagMatrix.cc, CMatrix.cc, CNDArray.cc, CRowVector.cc, CSparse.cc, DiagArray2.cc, MArray.cc, PermMatrix.cc, Sparse.cc, Sparse.h, chMatrix.cc, chNDArray.cc, dColVector.cc, dDiagMatrix.cc, dMatrix.cc, dNDArray.cc, dRowVector.cc, dSparse.cc, fCColVector.cc, fCDiagMatrix.cc, fCMatrix.cc, fCNDArray.cc, fCRowVector.cc, fColVector.cc, fDiagMatrix.cc, fMatrix.cc, fNDArray.cc, fRowVector.cc, idx-vector.cc, CmplxAEPBAL.cc, CmplxCHOL.cc, CmplxGEPBAL.cc, CmplxHESS.cc, CmplxLU.cc, CmplxQR.cc, CmplxSCHUR.cc, CmplxSVD.cc, DASPK.cc, EIG.cc, LSODE.cc, Quad.cc, SparseCmplxCHOL.cc, SparseCmplxLU.cc, SparseCmplxQR.cc, SparseQR.cc, SparsedbleCHOL.cc, SparsedbleLU.cc, base-lu.cc, bsxfun-defs.cc, dbleAEPBAL.cc, dbleCHOL.cc, dbleGEPBAL.cc, dbleHESS.cc, dbleLU.cc, dbleQR.cc, dbleSCHUR.cc, dbleSVD.cc, eigs-base.cc, fCmplxAEPBAL.cc, fCmplxCHOL.cc, fCmplxLU.cc, fCmplxQR.cc, fCmplxSCHUR.cc, fEIG.cc, floatAEPBAL.cc, floatCHOL.cc, floatGEPBAL.cc, floatHESS.cc, floatLU.cc, floatQR.cc, floatSCHUR.cc, floatSVD.cc, lo-specfun.cc, oct-fftw.cc, oct-rand.cc, oct-spparms.cc, sparse-base-chol.cc, sparse-dmsolve.cc, file-ops.cc, lo-sysdep.cc, mach-info.cc, oct-env.cc, oct-syscalls.cc, cmd-edit.cc, cmd-hist.cc, data-conv.cc, lo-ieee.cc, lo-regexp.cc, oct-base64.cc, oct-shlib.cc, pathsearch.cc, singleton-cleanup.cc, sparse-util.cc, unwind-prot.cc: Remove statements after call to handler that are no longer reachable. Place input validation first and immediately call handler if necessary. Change if/error_handler/else to if/error_handler and re-indent code.
author Rik <rik@octave.org>
date Sat, 23 Jan 2016 13:52:03 -0800
parents ff904ae0285b
children f7121e111991
comparison
equal deleted inserted replaced
21135:95da3bc8a281 21136:7cac4e7458f2
135 135
136 octave_idx_type 136 octave_idx_type
137 FloatEIG::init (const FloatMatrix& a, bool calc_ev) 137 FloatEIG::init (const FloatMatrix& a, bool calc_ev)
138 { 138 {
139 if (a.any_element_is_inf_or_nan ()) 139 if (a.any_element_is_inf_or_nan ())
140 { 140 (*current_liboctave_error_handler)
141 (*current_liboctave_error_handler) 141 ("EIG: matrix contains Inf or NaN values");
142 ("EIG: matrix contains Inf or NaN values");
143 return -1;
144 }
145 142
146 if (a.is_symmetric ()) 143 if (a.is_symmetric ())
147 return symmetric_init (a, calc_ev); 144 return symmetric_init (a, calc_ev);
148 145
149 octave_idx_type n = a.rows (); 146 octave_idx_type n = a.rows ();
150 147
151 if (n != a.cols ()) 148 if (n != a.cols ())
152 { 149 (*current_liboctave_error_handler) ("EIG requires square matrix");
153 (*current_liboctave_error_handler) ("EIG requires square matrix");
154 return -1;
155 }
156 150
157 octave_idx_type info = 0; 151 octave_idx_type info = 0;
158 152
159 FloatMatrix atmp = a; 153 FloatMatrix atmp = a;
160 float *tmp_data = atmp.fortran_vec (); 154 float *tmp_data = atmp.fortran_vec ();
180 n, tmp_data, n, pwr, pwi, dummy, 174 n, tmp_data, n, pwr, pwi, dummy,
181 idummy, pvr, n, &dummy_work, lwork, info 175 idummy, pvr, n, &dummy_work, lwork, info
182 F77_CHAR_ARG_LEN (1) 176 F77_CHAR_ARG_LEN (1)
183 F77_CHAR_ARG_LEN (1))); 177 F77_CHAR_ARG_LEN (1)));
184 178
185 if (info == 0) 179 if (info != 0)
180 (*current_liboctave_error_handler) ("sgeev workspace query failed");
181
182 lwork = static_cast<octave_idx_type> (dummy_work);
183 Array<float> work (dim_vector (lwork, 1));
184 float *pwork = work.fortran_vec ();
185
186 F77_XFCN (sgeev, SGEEV, (F77_CONST_CHAR_ARG2 ("N", 1),
187 F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
188 n, tmp_data, n, pwr, pwi, dummy,
189 idummy, pvr, n, pwork, lwork, info
190 F77_CHAR_ARG_LEN (1)
191 F77_CHAR_ARG_LEN (1)));
192
193 if (info < 0)
194 (*current_liboctave_error_handler) ("unrecoverable error in sgeev");
195
196 if (info > 0)
197 (*current_liboctave_error_handler) ("sgeev failed to converge");
198
199 lambda.resize (n);
200 v.resize (nvr, nvr);
201
202 for (octave_idx_type j = 0; j < n; j++)
186 { 203 {
187 lwork = static_cast<octave_idx_type> (dummy_work); 204 if (wi.elem (j) == 0.0)
188 Array<float> work (dim_vector (lwork, 1));
189 float *pwork = work.fortran_vec ();
190
191 F77_XFCN (sgeev, SGEEV, (F77_CONST_CHAR_ARG2 ("N", 1),
192 F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
193 n, tmp_data, n, pwr, pwi, dummy,
194 idummy, pvr, n, pwork, lwork, info
195 F77_CHAR_ARG_LEN (1)
196 F77_CHAR_ARG_LEN (1)));
197
198 if (info < 0)
199 { 205 {
200 (*current_liboctave_error_handler) ("unrecoverable error in sgeev"); 206 lambda.elem (j) = FloatComplex (wr.elem (j));
201 return info; 207 for (octave_idx_type i = 0; i < nvr; i++)
208 v.elem (i, j) = vr.elem (i, j);
202 } 209 }
203 210 else
204 if (info > 0)
205 { 211 {
206 (*current_liboctave_error_handler) ("sgeev failed to converge"); 212 if (j+1 >= n)
207 return info; 213 (*current_liboctave_error_handler) ("EIG: internal error");
208 } 214
209 215 lambda.elem (j) = FloatComplex (wr.elem (j), wi.elem (j));
210 lambda.resize (n); 216 lambda.elem (j+1) = FloatComplex (wr.elem (j+1), wi.elem (j+1));
211 v.resize (nvr, nvr); 217
212 218 for (octave_idx_type i = 0; i < nvr; i++)
213 for (octave_idx_type j = 0; j < n; j++)
214 {
215 if (wi.elem (j) == 0.0)
216 { 219 {
217 lambda.elem (j) = FloatComplex (wr.elem (j)); 220 float real_part = vr.elem (i, j);
218 for (octave_idx_type i = 0; i < nvr; i++) 221 float imag_part = vr.elem (i, j+1);
219 v.elem (i, j) = vr.elem (i, j); 222 v.elem (i, j) = FloatComplex (real_part, imag_part);
223 v.elem (i, j+1) = FloatComplex (real_part, -imag_part);
220 } 224 }
221 else 225 j++;
222 {
223 if (j+1 >= n)
224 {
225 (*current_liboctave_error_handler) ("EIG: internal error");
226 return -1;
227 }
228
229 lambda.elem (j) = FloatComplex (wr.elem (j), wi.elem (j));
230 lambda.elem (j+1) = FloatComplex (wr.elem (j+1), wi.elem (j+1));
231
232 for (octave_idx_type i = 0; i < nvr; i++)
233 {
234 float real_part = vr.elem (i, j);
235 float imag_part = vr.elem (i, j+1);
236 v.elem (i, j) = FloatComplex (real_part, imag_part);
237 v.elem (i, j+1) = FloatComplex (real_part, -imag_part);
238 }
239 j++;
240 }
241 } 226 }
242 } 227 }
243 else
244 (*current_liboctave_error_handler) ("sgeev workspace query failed");
245 228
246 return info; 229 return info;
247 } 230 }
248 231
249 octave_idx_type 232 octave_idx_type
250 FloatEIG::symmetric_init (const FloatMatrix& a, bool calc_ev) 233 FloatEIG::symmetric_init (const FloatMatrix& a, bool calc_ev)
251 { 234 {
252 octave_idx_type n = a.rows (); 235 octave_idx_type n = a.rows ();
253 236
254 if (n != a.cols ()) 237 if (n != a.cols ())
255 { 238 (*current_liboctave_error_handler) ("EIG requires square matrix");
256 (*current_liboctave_error_handler) ("EIG requires square matrix");
257 return -1;
258 }
259 239
260 octave_idx_type info = 0; 240 octave_idx_type info = 0;
261 241
262 FloatMatrix atmp = a; 242 FloatMatrix atmp = a;
263 float *tmp_data = atmp.fortran_vec (); 243 float *tmp_data = atmp.fortran_vec ();
272 F77_CONST_CHAR_ARG2 ("U", 1), 252 F77_CONST_CHAR_ARG2 ("U", 1),
273 n, tmp_data, n, pwr, &dummy_work, lwork, info 253 n, tmp_data, n, pwr, &dummy_work, lwork, info
274 F77_CHAR_ARG_LEN (1) 254 F77_CHAR_ARG_LEN (1)
275 F77_CHAR_ARG_LEN (1))); 255 F77_CHAR_ARG_LEN (1)));
276 256
277 if (info == 0) 257 if (info != 0)
278 {
279 lwork = static_cast<octave_idx_type> (dummy_work);
280 Array<float> work (dim_vector (lwork, 1));
281 float *pwork = work.fortran_vec ();
282
283 F77_XFCN (ssyev, SSYEV, (F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
284 F77_CONST_CHAR_ARG2 ("U", 1),
285 n, tmp_data, n, pwr, pwork, lwork, info
286 F77_CHAR_ARG_LEN (1)
287 F77_CHAR_ARG_LEN (1)));
288
289 if (info < 0)
290 {
291 (*current_liboctave_error_handler) ("unrecoverable error in ssyev");
292 return info;
293 }
294
295 if (info > 0)
296 {
297 (*current_liboctave_error_handler) ("ssyev failed to converge");
298 return info;
299 }
300
301 lambda = FloatComplexColumnVector (wr);
302 v = calc_ev ? FloatComplexMatrix (atmp) : FloatComplexMatrix ();
303 }
304 else
305 (*current_liboctave_error_handler) ("ssyev workspace query failed"); 258 (*current_liboctave_error_handler) ("ssyev workspace query failed");
259
260 lwork = static_cast<octave_idx_type> (dummy_work);
261 Array<float> work (dim_vector (lwork, 1));
262 float *pwork = work.fortran_vec ();
263
264 F77_XFCN (ssyev, SSYEV, (F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
265 F77_CONST_CHAR_ARG2 ("U", 1),
266 n, tmp_data, n, pwr, pwork, lwork, info
267 F77_CHAR_ARG_LEN (1)
268 F77_CHAR_ARG_LEN (1)));
269
270 if (info < 0)
271 (*current_liboctave_error_handler) ("unrecoverable error in ssyev");
272
273 if (info > 0)
274 (*current_liboctave_error_handler) ("ssyev failed to converge");
275
276 lambda = FloatComplexColumnVector (wr);
277 v = calc_ev ? FloatComplexMatrix (atmp) : FloatComplexMatrix ();
306 278
307 return info; 279 return info;
308 } 280 }
309 281
310 octave_idx_type 282 octave_idx_type
311 FloatEIG::init (const FloatComplexMatrix& a, bool calc_ev) 283 FloatEIG::init (const FloatComplexMatrix& a, bool calc_ev)
312 { 284 {
313 if (a.any_element_is_inf_or_nan ()) 285 if (a.any_element_is_inf_or_nan ())
314 { 286 (*current_liboctave_error_handler)
315 (*current_liboctave_error_handler) 287 ("EIG: matrix contains Inf or NaN values");
316 ("EIG: matrix contains Inf or NaN values");
317 return -1;
318 }
319 288
320 if (a.is_hermitian ()) 289 if (a.is_hermitian ())
321 return hermitian_init (a, calc_ev); 290 return hermitian_init (a, calc_ev);
322 291
323 octave_idx_type n = a.rows (); 292 octave_idx_type n = a.rows ();
324 293
325 if (n != a.cols ()) 294 if (n != a.cols ())
326 { 295 (*current_liboctave_error_handler) ("EIG requires square matrix");
327 (*current_liboctave_error_handler) ("EIG requires square matrix");
328 return -1;
329 }
330 296
331 octave_idx_type info = 0; 297 octave_idx_type info = 0;
332 298
333 FloatComplexMatrix atmp = a; 299 FloatComplexMatrix atmp = a;
334 FloatComplex *tmp_data = atmp.fortran_vec (); 300 FloatComplex *tmp_data = atmp.fortran_vec ();
355 n, tmp_data, n, pw, dummy, idummy, 321 n, tmp_data, n, pw, dummy, idummy,
356 pv, n, &dummy_work, lwork, prwork, info 322 pv, n, &dummy_work, lwork, prwork, info
357 F77_CHAR_ARG_LEN (1) 323 F77_CHAR_ARG_LEN (1)
358 F77_CHAR_ARG_LEN (1))); 324 F77_CHAR_ARG_LEN (1)));
359 325
360 if (info == 0) 326 if (info != 0)
361 {
362 lwork = static_cast<octave_idx_type> (dummy_work.real ());
363 Array<FloatComplex> work (dim_vector (lwork, 1));
364 FloatComplex *pwork = work.fortran_vec ();
365
366 F77_XFCN (cgeev, CGEEV, (F77_CONST_CHAR_ARG2 ("N", 1),
367 F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
368 n, tmp_data, n, pw, dummy, idummy,
369 pv, n, pwork, lwork, prwork, info
370 F77_CHAR_ARG_LEN (1)
371 F77_CHAR_ARG_LEN (1)));
372
373 if (info < 0)
374 {
375 (*current_liboctave_error_handler) ("unrecoverable error in cgeev");
376 return info;
377 }
378
379 if (info > 0)
380 {
381 (*current_liboctave_error_handler) ("cgeev failed to converge");
382 return info;
383 }
384
385 lambda = w;
386 v = vtmp;
387 }
388 else
389 (*current_liboctave_error_handler) ("cgeev workspace query failed"); 327 (*current_liboctave_error_handler) ("cgeev workspace query failed");
328
329 lwork = static_cast<octave_idx_type> (dummy_work.real ());
330 Array<FloatComplex> work (dim_vector (lwork, 1));
331 FloatComplex *pwork = work.fortran_vec ();
332
333 F77_XFCN (cgeev, CGEEV, (F77_CONST_CHAR_ARG2 ("N", 1),
334 F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
335 n, tmp_data, n, pw, dummy, idummy,
336 pv, n, pwork, lwork, prwork, info
337 F77_CHAR_ARG_LEN (1)
338 F77_CHAR_ARG_LEN (1)));
339
340 if (info < 0)
341 (*current_liboctave_error_handler) ("unrecoverable error in cgeev");
342
343 if (info > 0)
344 (*current_liboctave_error_handler) ("cgeev failed to converge");
345
346 lambda = w;
347 v = vtmp;
390 348
391 return info; 349 return info;
392 } 350 }
393 351
394 octave_idx_type 352 octave_idx_type
395 FloatEIG::hermitian_init (const FloatComplexMatrix& a, bool calc_ev) 353 FloatEIG::hermitian_init (const FloatComplexMatrix& a, bool calc_ev)
396 { 354 {
397 octave_idx_type n = a.rows (); 355 octave_idx_type n = a.rows ();
398 356
399 if (n != a.cols ()) 357 if (n != a.cols ())
400 { 358 (*current_liboctave_error_handler) ("EIG requires square matrix");
401 (*current_liboctave_error_handler) ("EIG requires square matrix");
402 return -1;
403 }
404 359
405 octave_idx_type info = 0; 360 octave_idx_type info = 0;
406 361
407 FloatComplexMatrix atmp = a; 362 FloatComplexMatrix atmp = a;
408 FloatComplex *tmp_data = atmp.fortran_vec (); 363 FloatComplex *tmp_data = atmp.fortran_vec ();
422 n, tmp_data, n, pwr, &dummy_work, lwork, 377 n, tmp_data, n, pwr, &dummy_work, lwork,
423 prwork, info 378 prwork, info
424 F77_CHAR_ARG_LEN (1) 379 F77_CHAR_ARG_LEN (1)
425 F77_CHAR_ARG_LEN (1))); 380 F77_CHAR_ARG_LEN (1)));
426 381
427 if (info == 0) 382 if (info != 0)
428 {
429 lwork = static_cast<octave_idx_type> (dummy_work.real ());
430 Array<FloatComplex> work (dim_vector (lwork, 1));
431 FloatComplex *pwork = work.fortran_vec ();
432
433 F77_XFCN (cheev, CHEEV, (F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
434 F77_CONST_CHAR_ARG2 ("U", 1),
435 n, tmp_data, n, pwr, pwork, lwork, prwork, info
436 F77_CHAR_ARG_LEN (1)
437 F77_CHAR_ARG_LEN (1)));
438
439 if (info < 0)
440 {
441 (*current_liboctave_error_handler) ("unrecoverable error in cheev");
442 return info;
443 }
444
445 if (info > 0)
446 {
447 (*current_liboctave_error_handler) ("cheev failed to converge");
448 return info;
449 }
450
451 lambda = FloatComplexColumnVector (wr);
452 v = calc_ev ? FloatComplexMatrix (atmp) : FloatComplexMatrix ();
453 }
454 else
455 (*current_liboctave_error_handler) ("cheev workspace query failed"); 383 (*current_liboctave_error_handler) ("cheev workspace query failed");
384
385 lwork = static_cast<octave_idx_type> (dummy_work.real ());
386 Array<FloatComplex> work (dim_vector (lwork, 1));
387 FloatComplex *pwork = work.fortran_vec ();
388
389 F77_XFCN (cheev, CHEEV, (F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
390 F77_CONST_CHAR_ARG2 ("U", 1),
391 n, tmp_data, n, pwr, pwork, lwork, prwork, info
392 F77_CHAR_ARG_LEN (1)
393 F77_CHAR_ARG_LEN (1)));
394
395 if (info < 0)
396 (*current_liboctave_error_handler) ("unrecoverable error in cheev");
397
398 if (info > 0)
399 (*current_liboctave_error_handler) ("cheev failed to converge");
400
401 lambda = FloatComplexColumnVector (wr);
402 v = calc_ev ? FloatComplexMatrix (atmp) : FloatComplexMatrix ();
456 403
457 return info; 404 return info;
458 } 405 }
459 406
460 octave_idx_type 407 octave_idx_type
461 FloatEIG::init (const FloatMatrix& a, const FloatMatrix& b, bool calc_ev) 408 FloatEIG::init (const FloatMatrix& a, const FloatMatrix& b, bool calc_ev)
462 { 409 {
463 if (a.any_element_is_inf_or_nan () || b.any_element_is_inf_or_nan ()) 410 if (a.any_element_is_inf_or_nan () || b.any_element_is_inf_or_nan ())
464 { 411 (*current_liboctave_error_handler)
465 (*current_liboctave_error_handler) 412 ("EIG: matrix contains Inf or NaN values");
466 ("EIG: matrix contains Inf or NaN values");
467 return -1;
468 }
469 413
470 octave_idx_type n = a.rows (); 414 octave_idx_type n = a.rows ();
471 octave_idx_type nb = b.rows (); 415 octave_idx_type nb = b.rows ();
472 416
473 if (n != a.cols () || nb != b.cols ()) 417 if (n != a.cols () || nb != b.cols ())
474 { 418 (*current_liboctave_error_handler) ("EIG requires square matrix");
475 (*current_liboctave_error_handler) ("EIG requires square matrix");
476 return -1;
477 }
478 419
479 if (n != nb) 420 if (n != nb)
480 { 421 (*current_liboctave_error_handler) ("EIG requires same size matrices");
481 (*current_liboctave_error_handler) ("EIG requires same size matrices");
482 return -1;
483 }
484 422
485 octave_idx_type info = 0; 423 octave_idx_type info = 0;
486 424
487 FloatMatrix tmp = b; 425 FloatMatrix tmp = b;
488 float *tmp_data = tmp.fortran_vec (); 426 float *tmp_data = tmp.fortran_vec ();
528 dummy, idummy, pvr, n, 466 dummy, idummy, pvr, n,
529 &dummy_work, lwork, info 467 &dummy_work, lwork, info
530 F77_CHAR_ARG_LEN (1) 468 F77_CHAR_ARG_LEN (1)
531 F77_CHAR_ARG_LEN (1))); 469 F77_CHAR_ARG_LEN (1)));
532 470
533 if (info == 0) 471 if (info != 0)
472 (*current_liboctave_error_handler) ("sggev workspace query failed");
473
474 lwork = static_cast<octave_idx_type> (dummy_work);
475 Array<float> work (dim_vector (lwork, 1));
476 float *pwork = work.fortran_vec ();
477
478 F77_XFCN (sggev, SGGEV, (F77_CONST_CHAR_ARG2 ("N", 1),
479 F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
480 n, atmp_data, n, btmp_data, n,
481 par, pai, pbeta,
482 dummy, idummy, pvr, n,
483 pwork, lwork, info
484 F77_CHAR_ARG_LEN (1)
485 F77_CHAR_ARG_LEN (1)));
486
487 if (info < 0)
488 (*current_liboctave_error_handler) ("unrecoverable error in sggev");
489
490 if (info > 0)
491 (*current_liboctave_error_handler) ("sggev failed to converge");
492
493 lambda.resize (n);
494 v.resize (nvr, nvr);
495
496 for (octave_idx_type j = 0; j < n; j++)
534 { 497 {
535 lwork = static_cast<octave_idx_type> (dummy_work); 498 if (ai.elem (j) == 0.0)
536 Array<float> work (dim_vector (lwork, 1));
537 float *pwork = work.fortran_vec ();
538
539 F77_XFCN (sggev, SGGEV, (F77_CONST_CHAR_ARG2 ("N", 1),
540 F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
541 n, atmp_data, n, btmp_data, n,
542 par, pai, pbeta,
543 dummy, idummy, pvr, n,
544 pwork, lwork, info
545 F77_CHAR_ARG_LEN (1)
546 F77_CHAR_ARG_LEN (1)));
547
548 if (info < 0)
549 { 499 {
550 (*current_liboctave_error_handler) ("unrecoverable error in sggev"); 500 lambda.elem (j) = FloatComplex (ar.elem (j) / beta.elem (j));
551 return info; 501 for (octave_idx_type i = 0; i < nvr; i++)
502 v.elem (i, j) = vr.elem (i, j);
552 } 503 }
553 504 else
554 if (info > 0)
555 { 505 {
556 (*current_liboctave_error_handler) ("sggev failed to converge"); 506 if (j+1 >= n)
557 return info; 507 (*current_liboctave_error_handler) ("EIG: internal error");
558 } 508
559 509 lambda.elem (j) = FloatComplex (ar.elem (j) / beta.elem (j),
560 lambda.resize (n); 510 ai.elem (j) / beta.elem (j));
561 v.resize (nvr, nvr); 511 lambda.elem (j+1) = FloatComplex (ar.elem (j+1) / beta.elem (j+1),
562 512 ai.elem (j+1) / beta.elem (j+1));
563 for (octave_idx_type j = 0; j < n; j++) 513
564 { 514 for (octave_idx_type i = 0; i < nvr; i++)
565 if (ai.elem (j) == 0.0)
566 { 515 {
567 lambda.elem (j) = FloatComplex (ar.elem (j) / beta.elem (j)); 516 float real_part = vr.elem (i, j);
568 for (octave_idx_type i = 0; i < nvr; i++) 517 float imag_part = vr.elem (i, j+1);
569 v.elem (i, j) = vr.elem (i, j); 518 v.elem (i, j) = FloatComplex (real_part, imag_part);
519 v.elem (i, j+1) = FloatComplex (real_part, -imag_part);
570 } 520 }
571 else 521 j++;
572 {
573 if (j+1 >= n)
574 {
575 (*current_liboctave_error_handler) ("EIG: internal error");
576 return -1;
577 }
578
579 lambda.elem (j) = FloatComplex (ar.elem (j) / beta.elem (j),
580 ai.elem (j) / beta.elem (j));
581 lambda.elem (j+1) = FloatComplex (ar.elem (j+1) / beta.elem (j+1),
582 ai.elem (j+1) / beta.elem (j+1));
583
584 for (octave_idx_type i = 0; i < nvr; i++)
585 {
586 float real_part = vr.elem (i, j);
587 float imag_part = vr.elem (i, j+1);
588 v.elem (i, j) = FloatComplex (real_part, imag_part);
589 v.elem (i, j+1) = FloatComplex (real_part, -imag_part);
590 }
591 j++;
592 }
593 } 522 }
594 } 523 }
595 else
596 (*current_liboctave_error_handler) ("sggev workspace query failed");
597 524
598 return info; 525 return info;
599 } 526 }
600 527
601 octave_idx_type 528 octave_idx_type
604 { 531 {
605 octave_idx_type n = a.rows (); 532 octave_idx_type n = a.rows ();
606 octave_idx_type nb = b.rows (); 533 octave_idx_type nb = b.rows ();
607 534
608 if (n != a.cols () || nb != b.cols ()) 535 if (n != a.cols () || nb != b.cols ())
609 { 536 (*current_liboctave_error_handler) ("EIG requires square matrix");
610 (*current_liboctave_error_handler) ("EIG requires square matrix");
611 return -1;
612 }
613 537
614 if (n != nb) 538 if (n != nb)
615 { 539 (*current_liboctave_error_handler) ("EIG requires same size matrices");
616 (*current_liboctave_error_handler) ("EIG requires same size matrices");
617 return -1;
618 }
619 540
620 octave_idx_type info = 0; 541 octave_idx_type info = 0;
621 542
622 FloatMatrix atmp = a; 543 FloatMatrix atmp = a;
623 float *atmp_data = atmp.fortran_vec (); 544 float *atmp_data = atmp.fortran_vec ();
637 btmp_data, n, 558 btmp_data, n,
638 pwr, &dummy_work, lwork, info 559 pwr, &dummy_work, lwork, info
639 F77_CHAR_ARG_LEN (1) 560 F77_CHAR_ARG_LEN (1)
640 F77_CHAR_ARG_LEN (1))); 561 F77_CHAR_ARG_LEN (1)));
641 562
642 if (info == 0) 563 if (info != 0)
643 {
644 lwork = static_cast<octave_idx_type> (dummy_work);
645 Array<float> work (dim_vector (lwork, 1));
646 float *pwork = work.fortran_vec ();
647
648 F77_XFCN (ssygv, SSYGV, (1, F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
649 F77_CONST_CHAR_ARG2 ("U", 1),
650 n, atmp_data, n,
651 btmp_data, n,
652 pwr, pwork, lwork, info
653 F77_CHAR_ARG_LEN (1)
654 F77_CHAR_ARG_LEN (1)));
655
656 if (info < 0)
657 {
658 (*current_liboctave_error_handler) ("unrecoverable error in ssygv");
659 return info;
660 }
661
662 if (info > 0)
663 {
664 (*current_liboctave_error_handler) ("ssygv failed to converge");
665 return info;
666 }
667
668 lambda = FloatComplexColumnVector (wr);
669 v = calc_ev ? FloatComplexMatrix (atmp) : FloatComplexMatrix ();
670 }
671 else
672 (*current_liboctave_error_handler) ("ssygv workspace query failed"); 564 (*current_liboctave_error_handler) ("ssygv workspace query failed");
565
566 lwork = static_cast<octave_idx_type> (dummy_work);
567 Array<float> work (dim_vector (lwork, 1));
568 float *pwork = work.fortran_vec ();
569
570 F77_XFCN (ssygv, SSYGV, (1, F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
571 F77_CONST_CHAR_ARG2 ("U", 1),
572 n, atmp_data, n,
573 btmp_data, n,
574 pwr, pwork, lwork, info
575 F77_CHAR_ARG_LEN (1)
576 F77_CHAR_ARG_LEN (1)));
577
578 if (info < 0)
579 (*current_liboctave_error_handler) ("unrecoverable error in ssygv");
580
581 if (info > 0)
582 (*current_liboctave_error_handler) ("ssygv failed to converge");
583
584 lambda = FloatComplexColumnVector (wr);
585 v = calc_ev ? FloatComplexMatrix (atmp) : FloatComplexMatrix ();
673 586
674 return info; 587 return info;
675 } 588 }
676 589
677 octave_idx_type 590 octave_idx_type
678 FloatEIG::init (const FloatComplexMatrix& a, const FloatComplexMatrix& b, 591 FloatEIG::init (const FloatComplexMatrix& a, const FloatComplexMatrix& b,
679 bool calc_ev) 592 bool calc_ev)
680 { 593 {
681 if (a.any_element_is_inf_or_nan () || b.any_element_is_inf_or_nan ()) 594 if (a.any_element_is_inf_or_nan () || b.any_element_is_inf_or_nan ())
682 { 595 (*current_liboctave_error_handler)
683 (*current_liboctave_error_handler) 596 ("EIG: matrix contains Inf or NaN values");
684 ("EIG: matrix contains Inf or NaN values");
685 return -1;
686 }
687 597
688 octave_idx_type n = a.rows (); 598 octave_idx_type n = a.rows ();
689 octave_idx_type nb = b.rows (); 599 octave_idx_type nb = b.rows ();
690 600
691 if (n != a.cols () || nb != b.cols ()) 601 if (n != a.cols () || nb != b.cols ())
692 { 602 (*current_liboctave_error_handler) ("EIG requires square matrix");
693 (*current_liboctave_error_handler) ("EIG requires square matrix");
694 return -1;
695 }
696 603
697 if (n != nb) 604 if (n != nb)
698 { 605 (*current_liboctave_error_handler) ("EIG requires same size matrices");
699 (*current_liboctave_error_handler) ("EIG requires same size matrices");
700 return -1;
701 }
702 606
703 octave_idx_type info = 0; 607 octave_idx_type info = 0;
704 608
705 FloatComplexMatrix tmp = b; 609 FloatComplexMatrix tmp = b;
706 FloatComplex *tmp_data = tmp.fortran_vec (); 610 FloatComplex *tmp_data = tmp.fortran_vec ();
746 palpha, pbeta, dummy, idummy, 650 palpha, pbeta, dummy, idummy,
747 pv, n, &dummy_work, lwork, prwork, info 651 pv, n, &dummy_work, lwork, prwork, info
748 F77_CHAR_ARG_LEN (1) 652 F77_CHAR_ARG_LEN (1)
749 F77_CHAR_ARG_LEN (1))); 653 F77_CHAR_ARG_LEN (1)));
750 654
751 if (info == 0) 655 if (info != 0)
752 {
753 lwork = static_cast<octave_idx_type> (dummy_work.real ());
754 Array<FloatComplex> work (dim_vector (lwork, 1));
755 FloatComplex *pwork = work.fortran_vec ();
756
757 F77_XFCN (cggev, CGGEV, (F77_CONST_CHAR_ARG2 ("N", 1),
758 F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
759 n, atmp_data, n, btmp_data, n,
760 palpha, pbeta, dummy, idummy,
761 pv, n, pwork, lwork, prwork, info
762 F77_CHAR_ARG_LEN (1)
763 F77_CHAR_ARG_LEN (1)));
764
765 if (info < 0)
766 {
767 (*current_liboctave_error_handler) ("unrecoverable error in cggev");
768 return info;
769 }
770
771 if (info > 0)
772 {
773 (*current_liboctave_error_handler) ("cggev failed to converge");
774 return info;
775 }
776
777 lambda.resize (n);
778
779 for (octave_idx_type j = 0; j < n; j++)
780 lambda.elem (j) = alpha.elem (j) / beta.elem (j);
781
782 v = vtmp;
783 }
784 else
785 (*current_liboctave_error_handler) ("cggev workspace query failed"); 656 (*current_liboctave_error_handler) ("cggev workspace query failed");
657
658 lwork = static_cast<octave_idx_type> (dummy_work.real ());
659 Array<FloatComplex> work (dim_vector (lwork, 1));
660 FloatComplex *pwork = work.fortran_vec ();
661
662 F77_XFCN (cggev, CGGEV, (F77_CONST_CHAR_ARG2 ("N", 1),
663 F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
664 n, atmp_data, n, btmp_data, n,
665 palpha, pbeta, dummy, idummy,
666 pv, n, pwork, lwork, prwork, info
667 F77_CHAR_ARG_LEN (1)
668 F77_CHAR_ARG_LEN (1)));
669
670 if (info < 0)
671 (*current_liboctave_error_handler) ("unrecoverable error in cggev");
672
673 if (info > 0)
674 (*current_liboctave_error_handler) ("cggev failed to converge");
675
676 lambda.resize (n);
677
678 for (octave_idx_type j = 0; j < n; j++)
679 lambda.elem (j) = alpha.elem (j) / beta.elem (j);
680
681 v = vtmp;
786 682
787 return info; 683 return info;
788 } 684 }
789 685
790 octave_idx_type 686 octave_idx_type
793 { 689 {
794 octave_idx_type n = a.rows (); 690 octave_idx_type n = a.rows ();
795 octave_idx_type nb = b.rows (); 691 octave_idx_type nb = b.rows ();
796 692
797 if (n != a.cols () || nb != b.cols ()) 693 if (n != a.cols () || nb != b.cols ())
798 { 694 (*current_liboctave_error_handler) ("EIG requires square matrix");
799 (*current_liboctave_error_handler) ("EIG requires square matrix");
800 return -1;
801 }
802 695
803 if (n != nb) 696 if (n != nb)
804 { 697 (*current_liboctave_error_handler) ("EIG requires same size matrices");
805 (*current_liboctave_error_handler) ("EIG requires same size matrices");
806 return -1;
807 }
808 698
809 octave_idx_type info = 0; 699 octave_idx_type info = 0;
810 700
811 FloatComplexMatrix atmp = a; 701 FloatComplexMatrix atmp = a;
812 FloatComplex *atmp_data = atmp.fortran_vec (); 702 FloatComplex *atmp_data = atmp.fortran_vec ();
831 pwr, &dummy_work, lwork, 721 pwr, &dummy_work, lwork,
832 prwork, info 722 prwork, info
833 F77_CHAR_ARG_LEN (1) 723 F77_CHAR_ARG_LEN (1)
834 F77_CHAR_ARG_LEN (1))); 724 F77_CHAR_ARG_LEN (1)));
835 725
836 if (info == 0) 726 if (info != 0)
837 {
838 lwork = static_cast<octave_idx_type> (dummy_work.real ());
839 Array<FloatComplex> work (dim_vector (lwork, 1));
840 FloatComplex *pwork = work.fortran_vec ();
841
842 F77_XFCN (chegv, CHEGV, (1, F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
843 F77_CONST_CHAR_ARG2 ("U", 1),
844 n, atmp_data, n,
845 btmp_data, n,
846 pwr, pwork, lwork, prwork, info
847 F77_CHAR_ARG_LEN (1)
848 F77_CHAR_ARG_LEN (1)));
849
850 if (info < 0)
851 {
852 (*current_liboctave_error_handler) ("unrecoverable error in zhegv");
853 return info;
854 }
855
856 if (info > 0)
857 {
858 (*current_liboctave_error_handler) ("zhegv failed to converge");
859 return info;
860 }
861
862 lambda = FloatComplexColumnVector (wr);
863 v = calc_ev ? FloatComplexMatrix (atmp) : FloatComplexMatrix ();
864 }
865 else
866 (*current_liboctave_error_handler) ("zhegv workspace query failed"); 727 (*current_liboctave_error_handler) ("zhegv workspace query failed");
867 728
729 lwork = static_cast<octave_idx_type> (dummy_work.real ());
730 Array<FloatComplex> work (dim_vector (lwork, 1));
731 FloatComplex *pwork = work.fortran_vec ();
732
733 F77_XFCN (chegv, CHEGV, (1, F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
734 F77_CONST_CHAR_ARG2 ("U", 1),
735 n, atmp_data, n,
736 btmp_data, n,
737 pwr, pwork, lwork, prwork, info
738 F77_CHAR_ARG_LEN (1)
739 F77_CHAR_ARG_LEN (1)));
740
741 if (info < 0)
742 (*current_liboctave_error_handler) ("unrecoverable error in zhegv");
743
744 if (info > 0)
745 (*current_liboctave_error_handler) ("zhegv failed to converge");
746
747 lambda = FloatComplexColumnVector (wr);
748 v = calc_ev ? FloatComplexMatrix (atmp) : FloatComplexMatrix ();
749
868 return info; 750 return info;
869 } 751 }