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