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