GNU Octave  3.8.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
floatSCHUR.cc
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 1994-2013 John W. Eaton
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 <iostream>
28 
29 #include "floatSCHUR.h"
30 #include "f77-fcn.h"
31 #include "lo-error.h"
32 
33 extern "C"
34 {
35  F77_RET_T
36  F77_FUNC (sgeesx, SGEESX) (F77_CONST_CHAR_ARG_DECL,
40  const octave_idx_type&, float*,
41  const octave_idx_type&, octave_idx_type&,
42  float*, float*, float*, const octave_idx_type&,
43  float&, float&, float*, const octave_idx_type&,
44  octave_idx_type*, const octave_idx_type&,
45  octave_idx_type*, octave_idx_type&
49 }
50 
51 static octave_idx_type
52 select_ana (const float& a, const float&)
53 {
54  return (a < 0.0);
55 }
56 
57 static octave_idx_type
58 select_dig (const float& a, const float& b)
59 {
60  return (hypot (a, b) < 1.0);
61 }
62 
64 FloatSCHUR::init (const FloatMatrix& a, const std::string& ord,
65  bool calc_unitary)
66 {
67  octave_idx_type a_nr = a.rows ();
68  octave_idx_type a_nc = a.cols ();
69 
70  if (a_nr != a_nc)
71  {
72  (*current_liboctave_error_handler) ("FloatSCHUR requires square matrix");
73  return -1;
74  }
75  else if (a_nr == 0)
76  {
77  schur_mat.clear ();
78  unitary_mat.clear ();
79  return 0;
80  }
81 
82  // Workspace requirements may need to be fixed if any of the
83  // following change.
84 
85  char jobvs;
86  char sense = 'N';
87  char sort = 'N';
88 
89  if (calc_unitary)
90  jobvs = 'V';
91  else
92  jobvs = 'N';
93 
94  char ord_char = ord.empty () ? 'U' : ord[0];
95 
96  if (ord_char == 'A' || ord_char == 'D' || ord_char == 'a' || ord_char == 'd')
97  sort = 'S';
98 
99  if (ord_char == 'A' || ord_char == 'a')
101  else if (ord_char == 'D' || ord_char == 'd')
103  else
104  selector = 0;
105 
106  octave_idx_type n = a_nc;
107  octave_idx_type lwork = 8 * n;
108  octave_idx_type liwork = 1;
109  octave_idx_type info;
110  octave_idx_type sdim;
111  float rconde;
112  float rcondv;
113 
114  schur_mat = a;
115 
116  if (calc_unitary)
117  unitary_mat.clear (n, n);
118 
119  float *s = schur_mat.fortran_vec ();
120  float *q = unitary_mat.fortran_vec ();
121 
122  Array<float> wr (dim_vector (n, 1));
123  float *pwr = wr.fortran_vec ();
124 
125  Array<float> wi (dim_vector (n, 1));
126  float *pwi = wi.fortran_vec ();
127 
128  Array<float> work (dim_vector (lwork, 1));
129  float *pwork = work.fortran_vec ();
130 
131  // BWORK is not referenced for the non-ordered Schur routine.
132  octave_idx_type ntmp = (ord_char == 'N' || ord_char == 'n') ? 0 : n;
133  Array<octave_idx_type> bwork (dim_vector (ntmp, 1));
134  octave_idx_type *pbwork = bwork.fortran_vec ();
135 
136  Array<octave_idx_type> iwork (dim_vector (liwork, 1));
137  octave_idx_type *piwork = iwork.fortran_vec ();
138 
139  F77_XFCN (sgeesx, SGEESX, (F77_CONST_CHAR_ARG2 (&jobvs, 1),
140  F77_CONST_CHAR_ARG2 (&sort, 1),
141  selector,
142  F77_CONST_CHAR_ARG2 (&sense, 1),
143  n, s, n, sdim, pwr, pwi, q, n, rconde, rcondv,
144  pwork, lwork, piwork, liwork, pbwork, info
145  F77_CHAR_ARG_LEN (1)
146  F77_CHAR_ARG_LEN (1)
147  F77_CHAR_ARG_LEN (1)));
148 
149  return info;
150 }
151 
153  : schur_mat (s), unitary_mat (u), selector (0)
154 {
155  octave_idx_type n = s.rows ();
156  if (s.columns () != n || u.rows () != n || u.columns () != n)
158  ("schur: inconsistent matrix dimensions");
159 }
160 
161 std::ostream&
162 operator << (std::ostream& os, const FloatSCHUR& a)
163 {
164  os << a.schur_matrix () << "\n";
165  os << a.unitary_matrix () << "\n";
166 
167  return os;
168 }