← Back to branch summary

~gabriel1984sibiu/octave/octave

« back to all changes in this revision

Viewing changes to liboctave/numeric/floatSCHUR.cc

  • Committer: Grevutiu Gabriel
  • Date: 2014-01-02 13:05:54 UTC
  • Revision ID: gabriel1984sibiu@gmail.com-20140102130554-3r7ivdjln1ni6kcg
New version (3.8.0) from upstream.

Show diffs side-by-side

added added

removed removed

Lines of Context:
liboctave/numeric/floatSCHUR.cc
 
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,
 
37
                             F77_CONST_CHAR_ARG_DECL,
 
38
                             FloatSCHUR::select_function,
 
39
                             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&
 
46
                             F77_CHAR_ARG_LEN_DECL
 
47
                             F77_CHAR_ARG_LEN_DECL
 
48
                             F77_CHAR_ARG_LEN_DECL);
 
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
 
 
63
octave_idx_type
 
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')
 
100
    selector = select_ana;
 
101
  else if (ord_char == 'D' || ord_char == 'd')
 
102
    selector = select_dig;
 
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
 
 
152
FloatSCHUR::FloatSCHUR (const FloatMatrix& s, const FloatMatrix& u)
 
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)
 
157
    (*current_liboctave_error_handler)
 
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
}