GNU Octave  4.2.1
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) 2016-2017 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 #if defined (HAVE_CONFIG_H)
24 # include "config.h"
25 #endif
26 
27 #include "defun.h"
28 #include "error.h"
29 #include "lo-lapack-proto.h"
30 #include "ovl.h"
31 
32 DEFUN (ordschur, args, ,
33  doc: /* -*- texinfo -*-
34 @deftypefn {} {[@var{UR}, @var{SR}] =} ordschur (@var{U}, @var{S}, @var{select})
35 Reorders the real Schur factorization (@var{U},@var{S}) obtained with the
36 @code{schur} function, so that selected eigenvalues appear in the upper left
37 diagonal blocks of the quasi triangular Schur matrix.
38 
39 The logical vector @var{select} specifies the selected eigenvalues as they
40 appear along @var{S}'s diagonal.
41 
42 For example, given the matrix @code{@var{A} = [1, 2; 3, 4]}, and its Schur
43 decomposition
44 
45 @example
46 [@var{U}, @var{S}] = schur (@var{A})
47 @end example
48 
49 @noindent
50 which returns
51 
52 @example
53 @group
54 @var{U} =
55 
56  -0.82456 -0.56577
57  0.56577 -0.82456
58 
59 @var{S} =
60 
61  -0.37228 -1.00000
62  0.00000 5.37228
63 
64 @end group
65 @end example
66 
67 It is possible to reorder the decomposition so that the positive eigenvalue
68 is in the upper left corner, by doing:
69 
70 @example
71 [@var{U}, @var{S}] = ordschur (@var{U}, @var{S}, [0,1])
72 @end example
73 
74 @seealso{schur}
75 @end deftypefn */)
76 {
77  if (args.length () != 3)
78  print_usage ();
79 
80  const Array<octave_idx_type> sel = args(2).octave_idx_type_vector_value ("ordschur: SELECT must be an array of integers");
81 
82  const octave_idx_type n = sel.numel ();
83 
84  const dim_vector dimU = args(0).dims ();
85  const dim_vector dimS = args(1).dims ();
86 
87  if (n != dimU(0))
88  error ("ordschur: SELECT must have same length as the sides of U and S");
89  else if (n != dimU(0) || n != dimS(0) || n != dimU(1) || n != dimS(1))
90  error ("ordschur: U and S must be square and of equal sizes");
91 
93 
94  const bool double_type = args(0).is_double_type ()
95  || args(1).is_double_type ();
96  const bool complex_type = args(0).is_complex_type ()
97  || args(1).is_complex_type ();
98 
99 #define PREPARE_ARGS(TYPE, TYPE_M, TYPE_COND) \
100  TYPE ## Matrix U = args(0).x ## TYPE_M ## _value ("ordschur: U and S must be real or complex floating point matrices"); \
101  TYPE ## Matrix S = args(1).x ## TYPE_M ## _value ("ordschur: U and S must be real or complex floating point matrices"); \
102  TYPE ## Matrix w (dim_vector (n, 1)); \
103  TYPE ## Matrix work (dim_vector (n, 1)); \
104  octave_idx_type m; \
105  octave_idx_type info; \
106  TYPE_COND cond1, cond2;
107 
108 #define PREPARE_OUTPUT() \
109  if (info != 0) \
110  error ("ordschur: trsen failed"); \
111  \
112  retval = ovl (U, S);
113 
114  if (double_type)
115  {
116  if (complex_type)
117  {
118  PREPARE_ARGS (Complex, complex_matrix, double)
119 
120  F77_XFCN (ztrsen, ztrsen,
121  (F77_CONST_CHAR_ARG ("N"), F77_CONST_CHAR_ARG ("V"),
122  sel.data (), n, F77_DBLE_CMPLX_ARG (S.fortran_vec ()), n,
123  F77_DBLE_CMPLX_ARG (U.fortran_vec ()), n,
124  F77_DBLE_CMPLX_ARG (w.fortran_vec ()), m, cond1, cond2,
125  F77_DBLE_CMPLX_ARG (work.fortran_vec ()), n,
126  info));
127 
129  }
130  else
131  {
132  PREPARE_ARGS (, matrix, double)
133  Matrix wi (dim_vector (n, 1));
134  Array<octave_idx_type> iwork (dim_vector (n, 1));
135 
136  F77_XFCN (dtrsen, dtrsen,
137  (F77_CONST_CHAR_ARG ("N"), F77_CONST_CHAR_ARG ("V"),
138  sel.data (), n, S.fortran_vec (), n, U.fortran_vec (), n,
139  w.fortran_vec (), wi.fortran_vec (), m, cond1, cond2,
140  work.fortran_vec (), n, iwork.fortran_vec (), n, info));
141 
142  PREPARE_OUTPUT ()
143  }
144  }
145  else
146  {
147  if (complex_type)
148  {
149  PREPARE_ARGS (FloatComplex, float_complex_matrix, float)
150 
151  F77_XFCN (ctrsen, ctrsen,
152  (F77_CONST_CHAR_ARG ("N"), F77_CONST_CHAR_ARG ("V"),
153  sel.data (), n, F77_CMPLX_ARG (S.fortran_vec ()), n,
154  F77_CMPLX_ARG (U.fortran_vec ()), n,
155  F77_CMPLX_ARG (w.fortran_vec ()), m, cond1, cond2,
156  F77_CMPLX_ARG (work.fortran_vec ()), n,
157  info));
158 
159  PREPARE_OUTPUT ()
160  }
161  else
162  {
163  PREPARE_ARGS (Float, float_matrix, float)
164  FloatMatrix wi (dim_vector (n, 1));
165  Array<octave_idx_type> iwork (dim_vector (n, 1));
166 
167  F77_XFCN (strsen, strsen,
168  (F77_CONST_CHAR_ARG ("N"), F77_CONST_CHAR_ARG ("V"),
169  sel.data (), n, S.fortran_vec (), n, U.fortran_vec (), n,
170  w.fortran_vec (), wi.fortran_vec (), m, cond1, cond2,
171  work.fortran_vec (), n, iwork.fortran_vec (), n, info));
172 
173  PREPARE_OUTPUT ()
174  }
175  }
176 
177 #undef PREPARE_ARGS
178 #undef PREPARE_OUTPUT
179 
180  return retval;
181 }
182 
183 /*
184 
185 %!test
186 %! A = [1, 2, 3, -2; 4, 5, 6, -5 ; 7, 8, 9, -5; 10, 11, 12, 4 ];
187 %! [U, T] = schur (A);
188 %! [US, TS] = ordschur (U, T, [ 0, 0, 1, 1 ]);
189 %! assert (US*TS*US', A, sqrt (eps));
190 %! assert (diag (T)(3:4), diag (TS)(1:2), sqrt (eps));
191 
192 %!test
193 %! A = [1, 2, 3, -2; 4, 5, 6, -5 ; 7, 8, 9, -5; 10, 11, 12, 4 ];
194 %! [U, T] = schur (A);
195 %! [US, TS] = ordschur (single (U), single (T), [ 0, 0, 1, 1 ]);
196 %! assert (US*TS*US', A, sqrt (eps ("single")));
197 %! assert (diag (T)(3:4), diag (TS)(1:2), sqrt (eps ("single")));
198 
199 %!test
200 %! A = [1, 2, 3, -2; 4, 5, 6, -5 ; 7, 8, 9, -5; 10, 11, 12, 4+3i ];
201 %! [U, T] = schur (A);
202 %! [US, TS] = ordschur (U, T, [ 0, 0, 1, 1 ]);
203 %! assert (US*TS*US', A, sqrt (eps));
204 %! assert (diag (T)(3:4), diag (TS)(1:2), sqrt (eps));
205 
206 %!test
207 %! A = [1, 2, 3, -2; 4, 5, 6, -5 ; 7, 8, 9, -5; 10, 11, 12, 4+3i ];
208 %! [U, T] = schur (A);
209 %! [US, TS] = ordschur (single (U), single (T), [ 0, 0, 1, 1 ]);
210 %! assert (US*TS*US', A, sqrt (eps ("single")));
211 %! assert (diag (T)(3:4), diag (TS)(1:2), sqrt (eps ("single")));
212 
213 */
#define PREPARE_ARGS(TYPE, TYPE_M, TYPE_COND)
OCTINTERP_API void print_usage(void)
Definition: defun.cc:52
octave_idx_type numel(void) const
Number of elements in the array.
Definition: Array.h:363
#define F77_DBLE_CMPLX_ARG(x)
Definition: f77-fcn.h:345
the sparsity preserving column transformation such that that defines the pivoting threshold can be given in which case it defines the then the first element defines the pivoting tolerance for the unsymmetric the values defined such that for full matrix
Definition: lu.cc:138
#define DEFUN(name, args_name, nargout_name, doc)
Definition: defun.h:46
void error(const char *fmt,...)
Definition: error.cc:570
#define F77_XFCN(f, F, args)
Definition: f77-fcn.h:52
#define PREPARE_OUTPUT()
JNIEnv void * args
Definition: ov-java.cc:67
nd deftypefn *octave_map m
Definition: ov-struct.cc:2058
std::complex< double > w(std::complex< double > z, double relerr=0)
const T * data(void) const
Definition: Array.h:582
octave_value retval
Definition: data.cc:6294
Definition: dMatrix.h:37
#define F77_CMPLX_ARG(x)
Definition: f77-fcn.h:339
static double wi[256]
Definition: randmtzig.cc:440
std::complex< float > FloatComplex
Definition: oct-cmplx.h:32
std::complex< double > Complex
Definition: oct-cmplx.h:31
const T * fortran_vec(void) const
Definition: Array.h:584
Vector representing the dimensions (size) of an Array.
Definition: dim-vector.h:87