GNU Octave  4.4.1
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
svd.cc
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 2016-2018 CarnĂ« Draug
4 Copyright (C) 1994-2018 John W. Eaton
5 
6 This file is part of Octave.
7 
8 Octave is free software: you can redistribute it and/or modify it
9 under the terms of the GNU General Public License as published by
10 the Free Software Foundation, either version 3 of the License, or
11 (at your option) any later version.
12 
13 Octave is distributed in the hope that it will be useful, but
14 WITHOUT ANY WARRANTY; without even the implied warranty of
15 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 GNU General Public License for more details.
17 
18 You should have received a copy of the GNU General Public License
19 along with Octave; see the file COPYING. If not, see
20 <https://www.gnu.org/licenses/>.
21 
22 */
23 
24 #if defined (HAVE_CONFIG_H)
25 # include "config.h"
26 #endif
27 
28 #include <cassert>
29 
30 #include <algorithm>
31 
32 #include "CMatrix.h"
33 #include "dDiagMatrix.h"
34 #include "dMatrix.h"
35 #include "fCMatrix.h"
36 #include "fDiagMatrix.h"
37 #include "fMatrix.h"
38 #include "lo-error.h"
39 #include "lo-lapack-proto.h"
40 #include "svd.h"
41 
42 namespace octave
43 {
44  namespace math
45  {
46  template <typename T>
47  T
49  {
50  if (m_type == svd::Type::sigma_only)
51  (*current_liboctave_error_handler)
52  ("svd: U not computed because type == svd::sigma_only");
53 
54  return left_sm;
55  }
56 
57  template <typename T>
58  T
60  {
61  if (m_type == svd::Type::sigma_only)
62  (*current_liboctave_error_handler)
63  ("svd: V not computed because type == svd::sigma_only");
64 
65  return right_sm;
66  }
67 
68  // GESVD specializations
69 
70 #define GESVD_REAL_STEP(f, F) \
71  F77_XFCN (f, F, (F77_CONST_CHAR_ARG2 (&jobu, 1), \
72  F77_CONST_CHAR_ARG2 (&jobv, 1), \
73  m, n, tmp_data, m1, s_vec, u, m1, vt, \
74  nrow_vt1, work.data (), lwork, info \
75  F77_CHAR_ARG_LEN (1) \
76  F77_CHAR_ARG_LEN (1)))
77 
78 #define GESVD_COMPLEX_STEP(f, F, CMPLX_ARG) \
79  F77_XFCN (f, F, (F77_CONST_CHAR_ARG2 (&jobu, 1), \
80  F77_CONST_CHAR_ARG2 (&jobv, 1), \
81  m, n, CMPLX_ARG (tmp_data), \
82  m1, s_vec, CMPLX_ARG (u), m1, \
83  CMPLX_ARG (vt), nrow_vt1, \
84  CMPLX_ARG (work.data ()), \
85  lwork, rwork.data (), info \
86  F77_CHAR_ARG_LEN (1) \
87  F77_CHAR_ARG_LEN (1)))
88 
89  // DGESVD
90  template<>
91  void
92  svd<Matrix>::gesvd (char& jobu, char& jobv, F77_INT m, F77_INT n,
93  double *tmp_data, F77_INT m1, double *s_vec,
94  double *u, double *vt, F77_INT nrow_vt1,
95  std::vector<double>& work, F77_INT& lwork,
96  F77_INT& info)
97  {
98  GESVD_REAL_STEP (dgesvd, DGESVD);
99 
100  lwork = static_cast<F77_INT> (work[0]);
101  work.reserve (lwork);
102 
103  GESVD_REAL_STEP (dgesvd, DGESVD);
104  }
105 
106  // SGESVD
107  template<>
108  void
109  svd<FloatMatrix>::gesvd (char& jobu, char& jobv, F77_INT m, F77_INT n,
110  float *tmp_data, F77_INT m1, float *s_vec,
111  float *u, float *vt, F77_INT nrow_vt1,
112  std::vector<float>& work, F77_INT& lwork,
113  F77_INT& info)
114  {
115  GESVD_REAL_STEP (sgesvd, SGESVD);
116 
117  lwork = static_cast<F77_INT> (work[0]);
118  work.reserve (lwork);
119 
120  GESVD_REAL_STEP (sgesvd, SGESVD);
121  }
122 
123  // ZGESVD
124  template<>
125  void
126  svd<ComplexMatrix>::gesvd (char& jobu, char& jobv, F77_INT m, F77_INT n,
127  Complex *tmp_data, F77_INT m1, double *s_vec,
128  Complex *u, Complex *vt, F77_INT nrow_vt1,
129  std::vector<Complex>& work, F77_INT& lwork,
130  F77_INT& info)
131  {
132  std::vector<double> rwork (5 * std::max (m, n));
133 
134  GESVD_COMPLEX_STEP (zgesvd, ZGESVD, F77_DBLE_CMPLX_ARG);
135 
136  lwork = static_cast<F77_INT> (work[0].real ());
137  work.reserve (lwork);
138 
139  GESVD_COMPLEX_STEP (zgesvd, ZGESVD, F77_DBLE_CMPLX_ARG);
140  }
141 
142  // CGESVD
143  template<>
144  void
145  svd<FloatComplexMatrix>::gesvd (char& jobu, char& jobv, F77_INT m,
146  F77_INT n, FloatComplex *tmp_data,
147  F77_INT m1, float *s_vec, FloatComplex *u,
148  FloatComplex *vt, F77_INT nrow_vt1,
149  std::vector<FloatComplex>& work,
150  F77_INT& lwork, F77_INT& info)
151  {
152  std::vector<float> rwork (5 * std::max (m, n));
153 
154  GESVD_COMPLEX_STEP (cgesvd, CGESVD, F77_CMPLX_ARG);
155 
156  lwork = static_cast<F77_INT> (work[0].real ());
157  work.reserve (lwork);
158 
159  GESVD_COMPLEX_STEP (cgesvd, CGESVD, F77_CMPLX_ARG);
160  }
161 
162 #undef GESVD_REAL_STEP
163 #undef GESVD_COMPLEX_STEP
164 
165  // GESDD specializations
166 
167 #define GESDD_REAL_STEP(f, F) \
168  F77_XFCN (f, F, (F77_CONST_CHAR_ARG2 (&jobz, 1), \
169  m, n, tmp_data, m1, s_vec, u, m1, vt, nrow_vt1, \
170  work.data (), lwork, iwork, info \
171  F77_CHAR_ARG_LEN (1)))
172 
173 #define GESDD_COMPLEX_STEP(f, F, CMPLX_ARG) \
174  F77_XFCN (f, F, (F77_CONST_CHAR_ARG2 (&jobz, 1), m, n, \
175  CMPLX_ARG (tmp_data), m1, \
176  s_vec, CMPLX_ARG (u), m1, \
177  CMPLX_ARG (vt), nrow_vt1, \
178  CMPLX_ARG (work.data ()), lwork, \
179  rwork.data (), iwork, info \
180  F77_CHAR_ARG_LEN (1)))
181 
182  // DGESDD
183  template<>
184  void
185  svd<Matrix>::gesdd (char& jobz, F77_INT m, F77_INT n, double *tmp_data,
186  F77_INT m1, double *s_vec, double *u, double *vt,
187  F77_INT nrow_vt1, std::vector<double>& work,
188  F77_INT& lwork, F77_INT *iwork, F77_INT& info)
189  {
190  GESDD_REAL_STEP (dgesdd, DGESDD);
191 
192  lwork = static_cast<F77_INT> (work[0]);
193  work.reserve (lwork);
194 
195  GESDD_REAL_STEP (dgesdd, DGESDD);
196  }
197 
198  // SGESDD
199  template<>
200  void
201  svd<FloatMatrix>::gesdd (char& jobz, F77_INT m, F77_INT n, float *tmp_data,
202  F77_INT m1, float *s_vec, float *u, float *vt,
203  F77_INT nrow_vt1, std::vector<float>& work,
204  F77_INT& lwork, F77_INT* iwork, F77_INT& info)
205  {
206  GESDD_REAL_STEP (sgesdd, SGESDD);
207 
208  lwork = static_cast<F77_INT> (work[0]);
209  work.reserve (lwork);
210 
211  GESDD_REAL_STEP (sgesdd, SGESDD);
212  }
213 
214  // ZGESDD
215  template<>
216  void
218  Complex *tmp_data, F77_INT m1, double *s_vec,
219  Complex *u, Complex *vt, F77_INT nrow_vt1,
220  std::vector<Complex>& work, F77_INT& lwork,
221  F77_INT *iwork, F77_INT& info)
222  {
223 
224  F77_INT min_mn = std::min (m, n);
225  F77_INT max_mn = std::max (m, n);
226 
227  F77_INT lrwork;
228  if (jobz == 'N')
229  lrwork = 7*min_mn;
230  else
231  lrwork = min_mn * std::max (5*min_mn+5, 2*max_mn+2*min_mn+1);
232 
233  std::vector<double> rwork (lrwork);
234 
235  GESDD_COMPLEX_STEP (zgesdd, ZGESDD, F77_DBLE_CMPLX_ARG);
236 
237  lwork = static_cast<F77_INT> (work[0].real ());
238  work.reserve (lwork);
239 
240  GESDD_COMPLEX_STEP (zgesdd, ZGESDD, F77_DBLE_CMPLX_ARG);
241  }
242 
243  // CGESDD
244  template<>
245  void
247  FloatComplex *tmp_data, F77_INT m1,
248  float *s_vec, FloatComplex *u,
249  FloatComplex *vt, F77_INT nrow_vt1,
250  std::vector<FloatComplex>& work,
251  F77_INT& lwork, F77_INT *iwork,
252  F77_INT& info)
253  {
254  F77_INT min_mn = std::min (m, n);
255  F77_INT max_mn = std::max (m, n);
256 
257  F77_INT lrwork;
258  if (jobz == 'N')
259  lrwork = 7*min_mn;
260  else
261  lrwork = min_mn * std::max (5*min_mn+5, 2*max_mn+2*min_mn+1);
262  std::vector<float> rwork (lrwork);
263 
264  GESDD_COMPLEX_STEP (cgesdd, CGESDD, F77_CMPLX_ARG);
265 
266  lwork = static_cast<F77_INT> (work[0].real ());
267  work.reserve (lwork);
268 
269  GESDD_COMPLEX_STEP (cgesdd, CGESDD, F77_CMPLX_ARG);
270  }
271 
272 #undef GESDD_REAL_STEP
273 #undef GESDD_COMPLEX_STEP
274 
275  template<typename T>
277  svd::Driver driver)
278  : m_type (type), m_driver (driver), left_sm (), sigma (), right_sm ()
279  {
280  F77_INT info;
281 
282  F77_INT m = to_f77_int (a.rows ());
283  F77_INT n = to_f77_int (a.cols ());
284 
285  if (m == 0 || n == 0)
286  {
287  switch (m_type)
288  {
289  case svd::Type::std:
290  left_sm = T (m, m, 0);
291  for (F77_INT i = 0; i < m; i++)
292  left_sm.xelem (i, i) = 1;
293  sigma = DM_T (m, n);
294  right_sm = T (n, n, 0);
295  for (F77_INT i = 0; i < n; i++)
296  right_sm.xelem (i, i) = 1;
297  break;
298 
299  case svd::Type::economy:
300  left_sm = T (m, 0, 0);
301  sigma = DM_T (0, 0);
302  right_sm = T (0, n, 0);
303  break;
304 
306  default:
307  sigma = DM_T (0, 1);
308  break;
309  }
310  return;
311  }
312 
313  T atmp = a;
314  P *tmp_data = atmp.fortran_vec ();
315 
316  F77_INT min_mn = (m < n ? m : n);
317 
318  char jobu = 'A';
319  char jobv = 'A';
320 
321  F77_INT ncol_u = m;
322  F77_INT nrow_vt = n;
323  F77_INT nrow_s = m;
324  F77_INT ncol_s = n;
325 
326  switch (m_type)
327  {
328  case svd::Type::economy:
329  jobu = jobv = 'S';
330  ncol_u = nrow_vt = nrow_s = ncol_s = min_mn;
331  break;
332 
334 
335  // Note: for this case, both jobu and jobv should be 'N', but there
336  // seems to be a bug in dgesvd from Lapack V2.0. To demonstrate the
337  // bug, set both jobu and jobv to 'N' and find the singular values of
338  // [eye(3), eye(3)]. The result is [-sqrt(2), -sqrt(2), -sqrt(2)].
339  //
340  // For Lapack 3.0, this problem seems to be fixed.
341 
342  jobu = jobv = 'N';
343  ncol_u = nrow_vt = 1;
344  break;
345 
346  default:
347  break;
348  }
349 
350  if (! (jobu == 'N' || jobu == 'O'))
351  left_sm.resize (m, ncol_u);
352 
353  P *u = left_sm.fortran_vec ();
354 
355  sigma.resize (nrow_s, ncol_s);
356  DM_P *s_vec = sigma.fortran_vec ();
357 
358  if (! (jobv == 'N' || jobv == 'O'))
359  right_sm.resize (nrow_vt, n);
360 
361  P *vt = right_sm.fortran_vec ();
362 
363  // Query _GESVD for the correct dimension of WORK.
364 
365  F77_INT lwork = -1;
366 
367  std::vector<P> work (1);
368 
369  F77_INT m1 = std::max (m, static_cast<F77_INT> (1));
370  F77_INT nrow_vt1 = std::max (nrow_vt, static_cast<F77_INT> (1));
371 
373  gesvd (jobu, jobv, m, n, tmp_data, m1, s_vec, u, vt, nrow_vt1,
374  work, lwork, info);
375  else if (m_driver == svd::Driver::GESDD)
376  {
377  assert (jobu == jobv);
378  char jobz = jobu;
379 
380  std::vector<F77_INT> iwork (8 * std::min (m, n));
381 
382  gesdd (jobz, m, n, tmp_data, m1, s_vec, u, vt, nrow_vt1,
383  work, lwork, iwork.data (), info);
384  }
385  else
386  (*current_liboctave_error_handler) ("svd: unknown driver");
387 
388  if (! (jobv == 'N' || jobv == 'O'))
389  right_sm = right_sm.hermitian ();
390  }
391 
392  // Instantiations we need.
393 
394  template class svd<Matrix>;
395 
396  template class svd<FloatMatrix>;
397 
398  template class svd<ComplexMatrix>;
399 
400  template class svd<FloatComplexMatrix>;
401  }
402 }
OCTAVE_NORETURN liboctave_error_handler current_liboctave_error_handler
Definition: lo-error.c:38
#define F77_DBLE_CMPLX_ARG(x)
Definition: f77-fcn.h:315
#define GESDD_COMPLEX_STEP(f, F, CMPLX_ARG)
Definition: svd.cc:173
DM_T sigma
Definition: svd.h:99
#define GESVD_COMPLEX_STEP(f, F, CMPLX_ARG)
Definition: svd.cc:78
svd::Type m_type
Definition: svd.h:95
T::element_type P
Definition: svd.h:92
u
Definition: lu.cc:138
svd(void)
Definition: svd.h:56
DM_T::element_type DM_P
Definition: svd.h:93
T right_singular_matrix(void) const
Definition: svd.cc:59
calling an anonymous function involves an overhead quite comparable to the overhead of an m file function Passing a handle to a built in function is because the interpreter is not involved in the internal loop For a
Definition: cellfun.cc:400
T::real_diag_matrix_type DM_T
Definition: svd.h:41
idx type
Definition: ov.cc:3114
#define GESDD_REAL_STEP(f, F)
Definition: svd.cc:167
void gesdd(char &jobz, octave_f77_int_type m, octave_f77_int_type n, P *tmp_data, octave_f77_int_type m1, DM_P *s_vec, P *u, P *vt, octave_f77_int_type nrow_vt1, std::vector< P > &work, octave_f77_int_type &lwork, octave_f77_int_type *iwork, octave_f77_int_type &info)
void gesvd(char &jobu, char &jobv, octave_f77_int_type m, octave_f77_int_type n, P *tmp_data, octave_f77_int_type m1, DM_P *s_vec, P *u, P *vt, octave_f77_int_type nrow_vt1, std::vector< P > &work, octave_f77_int_type &lwork, octave_f77_int_type &info)
charNDArray max(char d, const charNDArray &m)
Definition: chNDArray.cc:227
#define F77_CMPLX_ARG(x)
Definition: f77-fcn.h:309
octave_f77_int_type F77_INT
Definition: f77-fcn.h:305
#define GESVD_REAL_STEP(f, F)
Definition: svd.cc:70
for i
Definition: data.cc:5264
std::complex< float > FloatComplex
Definition: oct-cmplx.h:32
std::complex< double > Complex
Definition: oct-cmplx.h:31
charNDArray min(char d, const charNDArray &m)
Definition: chNDArray.cc:204
svd::Driver m_driver
Definition: svd.h:96
T left_singular_matrix(void) const
Definition: svd.cc:48