Mercurial > octave
annotate src/DLD-FUNCTIONS/schur.cc @ 10822:23d2378512a0
implement rsf2csf
author | Jaroslav Hajek <highegg@gmail.com> |
---|---|
date | Tue, 27 Jul 2010 11:26:43 +0200 |
parents | 9c9e07f5eb1c |
children | 1e6664326d32 |
rev | line source |
---|---|
2928 | 1 /* |
2 | |
8920 | 3 Copyright (C) 1996, 1997, 1999, 2000, 2004, 2005, 2006, 2007, 2008, 2009 |
7017 | 4 John W. Eaton |
2928 | 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 | |
7016 | 10 Free Software Foundation; either version 3 of the License, or (at your |
11 option) any later version. | |
2928 | 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 | |
7016 | 19 along with Octave; see the file COPYING. If not, see |
20 <http://www.gnu.org/licenses/>. | |
2928 | 21 |
22 */ | |
23 | |
24 #ifdef HAVE_CONFIG_H | |
25 #include <config.h> | |
26 #endif | |
27 | |
28 #include <string> | |
29 | |
30 #include "CmplxSCHUR.h" | |
31 #include "dbleSCHUR.h" | |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
32 #include "fCmplxSCHUR.h" |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
33 #include "floatSCHUR.h" |
2928 | 34 |
35 #include "defun-dld.h" | |
36 #include "error.h" | |
37 #include "gripes.h" | |
38 #include "oct-obj.h" | |
39 #include "utils.h" | |
40 | |
10617
9c9e07f5eb1c
make schur mark triangular matrices on output
Jaroslav Hajek <highegg@gmail.com>
parents:
10607
diff
changeset
|
41 template <class Matrix> |
9c9e07f5eb1c
make schur mark triangular matrices on output
Jaroslav Hajek <highegg@gmail.com>
parents:
10607
diff
changeset
|
42 static octave_value |
9c9e07f5eb1c
make schur mark triangular matrices on output
Jaroslav Hajek <highegg@gmail.com>
parents:
10607
diff
changeset
|
43 mark_upper_triangular (const Matrix& a) |
9c9e07f5eb1c
make schur mark triangular matrices on output
Jaroslav Hajek <highegg@gmail.com>
parents:
10607
diff
changeset
|
44 { |
9c9e07f5eb1c
make schur mark triangular matrices on output
Jaroslav Hajek <highegg@gmail.com>
parents:
10607
diff
changeset
|
45 octave_value retval = a; |
9c9e07f5eb1c
make schur mark triangular matrices on output
Jaroslav Hajek <highegg@gmail.com>
parents:
10607
diff
changeset
|
46 |
9c9e07f5eb1c
make schur mark triangular matrices on output
Jaroslav Hajek <highegg@gmail.com>
parents:
10607
diff
changeset
|
47 octave_idx_type n = a.rows (); |
9c9e07f5eb1c
make schur mark triangular matrices on output
Jaroslav Hajek <highegg@gmail.com>
parents:
10607
diff
changeset
|
48 assert (a.columns () == n); |
9c9e07f5eb1c
make schur mark triangular matrices on output
Jaroslav Hajek <highegg@gmail.com>
parents:
10607
diff
changeset
|
49 |
9c9e07f5eb1c
make schur mark triangular matrices on output
Jaroslav Hajek <highegg@gmail.com>
parents:
10607
diff
changeset
|
50 const typename Matrix::element_type zero = typename Matrix::element_type (); |
9c9e07f5eb1c
make schur mark triangular matrices on output
Jaroslav Hajek <highegg@gmail.com>
parents:
10607
diff
changeset
|
51 |
9c9e07f5eb1c
make schur mark triangular matrices on output
Jaroslav Hajek <highegg@gmail.com>
parents:
10607
diff
changeset
|
52 for (octave_idx_type i = 0; i < n; i++) |
9c9e07f5eb1c
make schur mark triangular matrices on output
Jaroslav Hajek <highegg@gmail.com>
parents:
10607
diff
changeset
|
53 if (a(i,i) == zero) |
9c9e07f5eb1c
make schur mark triangular matrices on output
Jaroslav Hajek <highegg@gmail.com>
parents:
10607
diff
changeset
|
54 return retval; |
9c9e07f5eb1c
make schur mark triangular matrices on output
Jaroslav Hajek <highegg@gmail.com>
parents:
10607
diff
changeset
|
55 |
9c9e07f5eb1c
make schur mark triangular matrices on output
Jaroslav Hajek <highegg@gmail.com>
parents:
10607
diff
changeset
|
56 retval.matrix_type (MatrixType::Upper); |
9c9e07f5eb1c
make schur mark triangular matrices on output
Jaroslav Hajek <highegg@gmail.com>
parents:
10607
diff
changeset
|
57 |
9c9e07f5eb1c
make schur mark triangular matrices on output
Jaroslav Hajek <highegg@gmail.com>
parents:
10607
diff
changeset
|
58 return retval; |
9c9e07f5eb1c
make schur mark triangular matrices on output
Jaroslav Hajek <highegg@gmail.com>
parents:
10607
diff
changeset
|
59 } |
9c9e07f5eb1c
make schur mark triangular matrices on output
Jaroslav Hajek <highegg@gmail.com>
parents:
10607
diff
changeset
|
60 |
2928 | 61 DEFUN_DLD (schur, args, nargout, |
3548 | 62 "-*- texinfo -*-\n\ |
3372 | 63 @deftypefn {Loadable Function} {@var{s} =} schur (@var{a})\n\ |
10607
f7501986e42d
make schur more Matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
10155
diff
changeset
|
64 @deftypefnx {Loadable Function} {@var{s} =} schur (@var{a}, \"complex\")\n\ |
3372 | 65 @deftypefnx {Loadable Function} {[@var{u}, @var{s}] =} schur (@var{a}, @var{opt})\n\ |
66 @cindex Schur decomposition\n\ | |
67 The Schur decomposition is used to compute eigenvalues of a\n\ | |
68 square matrix, and has applications in the solution of algebraic\n\ | |
69 Riccati equations in control (see @code{are} and @code{dare}).\n\ | |
70 @code{schur} always returns\n\ | |
71 @tex\n\ | |
72 $S = U^T A U$\n\ | |
73 @end tex\n\ | |
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
7814
diff
changeset
|
74 @ifnottex\n\ |
3372 | 75 @code{s = u' * a * u}\n\ |
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
7814
diff
changeset
|
76 @end ifnottex\n\ |
3372 | 77 where\n\ |
78 @tex\n\ | |
79 $U$\n\ | |
80 @end tex\n\ | |
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
7814
diff
changeset
|
81 @ifnottex\n\ |
3372 | 82 @code{u}\n\ |
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
7814
diff
changeset
|
83 @end ifnottex\n\ |
3372 | 84 is a unitary matrix\n\ |
85 @tex\n\ | |
86 ($U^T U$ is identity)\n\ | |
87 @end tex\n\ | |
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
7814
diff
changeset
|
88 @ifnottex\n\ |
3372 | 89 (@code{u'* u} is identity)\n\ |
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
7814
diff
changeset
|
90 @end ifnottex\n\ |
3372 | 91 and\n\ |
92 @tex\n\ | |
93 $S$\n\ | |
94 @end tex\n\ | |
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
7814
diff
changeset
|
95 @ifnottex\n\ |
3372 | 96 @code{s}\n\ |
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
7814
diff
changeset
|
97 @end ifnottex\n\ |
3372 | 98 is upper triangular. The eigenvalues of\n\ |
99 @tex\n\ | |
100 $A$ (and $S$)\n\ | |
101 @end tex\n\ | |
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
7814
diff
changeset
|
102 @ifnottex\n\ |
3372 | 103 @code{a} (and @code{s})\n\ |
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
7814
diff
changeset
|
104 @end ifnottex\n\ |
3372 | 105 are the diagonal elements of\n\ |
106 @tex\n\ | |
5555 | 107 $S$.\n\ |
3372 | 108 @end tex\n\ |
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
7814
diff
changeset
|
109 @ifnottex\n\ |
5555 | 110 @code{s}.\n\ |
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
7814
diff
changeset
|
111 @end ifnottex\n\ |
3372 | 112 If the matrix\n\ |
113 @tex\n\ | |
114 $A$\n\ | |
115 @end tex\n\ | |
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
7814
diff
changeset
|
116 @ifnottex\n\ |
3372 | 117 @code{a}\n\ |
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
7814
diff
changeset
|
118 @end ifnottex\n\ |
3372 | 119 is real, then the real Schur decomposition is computed, in which the\n\ |
120 matrix\n\ | |
121 @tex\n\ | |
122 $U$\n\ | |
123 @end tex\n\ | |
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
7814
diff
changeset
|
124 @ifnottex\n\ |
3372 | 125 @code{u}\n\ |
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
7814
diff
changeset
|
126 @end ifnottex\n\ |
3372 | 127 is orthogonal and\n\ |
128 @tex\n\ | |
129 $S$\n\ | |
130 @end tex\n\ | |
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
7814
diff
changeset
|
131 @ifnottex\n\ |
3372 | 132 @code{s}\n\ |
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
7814
diff
changeset
|
133 @end ifnottex\n\ |
3372 | 134 is block upper triangular\n\ |
135 with blocks of size at most\n\ | |
136 @tex\n\ | |
137 $2\\times 2$\n\ | |
138 @end tex\n\ | |
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
7814
diff
changeset
|
139 @ifnottex\n\ |
3372 | 140 @code{2 x 2}\n\ |
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
7814
diff
changeset
|
141 @end ifnottex\n\ |
3600 | 142 along the diagonal. The diagonal elements of\n\ |
3372 | 143 @tex\n\ |
144 $S$\n\ | |
145 @end tex\n\ | |
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
7814
diff
changeset
|
146 @ifnottex\n\ |
3372 | 147 @code{s}\n\ |
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
7814
diff
changeset
|
148 @end ifnottex\n\ |
3372 | 149 (or the eigenvalues of the\n\ |
150 @tex\n\ | |
151 $2\\times 2$\n\ | |
152 @end tex\n\ | |
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
7814
diff
changeset
|
153 @ifnottex\n\ |
3372 | 154 @code{2 x 2}\n\ |
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
7814
diff
changeset
|
155 @end ifnottex\n\ |
3372 | 156 blocks, when\n\ |
157 appropriate) are the eigenvalues of\n\ | |
158 @tex\n\ | |
159 $A$\n\ | |
160 @end tex\n\ | |
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
7814
diff
changeset
|
161 @ifnottex\n\ |
3372 | 162 @code{a}\n\ |
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
7814
diff
changeset
|
163 @end ifnottex\n\ |
3372 | 164 and\n\ |
165 @tex\n\ | |
166 $S$.\n\ | |
167 @end tex\n\ | |
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
7814
diff
changeset
|
168 @ifnottex\n\ |
3372 | 169 @code{s}.\n\ |
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
7814
diff
changeset
|
170 @end ifnottex\n\ |
2928 | 171 \n\ |
10607
f7501986e42d
make schur more Matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
10155
diff
changeset
|
172 A complex decomposition may be forced by passing \"complex\" as @var{opt}.\n\ |
f7501986e42d
make schur more Matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
10155
diff
changeset
|
173 \n\ |
3372 | 174 The eigenvalues are optionally ordered along the diagonal according to\n\ |
175 the value of @code{opt}. @code{opt = \"a\"} indicates that all\n\ | |
176 eigenvalues with negative real parts should be moved to the leading\n\ | |
177 block of\n\ | |
178 @tex\n\ | |
179 $S$\n\ | |
180 @end tex\n\ | |
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
7814
diff
changeset
|
181 @ifnottex\n\ |
3372 | 182 @code{s}\n\ |
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
7814
diff
changeset
|
183 @end ifnottex\n\ |
3372 | 184 (used in @code{are}), @code{opt = \"d\"} indicates that all eigenvalues\n\ |
185 with magnitude less than one should be moved to the leading block of\n\ | |
186 @tex\n\ | |
187 $S$\n\ | |
188 @end tex\n\ | |
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
7814
diff
changeset
|
189 @ifnottex\n\ |
3372 | 190 @code{s}\n\ |
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
7814
diff
changeset
|
191 @end ifnottex\n\ |
3372 | 192 (used in @code{dare}), and @code{opt = \"u\"}, the default, indicates that\n\ |
193 no ordering of eigenvalues should occur. The leading\n\ | |
194 @tex\n\ | |
195 $k$\n\ | |
196 @end tex\n\ | |
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
7814
diff
changeset
|
197 @ifnottex\n\ |
3372 | 198 @code{k}\n\ |
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
7814
diff
changeset
|
199 @end ifnottex\n\ |
3372 | 200 columns of\n\ |
201 @tex\n\ | |
202 $U$\n\ | |
203 @end tex\n\ | |
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
7814
diff
changeset
|
204 @ifnottex\n\ |
3372 | 205 @code{u}\n\ |
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
7814
diff
changeset
|
206 @end ifnottex\n\ |
3372 | 207 always span the\n\ |
208 @tex\n\ | |
209 $A$-invariant\n\ | |
210 @end tex\n\ | |
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
7814
diff
changeset
|
211 @ifnottex\n\ |
3372 | 212 @code{a}-invariant\n\ |
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
7814
diff
changeset
|
213 @end ifnottex\n\ |
3372 | 214 subspace corresponding to the\n\ |
215 @tex\n\ | |
216 $k$\n\ | |
217 @end tex\n\ | |
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
7814
diff
changeset
|
218 @ifnottex\n\ |
3372 | 219 @code{k}\n\ |
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
7814
diff
changeset
|
220 @end ifnottex\n\ |
3372 | 221 leading eigenvalues of\n\ |
222 @tex\n\ | |
223 $S$.\n\ | |
224 @end tex\n\ | |
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
7814
diff
changeset
|
225 @ifnottex\n\ |
3372 | 226 @code{s}.\n\ |
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
7814
diff
changeset
|
227 @end ifnottex\n\ |
3372 | 228 @end deftypefn") |
2928 | 229 { |
230 octave_value_list retval; | |
231 | |
232 int nargin = args.length (); | |
233 | |
234 if (nargin < 1 || nargin > 2 || nargout > 2) | |
235 { | |
5823 | 236 print_usage (); |
2928 | 237 return retval; |
238 } | |
239 | |
240 octave_value arg = args(0); | |
241 | |
3523 | 242 std::string ord; |
2928 | 243 |
244 if (nargin == 2) | |
245 { | |
246 ord = args(1).string_value (); | |
247 | |
248 if (error_state) | |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9758
diff
changeset
|
249 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9758
diff
changeset
|
250 error ("schur: expecting string as second argument"); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9758
diff
changeset
|
251 return retval; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9758
diff
changeset
|
252 } |
2928 | 253 } |
254 | |
10607
f7501986e42d
make schur more Matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
10155
diff
changeset
|
255 bool force_complex = false; |
2928 | 256 |
10607
f7501986e42d
make schur more Matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
10155
diff
changeset
|
257 if (ord == "complex") |
f7501986e42d
make schur more Matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
10155
diff
changeset
|
258 { |
f7501986e42d
make schur more Matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
10155
diff
changeset
|
259 force_complex = true; |
f7501986e42d
make schur more Matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
10155
diff
changeset
|
260 ord = std::string (); |
f7501986e42d
make schur more Matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
10155
diff
changeset
|
261 } |
f7501986e42d
make schur more Matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
10155
diff
changeset
|
262 else |
2928 | 263 { |
10607
f7501986e42d
make schur more Matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
10155
diff
changeset
|
264 char ord_char = ord.empty () ? 'U' : ord[0]; |
f7501986e42d
make schur more Matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
10155
diff
changeset
|
265 |
f7501986e42d
make schur more Matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
10155
diff
changeset
|
266 if (ord_char != 'U' && ord_char != 'A' && ord_char != 'D' |
f7501986e42d
make schur more Matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
10155
diff
changeset
|
267 && ord_char != 'u' && ord_char != 'a' && ord_char != 'd') |
f7501986e42d
make schur more Matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
10155
diff
changeset
|
268 { |
f7501986e42d
make schur more Matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
10155
diff
changeset
|
269 warning ("schur: incorrect ordered schur argument `%c'", |
f7501986e42d
make schur more Matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
10155
diff
changeset
|
270 ord.c_str ()); |
f7501986e42d
make schur more Matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
10155
diff
changeset
|
271 return retval; |
f7501986e42d
make schur more Matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
10155
diff
changeset
|
272 } |
2928 | 273 } |
274 | |
5275 | 275 octave_idx_type nr = arg.rows (); |
276 octave_idx_type nc = arg.columns (); | |
2928 | 277 |
278 if (nr != nc) | |
279 { | |
280 gripe_square_matrix_required ("schur"); | |
281 return retval; | |
282 } | |
283 | |
10607
f7501986e42d
make schur more Matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
10155
diff
changeset
|
284 if (! arg.is_numeric_type ()) |
f7501986e42d
make schur more Matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
10155
diff
changeset
|
285 gripe_wrong_type_arg ("schur", arg); |
f7501986e42d
make schur more Matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
10155
diff
changeset
|
286 else if (arg.is_single_type ()) |
2928 | 287 { |
10607
f7501986e42d
make schur more Matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
10155
diff
changeset
|
288 if (! force_complex && arg.is_real_type ()) |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9758
diff
changeset
|
289 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9758
diff
changeset
|
290 FloatMatrix tmp = arg.float_matrix_value (); |
2928 | 291 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9758
diff
changeset
|
292 if (! error_state) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9758
diff
changeset
|
293 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9758
diff
changeset
|
294 if (nargout == 0 || nargout == 1) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9758
diff
changeset
|
295 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9758
diff
changeset
|
296 FloatSCHUR result (tmp, ord, false); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9758
diff
changeset
|
297 retval(0) = result.schur_matrix (); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9758
diff
changeset
|
298 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9758
diff
changeset
|
299 else |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9758
diff
changeset
|
300 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9758
diff
changeset
|
301 FloatSCHUR result (tmp, ord, true); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9758
diff
changeset
|
302 retval(1) = result.schur_matrix (); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9758
diff
changeset
|
303 retval(0) = result.unitary_matrix (); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9758
diff
changeset
|
304 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9758
diff
changeset
|
305 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9758
diff
changeset
|
306 } |
10607
f7501986e42d
make schur more Matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
10155
diff
changeset
|
307 else |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9758
diff
changeset
|
308 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9758
diff
changeset
|
309 FloatComplexMatrix ctmp = arg.float_complex_matrix_value (); |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
310 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9758
diff
changeset
|
311 if (! error_state) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9758
diff
changeset
|
312 { |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
313 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9758
diff
changeset
|
314 if (nargout == 0 || nargout == 1) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9758
diff
changeset
|
315 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9758
diff
changeset
|
316 FloatComplexSCHUR result (ctmp, ord, false); |
10617
9c9e07f5eb1c
make schur mark triangular matrices on output
Jaroslav Hajek <highegg@gmail.com>
parents:
10607
diff
changeset
|
317 retval(0) = mark_upper_triangular (result.schur_matrix ()); |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9758
diff
changeset
|
318 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9758
diff
changeset
|
319 else |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9758
diff
changeset
|
320 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9758
diff
changeset
|
321 FloatComplexSCHUR result (ctmp, ord, true); |
10617
9c9e07f5eb1c
make schur mark triangular matrices on output
Jaroslav Hajek <highegg@gmail.com>
parents:
10607
diff
changeset
|
322 retval(1) = mark_upper_triangular (result.schur_matrix ()); |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9758
diff
changeset
|
323 retval(0) = result.unitary_matrix (); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9758
diff
changeset
|
324 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9758
diff
changeset
|
325 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9758
diff
changeset
|
326 } |
2928 | 327 } |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
328 else |
2928 | 329 { |
10607
f7501986e42d
make schur more Matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
10155
diff
changeset
|
330 if (! force_complex && arg.is_real_type ()) |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9758
diff
changeset
|
331 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9758
diff
changeset
|
332 Matrix tmp = arg.matrix_value (); |
2928 | 333 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9758
diff
changeset
|
334 if (! error_state) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9758
diff
changeset
|
335 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9758
diff
changeset
|
336 if (nargout == 0 || nargout == 1) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9758
diff
changeset
|
337 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9758
diff
changeset
|
338 SCHUR result (tmp, ord, false); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9758
diff
changeset
|
339 retval(0) = result.schur_matrix (); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9758
diff
changeset
|
340 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9758
diff
changeset
|
341 else |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9758
diff
changeset
|
342 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9758
diff
changeset
|
343 SCHUR result (tmp, ord, true); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9758
diff
changeset
|
344 retval(1) = result.schur_matrix (); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9758
diff
changeset
|
345 retval(0) = result.unitary_matrix (); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9758
diff
changeset
|
346 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9758
diff
changeset
|
347 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9758
diff
changeset
|
348 } |
10607
f7501986e42d
make schur more Matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
10155
diff
changeset
|
349 else |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9758
diff
changeset
|
350 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9758
diff
changeset
|
351 ComplexMatrix ctmp = arg.complex_matrix_value (); |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
352 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9758
diff
changeset
|
353 if (! error_state) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9758
diff
changeset
|
354 { |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
355 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9758
diff
changeset
|
356 if (nargout == 0 || nargout == 1) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9758
diff
changeset
|
357 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9758
diff
changeset
|
358 ComplexSCHUR result (ctmp, ord, false); |
10617
9c9e07f5eb1c
make schur mark triangular matrices on output
Jaroslav Hajek <highegg@gmail.com>
parents:
10607
diff
changeset
|
359 retval(0) = mark_upper_triangular (result.schur_matrix ()); |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9758
diff
changeset
|
360 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9758
diff
changeset
|
361 else |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9758
diff
changeset
|
362 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9758
diff
changeset
|
363 ComplexSCHUR result (ctmp, ord, true); |
10617
9c9e07f5eb1c
make schur mark triangular matrices on output
Jaroslav Hajek <highegg@gmail.com>
parents:
10607
diff
changeset
|
364 retval(1) = mark_upper_triangular (result.schur_matrix ()); |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9758
diff
changeset
|
365 retval(0) = result.unitary_matrix (); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9758
diff
changeset
|
366 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9758
diff
changeset
|
367 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9758
diff
changeset
|
368 } |
2928 | 369 } |
370 | |
371 return retval; | |
372 } | |
373 | |
374 /* | |
7814
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
375 |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
376 %!test |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
377 %! a = [1, 2, 3; 4, 5, 9; 7, 8, 6]; |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
378 %! [u, s] = schur (a); |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
379 %! assert(u' * a * u, s, sqrt (eps)); |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
380 |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
381 %!test |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
382 %! a = single([1, 2, 3; 4, 5, 9; 7, 8, 6]); |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
383 %! [u, s] = schur (a); |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
384 %! assert(u' * a * u, s, sqrt (eps('single'))); |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
385 |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
386 %!test |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
387 %! fail("schur ([1, 2; 3, 4], 2)","warning"); |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
388 |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
389 %!error <Invalid call to schur.*> schur (); |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
390 %!error schur ([1, 2, 3; 4, 5, 6]); |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
391 |
2928 | 392 */ |
10822 | 393 |
394 DEFUN_DLD (rsf2csf, args, nargout, | |
395 "-*- texinfo -*-\n\ | |
396 @deftypefn {Function File} {[@var{U}, @var{T}] =} rsf2csf (@var{UR}, @var{TR})\n\ | |
397 Converts a real, upper quasi-triangular Schur form @var{TR} to a complex,\n\ | |
398 upper triangular Schur form @var{T}.\n\ | |
399 \n\ | |
400 Note that the following relations hold: \n\ | |
401 \n\ | |
402 @code{@var{UR} * @var{TR} * @var{UR}' = @var{U} * @var{T} * @var{U}'} and\n\ | |
403 @code{@var{U}' * @var{U}} is identity.\n\ | |
404 \n\ | |
405 Note also that U and T are not unique.\n\ | |
406 \n\ | |
407 @end deftypefn") | |
408 { | |
409 octave_value_list retval; | |
410 | |
411 if (args.length () == 2 && nargout <= 2) | |
412 { | |
413 if (! args(0).is_numeric_type ()) | |
414 gripe_wrong_type_arg ("rsf2csf", args(0)); | |
415 else if (! args(1).is_numeric_type ()) | |
416 gripe_wrong_type_arg ("rsf2csf", args(1)); | |
417 else if (args(0).is_complex_type () || args(1).is_complex_type ()) | |
418 error ("rsf2csf: both matrices should be real"); | |
419 else | |
420 { | |
421 | |
422 if (args(0).is_single_type () || args(1).is_single_type ()) | |
423 { | |
424 FloatMatrix u = args(0).float_matrix_value (); | |
425 FloatMatrix t = args(1).float_matrix_value (); | |
426 if (! error_state) | |
427 { | |
428 FloatComplexSCHUR cs (FloatSCHUR (t, u)); | |
429 | |
430 retval(1) = cs.schur_matrix (); | |
431 retval(0) = cs.unitary_matrix (); | |
432 } | |
433 } | |
434 else | |
435 { | |
436 Matrix u = args(0).matrix_value (); | |
437 Matrix t = args(1).matrix_value (); | |
438 if (! error_state) | |
439 { | |
440 ComplexSCHUR cs (SCHUR (t, u)); | |
441 | |
442 retval(1) = cs.schur_matrix (); | |
443 retval(0) = cs.unitary_matrix (); | |
444 } | |
445 } | |
446 } | |
447 } | |
448 else | |
449 print_usage (); | |
450 | |
451 return retval; | |
452 } | |
453 | |
454 /* | |
455 | |
456 %!test | |
457 %! A = [1, 1, 1, 2; 1, 2, 1, 1; 1, 1, 3, 1; -2, 1, 1, 1]; | |
458 %! [u, t] = schur (A); | |
459 %! [U, T] = rsf2csf (u, t); | |
460 %! assert (norm (u * t * u' - U * T * U'), 0, 1e-12) | |
461 %! assert (norm (A - U * T * U'), 0, 1e-12) | |
462 | |
463 %!test | |
464 %! A = rand (10); | |
465 %! [u, t] = schur (A); | |
466 %! [U, T] = rsf2csf (u, t); | |
467 %! assert (norm (tril (T, -1)), 0) | |
468 %! assert (norm (U * U'), 1, 1e-14) | |
469 | |
470 */ |