octave-lyh: liboctave/numeric/EIG.cc annotate

annotate liboctave/numeric/EIG.cc @ 17535:c12c688a35ed default tip lyh

Fix warnings

author	LYH <lyh.kernel@gmail.com>
date	Fri, 27 Sep 2013 17:43:27 +0800
parents	648dabbb4c6b
children

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

Mercurial > hg > octave-lyh

annotate liboctave/numeric/EIG.cc @ 17535:c12c688a35ed default tip lyh