GNU Octave  4.0.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
ordschur.cc
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 2015 S├ębastien Villemot <sebastien@debian.org>
4 
5 This file is part of Octave.
6 
7 Octave is free software; you can redistribute it and/or modify it
8 under the terms of the GNU General Public License as published by the
9 Free Software Foundation; either version 3 of the License, or (at your
10 option) any later version.
11 
12 Octave is distributed in the hope that it will be useful, but WITHOUT
13 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
14 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
15 for more details.
16 
17 You should have received a copy of the GNU General Public License
18 along with Octave; see the file COPYING. If not, see
19 <http://www.gnu.org/licenses/>.
20 
21 */
22 
23 #ifdef HAVE_CONFIG_H
24 #include <config.h>
25 #endif
26 
27 #include "defun.h"
28 #include "error.h"
29 #include "oct-obj.h"
30 #include "f77-fcn.h"
31 
32 extern "C"
33 {
34  F77_RET_T
36  const octave_idx_type*, const octave_idx_type&,
37  double*, const octave_idx_type&, double*, const octave_idx_type&,
38  double*, double*, octave_idx_type&, double&, double&, double*,
39  const octave_idx_type&, octave_idx_type*,
40  const octave_idx_type&, octave_idx_type&);
41 
42  F77_RET_T
44  const octave_idx_type*, const octave_idx_type&,
45  Complex*, const octave_idx_type&, Complex*, const octave_idx_type&,
46  Complex*, octave_idx_type&, double&, double&, Complex*,
47  const octave_idx_type&, octave_idx_type &);
48 
49  F77_RET_T
51  const octave_idx_type*, const octave_idx_type&,
52  float*, const octave_idx_type&, float*, const octave_idx_type&,
53  float*, float*, octave_idx_type&, float&, float&, float*,
54  const octave_idx_type&, octave_idx_type*,
55  const octave_idx_type&, octave_idx_type&);
56 
57  F77_RET_T
59  const octave_idx_type*, const octave_idx_type&,
60  FloatComplex*, const octave_idx_type&, FloatComplex*, const octave_idx_type&,
61  FloatComplex*, octave_idx_type&, float&, float&, FloatComplex*,
62  const octave_idx_type&, octave_idx_type &);
63 }
64 
65 DEFUN (ordschur, args, ,
66  "-*- texinfo -*-\n\
67 @deftypefn {Loadable Function} {[@var{UR}, @var{SR}] =} ordschur (@var{U}, @var{S}, @var{select})\n\
68 Reorders the real Schur factorization (@var{U},@var{S}) obtained with the\n\
69 @code{schur} function, so that selected eigenvalues appear in the upper left\n\
70 diagonal blocks of the quasi triangular Schur matrix.\n\
71 \n\
72 The logical vector @var{select} specifies the selected eigenvalues as they\n\
73 appear along @var{S}'s diagonal.\n\
74 \n\
75 For example, given the matrix @code{@var{A} = [1, 2; 3, 4]}, and its Schur\n\
76 decomposition\n\
77 \n\
78 @example\n\
79 [@var{U}, @var{S}] = schur (@var{A})\n\
80 @end example\n\
81 \n\
82 @noindent\n\
83 which returns\n\
84 \n\
85 @example\n\
86 @group\n\
87 @var{U} =\n\
88 \n\
89  -0.82456 -0.56577\n\
90  0.56577 -0.82456\n\
91 \n\
92 @var{S} =\n\
93 \n\
94  -0.37228 -1.00000\n\
95  0.00000 5.37228\n\
96 \n\
97 @end group\n\
98 @end example\n\
99 \n\
100 It is possible to reorder the decomposition so that the positive eigenvalue\n\
101 is in the upper left corner, by doing:\n\
102 \n\
103 @example\n\
104 [@var{U}, @var{S}] = ordschur (@var{U}, @var{S}, [0,1])\n\
105 @end example\n\
106 \n\
107 @seealso{schur}\n\
108 @end deftypefn")
109 {
110  const octave_idx_type nargin = args.length ();
111  octave_value_list retval;
112 
113  if (nargin != 3)
114  {
115  print_usage ();
116  return retval;
117  }
118 
119  const Array<octave_idx_type> sel = args(2).octave_idx_type_vector_value ();
120  if (error_state)
121  {
122  error ("ordschur: SELECT must be an array of integers");
123  return retval;
124  }
125  const octave_idx_type n = sel.numel ();
126 
127  const dim_vector dimU = args(0).dims ();
128  const dim_vector dimS = args(1).dims ();
129  if (n != dimU(0))
130  {
131  error ("ordschur: SELECT must have same length as the sides of U and S");
132  return retval;
133  }
134  else if (n != dimU(0) || n != dimS(0) || n != dimU(1) || n != dimS(1))
135  {
136  error ("ordschur: U and S must be square and of equal sizes");
137  return retval;
138  }
139 
140  const bool double_type = args(0).is_double_type ()
141  || args(1).is_double_type ();
142  const bool complex_type = args(0).is_complex_type ()
143  || args(1).is_complex_type ();
144 
145 #define PREPARE_ARGS(TYPE, TYPE_M, TYPE_COND) \
146  TYPE ## Matrix U = args(0).TYPE_M ## _value (); \
147  TYPE ## Matrix S = args(1).TYPE_M ## _value (); \
148  if (error_state) \
149  { \
150  error ("ordschur: U and S must be real or complex floating point matrices"); \
151  return retval; \
152  } \
153  TYPE ## Matrix w (dim_vector (n, 1)); \
154  TYPE ## Matrix work (dim_vector (n, 1)); \
155  octave_idx_type m; \
156  octave_idx_type info; \
157  TYPE_COND cond1, cond2;
158 
159 #define PREPARE_OUTPUT()\
160  if (info != 0) \
161  { \
162  error ("ordschur: trsen failed"); \
163  return retval; \
164  } \
165  retval(0) = U; \
166  retval(1) = S;
167 
168  if (double_type)
169  {
170  if (complex_type)
171  {
172  PREPARE_ARGS (Complex, complex_matrix, double)
173 
174  F77_XFCN (ztrsen, ztrsen,
176  sel.data (), n, S.fortran_vec (), n, U.fortran_vec (), n,
177  w.fortran_vec (), m, cond1, cond2, work.fortran_vec (), n,
178  info));
180  }
181  else
182  {
183  PREPARE_ARGS (, matrix, double)
184  Matrix wi (dim_vector (n, 1));
185  Array<octave_idx_type> iwork (dim_vector (n, 1));
186 
187  F77_XFCN (dtrsen, dtrsen,
189  sel.data (), n, S.fortran_vec (), n, U.fortran_vec (), n,
190  w.fortran_vec (), wi.fortran_vec (), m, cond1, cond2,
191  work.fortran_vec (), n, iwork.fortran_vec (), n, info));
192  PREPARE_OUTPUT ()
193  }
194  }
195  else
196  {
197  if (complex_type)
198  {
199  PREPARE_ARGS (FloatComplex, float_complex_matrix, float)
200 
201  F77_XFCN (ctrsen, ctrsen,
203  sel.data (), n, S.fortran_vec (), n, U.fortran_vec (), n,
204  w.fortran_vec (), m, cond1, cond2, work.fortran_vec (), n,
205  info));
206  PREPARE_OUTPUT ()
207  }
208  else
209  {
210  PREPARE_ARGS (Float, float_matrix, float)
211  FloatMatrix wi (dim_vector (n, 1));
212  Array<octave_idx_type> iwork (dim_vector (n, 1));
213 
214  F77_XFCN (strsen, strsen,
216  sel.data (), n, S.fortran_vec (), n, U.fortran_vec (), n,
217  w.fortran_vec (), wi.fortran_vec (), m, cond1, cond2,
218  work.fortran_vec (), n, iwork.fortran_vec (), n, info));
219  PREPARE_OUTPUT ()
220  }
221  }
222 
223 #undef PREPARE_ARGS
224 #undef PREPARE_OUTPUT
225 
226  return retval;
227 }
228 
229 /*
230 
231 %!test
232 %! A = [1, 2, 3, -2; 4, 5, 6, -5 ; 7, 8, 9, -5; 10, 11, 12, 4 ];
233 %! [U, T] = schur (A);
234 %! [US, TS] = ordschur (U, T, [ 0, 0, 1, 1 ]);
235 %! assert (US*TS*US', A, sqrt (eps))
236 %! assert (diag (T)(3:4), diag (TS)(1:2), sqrt (eps))
237 
238 %!test
239 %! A = [1, 2, 3, -2; 4, 5, 6, -5 ; 7, 8, 9, -5; 10, 11, 12, 4 ];
240 %! [U, T] = schur (A);
241 %! [US, TS] = ordschur (single (U), single (T), [ 0, 0, 1, 1 ]);
242 %! assert (US*TS*US', A, sqrt (eps ("single")))
243 %! assert (diag (T)(3:4), diag (TS)(1:2), sqrt (eps ("single")))
244 
245 %!test
246 %! A = [1, 2, 3, -2; 4, 5, 6, -5 ; 7, 8, 9, -5; 10, 11, 12, 4+3i ];
247 %! [U, T] = schur (A);
248 %! [US, TS] = ordschur (U, T, [ 0, 0, 1, 1 ]);
249 %! assert (US*TS*US', A, sqrt (eps))
250 %! assert (diag (T)(3:4), diag (TS)(1:2), sqrt (eps))
251 
252 %!test
253 %! A = [1, 2, 3, -2; 4, 5, 6, -5 ; 7, 8, 9, -5; 10, 11, 12, 4+3i ];
254 %! [U, T] = schur (A);
255 %! [US, TS] = ordschur (single (U), single (T), [ 0, 0, 1, 1 ]);
256 %! assert (US*TS*US', A, sqrt (eps ("single")))
257 %! assert (diag (T)(3:4), diag (TS)(1:2), sqrt (eps ("single")))
258 
259 */
F77_RET_T F77_CONST_CHAR_ARG_DECL
Definition: ordschur.cc:35
#define PREPARE_ARGS(TYPE, TYPE_M, TYPE_COND)
OCTINTERP_API void print_usage(void)
Definition: defun.cc:51
octave_idx_type numel(void) const
Number of elements in the array.
Definition: Array.h:275
#define DEFUN(name, args_name, nargout_name, doc)
Definition: defun.h:44
void error(const char *fmt,...)
Definition: error.cc:476
#define F77_XFCN(f, F, args)
Definition: f77-fcn.h:51
#define PREPARE_OUTPUT()
F77_RET_T F77_FUNC(dtrsen, DTRSEN)(F77_CONST_CHAR_ARG_DECL
std::complex< double > w(std::complex< double > z, double relerr=0)
const T * data(void) const
Definition: Array.h:479
int error_state
Definition: error.cc:101
#define F77_RET_T
Definition: f77-fcn.h:264
Definition: dMatrix.h:35
#define F77_CONST_CHAR_ARG(x)
Definition: f77-fcn.h:249
static double wi[256]
Definition: randmtzig.c:443
std::complex< float > FloatComplex
Definition: oct-cmplx.h:30
std::complex< double > Complex
Definition: oct-cmplx.h:29
const T * fortran_vec(void) const
Definition: Array.h:481