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