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
floatSVD.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 "floatSVD.h"
30 #include "f77-fcn.h"
31 #include "oct-locbuf.h"
32 
33 extern "C"
34 {
35  F77_RET_T
36  F77_FUNC (sgesvd, SGESVD) (F77_CONST_CHAR_ARG_DECL,
38  const octave_idx_type&, const octave_idx_type&,
39  float*, const octave_idx_type&, float*,
40  float*, const octave_idx_type&, float*,
41  const octave_idx_type&, float*,
42  const octave_idx_type&, octave_idx_type&
45 
46  F77_RET_T
47  F77_FUNC (sgesdd, SGESDD) (F77_CONST_CHAR_ARG_DECL,
48  const octave_idx_type&, const octave_idx_type&,
49  float*, const octave_idx_type&, float*,
50  float*, const octave_idx_type&, float*,
51  const octave_idx_type&, float*,
52  const octave_idx_type&, octave_idx_type *,
53  octave_idx_type&
55 }
56 
59 {
61  {
62  (*current_liboctave_error_handler)
63  ("FloatSVD: U not computed because type == SVD::sigma_only");
64  return FloatMatrix ();
65  }
66  else
67  return left_sm;
68 }
69 
72 {
74  {
75  (*current_liboctave_error_handler)
76  ("FloatSVD: V not computed because type == SVD::sigma_only");
77  return FloatMatrix ();
78  }
79  else
80  return right_sm;
81 }
82 
84 FloatSVD::init (const FloatMatrix& a, SVD::type svd_type,
85  SVD::driver svd_driver)
86 {
87  octave_idx_type info;
88 
89  octave_idx_type m = a.rows ();
90  octave_idx_type n = a.cols ();
91 
92  FloatMatrix atmp = a;
93  float *tmp_data = atmp.fortran_vec ();
94 
95  octave_idx_type min_mn = m < n ? m : n;
96 
97  char jobu = 'A';
98  char jobv = 'A';
99 
100  octave_idx_type ncol_u = m;
101  octave_idx_type nrow_vt = n;
102  octave_idx_type nrow_s = m;
103  octave_idx_type ncol_s = n;
104 
105  switch (svd_type)
106  {
107  case SVD::economy:
108  jobu = jobv = 'S';
109  ncol_u = nrow_vt = nrow_s = ncol_s = min_mn;
110  break;
111 
112  case SVD::sigma_only:
113 
114  // Note: for this case, both jobu and jobv should be 'N', but
115  // there seems to be a bug in dgesvd from Lapack V2.0. To
116  // demonstrate the bug, set both jobu and jobv to 'N' and find
117  // the singular values of [eye(3), eye(3)]. The result is
118  // [-sqrt(2), -sqrt(2), -sqrt(2)].
119  //
120  // For Lapack 3.0, this problem seems to be fixed.
121 
122  jobu = jobv = 'N';
123  ncol_u = nrow_vt = 1;
124  break;
125 
126  default:
127  break;
128  }
129 
130  type_computed = svd_type;
131 
132  if (! (jobu == 'N' || jobu == 'O'))
133  left_sm.resize (m, ncol_u);
134 
135  float *u = left_sm.fortran_vec ();
136 
137  sigma.resize (nrow_s, ncol_s);
138  float *s_vec = sigma.fortran_vec ();
139 
140  if (! (jobv == 'N' || jobv == 'O'))
141  right_sm.resize (nrow_vt, n);
142 
143  float *vt = right_sm.fortran_vec ();
144 
145  // Query SGESVD for the correct dimension of WORK.
146 
147  octave_idx_type lwork = -1;
148 
149  Array<float> work (dim_vector (1, 1));
150 
151  octave_idx_type one = 1;
152  octave_idx_type m1 = std::max (m, one);
153  octave_idx_type nrow_vt1 = std::max (nrow_vt, one);
154 
155  if (svd_driver == SVD::GESVD)
156  {
157  F77_XFCN (sgesvd, SGESVD, (F77_CONST_CHAR_ARG2 (&jobu, 1),
158  F77_CONST_CHAR_ARG2 (&jobv, 1),
159  m, n, tmp_data, m1, s_vec, u, m1, vt,
160  nrow_vt1, work.fortran_vec (), lwork, info
161  F77_CHAR_ARG_LEN (1)
162  F77_CHAR_ARG_LEN (1)));
163 
164  lwork = static_cast<octave_idx_type> (work(0));
165  work.resize (dim_vector (lwork, 1));
166 
167  F77_XFCN (sgesvd, SGESVD, (F77_CONST_CHAR_ARG2 (&jobu, 1),
168  F77_CONST_CHAR_ARG2 (&jobv, 1),
169  m, n, tmp_data, m1, s_vec, u, m1, vt,
170  nrow_vt1, work.fortran_vec (), lwork, info
171  F77_CHAR_ARG_LEN (1)
172  F77_CHAR_ARG_LEN (1)));
173 
174  }
175  else if (svd_driver == SVD::GESDD)
176  {
177  assert (jobu == jobv);
178  char jobz = jobu;
179  OCTAVE_LOCAL_BUFFER (octave_idx_type, iwork, 8*min_mn);
180 
181  F77_XFCN (sgesdd, SGESDD, (F77_CONST_CHAR_ARG2 (&jobz, 1),
182  m, n, tmp_data, m1, s_vec, u, m1, vt, nrow_vt1,
183  work.fortran_vec (), lwork, iwork, info
184  F77_CHAR_ARG_LEN (1)));
185 
186  lwork = static_cast<octave_idx_type> (work(0));
187  work.resize (dim_vector (lwork, 1));
188 
189  F77_XFCN (sgesdd, SGESDD, (F77_CONST_CHAR_ARG2 (&jobz, 1),
190  m, n, tmp_data, m1, s_vec, u, m1, vt, nrow_vt1,
191  work.fortran_vec (), lwork, iwork, info
192  F77_CHAR_ARG_LEN (1)));
193 
194  }
195  else
196  assert (0); // impossible
197 
198  if (! (jobv == 'N' || jobv == 'O'))
200 
201  return info;
202 }
203 
204 std::ostream&
205 operator << (std::ostream& os, const FloatSVD& a)
206 {
207  os << a.left_singular_matrix () << "\n";
208  os << a.singular_values () << "\n";
209  os << a.right_singular_matrix () << "\n";
210 
211  return os;
212 }