Mercurial > octave-nkf

annotate liboctave/EIG.cc @ 9636:74be4b7273e4 ss-3-3-50

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