Mercurial > octave-nkf

comparison src/qzval.cc @ 44:f3215b07c171

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 1993-08-10 22:33:20 by jwe] Initial revision
author jwe
date Tue, 10 Aug 1993 22:33:20 +0000
parents
children ed620db95182
comparison
equal deleted inserted replaced
43:e399beacf758 44:f3215b07c171
1 // tc-qzval.cc -*- C++ -*-
2 /*
3
4 Copyright (C) 1993 John W. Eaton
5
6 This file is part of Octave.
7
8 Octave is free software; you can redistribute it and/or modify it
9 under the terms of the GNU General Public License as published by the
10 Free Software Foundation; either version 2, or (at your option) any
11 later version.
12
13 Octave is distributed in the hope that it will be useful, but WITHOUT
14 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
15 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
16 for more details.
17
18 You should have received a copy of the GNU General Public License
19 along with Octave; see the file COPYING. If not, write to the Free
20 Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
21
22 */
23
24 // Written by A. S. Hodel <scotte@eng.auburn.edu>
25
26 #ifdef __GNUG__
27 #pragma implementation
28 #endif
29
30 #include "Matrix.h"
31
32 #include "tree-const.h"
33 #include "user-prefs.h"
34 #include "gripes.h"
35 #include "error.h"
36
37 int F77_FCN (qzhes) (const int*, const int*, double*, double*, const
38 long*, double*);
39
40 int F77_FCN (qzit) (const int*, const int*, double*, double*, const
41 double*, const long*, double*, int*);
42
43 int F77_FCN (qzval) (const int*, const int*, double*, double*,
44 double*, double*, double*, const long*, double*);
45
46 // XXX FIXME XXX
47 extern int empty_arg (tree_constant&);
48 extern tree_constant* empty_tree (int , char*);
49 extern ComplexMatrix ComplexMatrixLoad (tree_constant&);
50
51 #ifdef WITH_DLD
52 tree_constant *
53 builtin_qzvalue_2 (tree_constant *args, int nargin, int nargout)
54 {
55 return qzvalue (args, nargin, nargout);
56 }
57 #endif
58
59 tree_constant *
60 qzvalue (tree_constant *args, int nargin, int nargout)
61 {
62 tree_constant *retval = NULL_TREE_CONST;
63
64 tree_constant arga = args[1].make_numeric ();
65 tree_constant argb = args[2].make_numeric();
66
67 if (empty_arg (arga) || empty_arg (argb))
68 retval = empty_tree (nargout, "qzvalue");
69 else
70 {
71
72 // Arguments are not empty, so check for correct dimensions.
73
74 int a_rows = arga.rows();
75 int a_cols = arga.columns();
76 int b_rows = argb.rows();
77 int b_cols = argb.columns();
78
79 if ((a_rows != a_cols) || (b_rows != b_cols))
80 {
81 gripe_square_matrix_required ("qzvalue: first two parameters:");
82 return retval;
83 }
84 else if (a_rows != b_rows)
85 {
86 gripe_nonconformant ();
87 return retval;
88 }
89
90 // Dimensions look o.k., let's solve the problem.
91
92 retval = new tree_constant[nargout+1];
93
94 if (arga.is_complex_type () || argb.is_complex_type ())
95 error ("qzvalue: cannot yet do complex matrix arguments\n");
96 else
97 {
98
99 // Do everything in real arithmetic.
100
101 Matrix jnk (a_rows, a_rows, 0.0);
102
103 ColumnVector alfr (a_rows);
104 ColumnVector alfi (a_rows);
105 ColumnVector beta (a_rows);
106
107 long matz = 0;
108 int info;
109
110 // XXX FIXME ??? XXX
111 double eps = DBL_EPSILON;
112
113 Matrix ca = arga.to_matrix ();
114 Matrix cb = argb.to_matrix ();
115
116 // Use EISPACK qz functions.
117
118 F77_FCN (qzhes) (&a_rows, &a_rows, ca.fortran_vec (),
119 cb.fortran_vec (), &matz, jnk.fortran_vec ());
120
121 F77_FCN (qzit) (&a_rows, &a_rows, ca.fortran_vec (),
122 cb.fortran_vec (), &eps, &matz,
123 jnk.fortran_vec (), &info);
124
125 if (info)
126 error ("qzvalue: trouble in qzit, info = %d", info);
127
128 F77_FCN (qzval) (&a_rows, &a_rows, ca.fortran_vec (),
129 cb.fortran_vec (), alfr.fortran_vec (),
130 alfi.fortran_vec (), beta.fortran_vec (),
131 &matz, jnk.fortran_vec ());
132
133 // Count and extract finite generalized eigenvalues.
134
135 int i, cnt;
136 Complex Im (0, 1);
137 for (i = 0, cnt = 0; i < a_rows; i++)
138 if (beta (i) != 0)
139 cnt++;
140
141 ComplexColumnVector cx (cnt, 0.0);
142
143 for (i = 0; i < a_rows; i++)
144 {
145 if (beta (i) != 0)
146 {
147
148 // Finite generalized eigenvalue.
149
150 cnt--;
151 cx (cnt) = (alfr (i) + Im * alfi (i)) / beta (i);
152 }
153 }
154 retval[0] = tree_constant (cx);
155 }
156 }
157 return retval;
158 }
159
160 /*
161 ;;; Local Variables: ***
162 ;;; mode: C++ ***
163 ;;; page-delimiter: "^/\\*" ***
164 ;;; End: ***
165 */