← Back to branch summary

~ubuntu-branches/debian/sid/octave3.0/sid

« back to all changes in this revision

Viewing changes to src/DLD-FUNCTIONS/schur.cc

  • Committer: Bazaar Package Importer
  • Author(s): Rafael Laboissiere
  • Date: 2007-12-23 16:04:15 UTC
  • Revision ID: james.westby@ubuntu.com-20071223160415-n4gk468dihy22e9v
Tags: upstream-3.0.0
Import upstream version 3.0.0

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/DLD-FUNCTIONS/schur.cc
 
1
/*
 
2
 
 
3
Copyright (C) 1996, 1997, 1999, 2000, 2004, 2005, 2006, 2007
 
4
              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 3 of the License, or (at your
 
11
option) any 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, see
 
20
<http://www.gnu.org/licenses/>.
 
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"
 
32
 
 
33
#include "defun-dld.h"
 
34
#include "error.h"
 
35
#include "gripes.h"
 
36
#include "oct-obj.h"
 
37
#include "utils.h"
 
38
 
 
39
DEFUN_DLD (schur, args, nargout,
 
40
  "-*- texinfo -*-\n\
 
41
@deftypefn {Loadable Function} {@var{s} =} schur (@var{a})\n\
 
42
@deftypefnx {Loadable Function} {[@var{u}, @var{s}] =} schur (@var{a}, @var{opt})\n\
 
43
@cindex Schur decomposition\n\
 
44
The Schur decomposition is used to compute eigenvalues of a\n\
 
45
square matrix, and has applications in the solution of algebraic\n\
 
46
Riccati equations in control (see @code{are} and @code{dare}).\n\
 
47
@code{schur} always returns\n\
 
48
@iftex\n\
 
49
@tex\n\
 
50
$S = U^T A U$\n\
 
51
@end tex\n\
 
52
@end iftex\n\
 
53
@ifinfo\n\
 
54
@code{s = u' * a * u}\n\
 
55
@end ifinfo\n\
 
56
where\n\
 
57
@iftex\n\
 
58
@tex\n\
 
59
$U$\n\
 
60
@end tex\n\
 
61
@end iftex\n\
 
62
@ifinfo\n\
 
63
@code{u}\n\
 
64
@end ifinfo\n\
 
65
 is a unitary matrix\n\
 
66
@iftex\n\
 
67
@tex\n\
 
68
($U^T U$ is identity)\n\
 
69
@end tex\n\
 
70
@end iftex\n\
 
71
@ifinfo\n\
 
72
(@code{u'* u} is identity)\n\
 
73
@end ifinfo\n\
 
74
and\n\
 
75
@iftex\n\
 
76
@tex\n\
 
77
$S$\n\
 
78
@end tex\n\
 
79
@end iftex\n\
 
80
@ifinfo\n\
 
81
@code{s}\n\
 
82
@end ifinfo\n\
 
83
is upper triangular.  The eigenvalues of\n\
 
84
@iftex\n\
 
85
@tex\n\
 
86
$A$ (and $S$)\n\
 
87
@end tex\n\
 
88
@end iftex\n\
 
89
@ifinfo\n\
 
90
@code{a} (and @code{s})\n\
 
91
@end ifinfo\n\
 
92
are the diagonal elements of\n\
 
93
@iftex\n\
 
94
@tex\n\
 
95
$S$.\n\
 
96
@end tex\n\
 
97
@end iftex\n\
 
98
@ifinfo\n\
 
99
@code{s}.\n\
 
100
@end ifinfo\n\
 
101
If the matrix\n\
 
102
@iftex\n\
 
103
@tex\n\
 
104
$A$\n\
 
105
@end tex\n\
 
106
@end iftex\n\
 
107
@ifinfo\n\
 
108
@code{a}\n\
 
109
@end ifinfo\n\
 
110
is real, then the real Schur decomposition is computed, in which the\n\
 
111
matrix\n\
 
112
@iftex\n\
 
113
@tex\n\
 
114
$U$\n\
 
115
@end tex\n\
 
116
@end iftex\n\
 
117
@ifinfo\n\
 
118
@code{u}\n\
 
119
@end ifinfo\n\
 
120
is orthogonal and\n\
 
121
@iftex\n\
 
122
@tex\n\
 
123
$S$\n\
 
124
@end tex\n\
 
125
@end iftex\n\
 
126
@ifinfo\n\
 
127
@code{s}\n\
 
128
@end ifinfo\n\
 
129
is block upper triangular\n\
 
130
with blocks of size at most\n\
 
131
@iftex\n\
 
132
@tex\n\
 
133
$2\\times 2$\n\
 
134
@end tex\n\
 
135
@end iftex\n\
 
136
@ifinfo\n\
 
137
@code{2 x 2}\n\
 
138
@end ifinfo\n\
 
139
along the diagonal.  The diagonal elements of\n\
 
140
@iftex\n\
 
141
@tex\n\
 
142
$S$\n\
 
143
@end tex\n\
 
144
@end iftex\n\
 
145
@ifinfo\n\
 
146
@code{s}\n\
 
147
@end ifinfo\n\
 
148
(or the eigenvalues of the\n\
 
149
@iftex\n\
 
150
@tex\n\
 
151
$2\\times 2$\n\
 
152
@end tex\n\
 
153
@end iftex\n\
 
154
@ifinfo\n\
 
155
@code{2 x 2}\n\
 
156
@end ifinfo\n\
 
157
blocks, when\n\
 
158
appropriate) are the eigenvalues of\n\
 
159
@iftex\n\
 
160
@tex\n\
 
161
$A$\n\
 
162
@end tex\n\
 
163
@end iftex\n\
 
164
@ifinfo\n\
 
165
@code{a}\n\
 
166
@end ifinfo\n\
 
167
and\n\
 
168
@iftex\n\
 
169
@tex\n\
 
170
$S$.\n\
 
171
@end tex\n\
 
172
@end iftex\n\
 
173
@ifinfo\n\
 
174
@code{s}.\n\
 
175
@end ifinfo\n\
 
176
\n\
 
177
The eigenvalues are optionally ordered along the diagonal according to\n\
 
178
the value of @code{opt}.  @code{opt = \"a\"} indicates that all\n\
 
179
eigenvalues with negative real parts should be moved to the leading\n\
 
180
block of\n\
 
181
@iftex\n\
 
182
@tex\n\
 
183
$S$\n\
 
184
@end tex\n\
 
185
@end iftex\n\
 
186
@ifinfo\n\
 
187
@code{s}\n\
 
188
@end ifinfo\n\
 
189
(used in @code{are}), @code{opt = \"d\"} indicates that all eigenvalues\n\
 
190
with magnitude less than one should be moved to the leading block of\n\
 
191
@iftex\n\
 
192
@tex\n\
 
193
$S$\n\
 
194
@end tex\n\
 
195
@end iftex\n\
 
196
@ifinfo\n\
 
197
@code{s}\n\
 
198
@end ifinfo\n\
 
199
(used in @code{dare}), and @code{opt = \"u\"}, the default, indicates that\n\
 
200
no ordering of eigenvalues should occur.  The leading\n\
 
201
@iftex\n\
 
202
@tex\n\
 
203
$k$\n\
 
204
@end tex\n\
 
205
@end iftex\n\
 
206
@ifinfo\n\
 
207
@code{k}\n\
 
208
@end ifinfo\n\
 
209
columns of\n\
 
210
@iftex\n\
 
211
@tex\n\
 
212
$U$\n\
 
213
@end tex\n\
 
214
@end iftex\n\
 
215
@ifinfo\n\
 
216
@code{u}\n\
 
217
@end ifinfo\n\
 
218
always span the\n\
 
219
@iftex\n\
 
220
@tex\n\
 
221
$A$-invariant\n\
 
222
@end tex\n\
 
223
@end iftex\n\
 
224
@ifinfo\n\
 
225
@code{a}-invariant\n\
 
226
@end ifinfo\n\
 
227
subspace corresponding to the\n\
 
228
@iftex\n\
 
229
@tex\n\
 
230
$k$\n\
 
231
@end tex\n\
 
232
@end iftex\n\
 
233
@ifinfo\n\
 
234
@code{k}\n\
 
235
@end ifinfo\n\
 
236
leading eigenvalues of\n\
 
237
@iftex\n\
 
238
@tex\n\
 
239
$S$.\n\
 
240
@end tex\n\
 
241
@end iftex\n\
 
242
@ifinfo\n\
 
243
@code{s}.\n\
 
244
@end ifinfo\n\
 
245
@end deftypefn")
 
246
{
 
247
  octave_value_list retval;
 
248
 
 
249
  int nargin = args.length ();
 
250
 
 
251
  if (nargin < 1 || nargin > 2 || nargout > 2)
 
252
    {
 
253
      print_usage ();
 
254
      return retval;
 
255
    }
 
256
 
 
257
  octave_value arg = args(0);
 
258
 
 
259
  std::string ord;
 
260
 
 
261
  if (nargin == 2)
 
262
    {
 
263
      ord = args(1).string_value (); 
 
264
 
 
265
      if (error_state)
 
266
        {
 
267
          error ("schur: expecting string as second argument");
 
268
          return retval;
 
269
        }
 
270
    }
 
271
 
 
272
  char ord_char = ord.empty () ? 'U' : ord[0];
 
273
 
 
274
  if (ord_char != 'U' && ord_char != 'A' && ord_char != 'D'
 
275
      && ord_char != 'u' && ord_char != 'a' && ord_char != 'd')
 
276
    {
 
277
      warning ("schur: incorrect ordered schur argument `%c'",
 
278
               ord.c_str ());
 
279
      return retval;
 
280
    }
 
281
 
 
282
  octave_idx_type nr = arg.rows ();
 
283
  octave_idx_type nc = arg.columns ();
 
284
 
 
285
  int arg_is_empty = empty_arg ("schur", nr, nc);
 
286
 
 
287
  if (arg_is_empty < 0)
 
288
    return retval;
 
289
  else if (arg_is_empty > 0)
 
290
    return octave_value_list (2, Matrix ());
 
291
 
 
292
  if (nr != nc)
 
293
    {
 
294
      gripe_square_matrix_required ("schur");
 
295
      return retval;
 
296
    }
 
297
 
 
298
  if (arg.is_real_type ())
 
299
    {
 
300
      Matrix tmp = arg.matrix_value ();
 
301
 
 
302
      if (! error_state)
 
303
        {
 
304
          if (nargout == 0 || nargout == 1)
 
305
            {
 
306
              SCHUR result (tmp, ord, false);
 
307
              retval(0) = result.schur_matrix ();
 
308
            }
 
309
          else
 
310
            {
 
311
              SCHUR result (tmp, ord, true);
 
312
              retval(1) = result.schur_matrix ();
 
313
              retval(0) = result.unitary_matrix ();
 
314
            }
 
315
        }
 
316
    }
 
317
  else if (arg.is_complex_type ())
 
318
    {
 
319
      ComplexMatrix ctmp = arg.complex_matrix_value ();
 
320
 
 
321
      if (! error_state)
 
322
        {
 
323
 
 
324
          if (nargout == 0 || nargout == 1)
 
325
            {
 
326
              ComplexSCHUR result (ctmp, ord, false);
 
327
              retval(0) = result.schur_matrix ();
 
328
            }
 
329
          else
 
330
            {
 
331
              ComplexSCHUR result (ctmp, ord, true);
 
332
              retval(1) = result.schur_matrix ();
 
333
              retval(0) = result.unitary_matrix ();
 
334
            }
 
335
        }
 
336
    }    
 
337
  else
 
338
    {
 
339
      gripe_wrong_type_arg ("schur", arg);
 
340
    }
 
341
 
 
342
  return retval; 
 
343
}
 
344
 
 
345
/*
 
346
;;; Local Variables: ***
 
347
;;; mode: C++ ***
 
348
;;; End: ***
 
349
*/