Mercurial > hg > octave-nkf

comparison liboctave/dbleSCHUR.cc @ 457:3d4b4f0fa5ba

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 1994-06-06 00:33:33 by jwe] Initial revision
author jwe
date Mon, 06 Jun 1994 00:33:51 +0000
parents
children 883197c5ad75
comparison
equal deleted inserted replaced
456:a1b3aae0fbc3 457:3d4b4f0fa5ba
1 // -*- C++ -*-
2 /*
3
4 Copyright (C) 1992, 1993, 1994 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 #ifdef HAVE_CONFIG_H
25 #include "config.h"
26 #endif
27
28 #if defined (__GNUG__)
29 #pragma implementation
30 #endif
31
32 #include "dbleSCHUR.h"
33 #include "mx-inlines.cc"
34 #include "lo-error.h"
35 #include "f77-uscore.h"
36
37 extern "C"
38 {
39 int F77_FCN (dgeesx) (const char*, const char*,
40 int (*)(double*, double*), const char*,
41 const int*, double*, const int*, int*, double*,
42 double*, double*, const int*, double*, double*,
43 double*, const int*, int*, const int*, int*,
44 int*, long, long);
45 }
46
47 static int
48 select_ana (double *a, double *b)
49 {
50 return (*a < 0.0);
51 }
52
53 static int
54 select_dig (double *a, double *b)
55 {
56 return (hypot (*a, *b) < 1.0);
57 }
58
59 int
60 SCHUR::init (const Matrix& a, const char *ord)
61 {
62 int a_nr = a.rows ();
63 int a_nc = a.cols ();
64 if (a_nr != a_nc)
65 {
66 (*current_liboctave_error_handler) ("SCHUR requires square matrix");
67 return -1;
68 }
69
70 char jobvs = 'V';
71 char sort;
72
73 if (*ord == 'A' || *ord == 'D' || *ord == 'a' || *ord == 'd')
74 sort = 'S';
75 else
76 sort = 'N';
77
78 char sense = 'N';
79
80 int n = a_nc;
81 int lwork = 8 * n;
82 int liwork = 1;
83 int info;
84 int sdim;
85 double rconde;
86 double rcondv;
87
88 double *s = dup (a.data (), a.length ());
89
90 double *wr = new double [n];
91 double *wi = new double [n];
92 double *q = new double [n*n];
93 double *work = new double [lwork];
94
95 // These are not referenced for the non-ordered Schur routine.
96
97 int *iwork = (int *) NULL;
98 int *bwork = (int *) NULL;
99 if (*ord == 'A' || *ord == 'D' || *ord == 'a' || *ord == 'd')
100 {
101 iwork = new int [liwork];
102 bwork = new int [n];
103 }
104
105 if (*ord == 'A' || *ord == 'a')
106 {
107 F77_FCN (dgeesx) (&jobvs, &sort, select_ana, &sense, &n, s, &n,
108 &sdim, wr, wi, q, &n, &rconde, &rcondv, work,
109 &lwork, iwork, &liwork, bwork, &info, 1L, 1L);
110 }
111 else if (*ord == 'D' || *ord == 'd')
112 {
113 F77_FCN (dgeesx) (&jobvs, &sort, select_dig, &sense, &n, s, &n,
114 &sdim, wr, wi, q, &n, &rconde, &rcondv, work,
115 &lwork, iwork, &liwork, bwork, &info, 1L, 1L);
116
117 }
118 else
119 {
120 F77_FCN (dgeesx) (&jobvs, &sort, (void *) 0, &sense, &n, s,
121 &n, &sdim, wr, wi, q, &n, &rconde, &rcondv,
122 work, &lwork, iwork, &liwork, bwork, &info,
123 1L, 1L);
124 }
125
126 schur_mat = Matrix (s, n, n);
127 unitary_mat = Matrix (q, n, n);
128
129 delete [] wr;
130 delete [] wi;
131 delete [] work;
132 delete [] iwork;
133 delete [] bwork;
134
135 return info;
136 }
137
138 ostream&
139 operator << (ostream& os, const SCHUR& a)
140 {
141 os << a.schur_matrix () << "\n";
142 os << a.unitary_matrix () << "\n";
143
144 return os;
145 }
146
147 /*
148 ;;; Local Variables: ***
149 ;;; mode: C++ ***
150 ;;; page-delimiter: "^/\\*" ***
151 ;;; End: ***
152 */