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
CmplxSVD.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 "CmplxSVD.h"
28 #include "f77-fcn.h"
29 #include "lo-error.h"
30 #include "oct-locbuf.h"
31 
32 extern "C"
33 {
34  F77_RET_T
35  F77_FUNC (zgesvd, ZGESVD) (F77_CONST_CHAR_ARG_DECL,
37  const octave_idx_type&, const octave_idx_type&,
38  Complex*, const octave_idx_type&,
39  double*, Complex*, const octave_idx_type&,
40  Complex*, const octave_idx_type&, Complex*,
41  const octave_idx_type&, double*, octave_idx_type&
44 
45  F77_RET_T
46  F77_FUNC (zgesdd, ZGESDD) (F77_CONST_CHAR_ARG_DECL,
47  const octave_idx_type&, const octave_idx_type&,
48  Complex*, const octave_idx_type&,
49  double*, Complex*, const octave_idx_type&,
50  Complex*, const octave_idx_type&, Complex*,
51  const octave_idx_type&, double*,
52  octave_idx_type *, octave_idx_type&
54 }
55 
58 {
60  {
61  (*current_liboctave_error_handler)
62  ("ComplexSVD: U not computed because type == SVD::sigma_only");
63  return ComplexMatrix ();
64  }
65  else
66  return left_sm;
67 }
68 
71 {
73  {
74  (*current_liboctave_error_handler)
75  ("ComplexSVD: V not computed because type == SVD::sigma_only");
76  return ComplexMatrix ();
77  }
78  else
79  return right_sm;
80 }
81 
84  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  ComplexMatrix atmp = a;
92  Complex *tmp_data = atmp.fortran_vec ();
93 
94  octave_idx_type min_mn = m < n ? m : n;
95  octave_idx_type max_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  Complex *u = left_sm.fortran_vec ();
136 
137  sigma.resize (nrow_s, ncol_s);
138  double *s_vec = sigma.fortran_vec ();
139 
140  if (! (jobv == 'N' || jobv == 'O'))
141  right_sm.resize (nrow_vt, n);
142 
143  Complex *vt = right_sm.fortran_vec ();
144 
145  // Query ZGESVD for the correct dimension of WORK.
146 
147  octave_idx_type lwork = -1;
148 
149  Array<Complex> 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  octave_idx_type lrwork = 5*max_mn;
158  Array<double> rwork (dim_vector (lrwork, 1));
159 
160  F77_XFCN (zgesvd, ZGESVD, (F77_CONST_CHAR_ARG2 (&jobu, 1),
161  F77_CONST_CHAR_ARG2 (&jobv, 1),
162  m, n, tmp_data, m1, s_vec, u, m1, vt,
163  nrow_vt1, work.fortran_vec (), lwork,
164  rwork.fortran_vec (), info
165  F77_CHAR_ARG_LEN (1)
166  F77_CHAR_ARG_LEN (1)));
167 
168  lwork = static_cast<octave_idx_type> (work(0).real ());
169  work.resize (dim_vector (lwork, 1));
170 
171  F77_XFCN (zgesvd, ZGESVD, (F77_CONST_CHAR_ARG2 (&jobu, 1),
172  F77_CONST_CHAR_ARG2 (&jobv, 1),
173  m, n, tmp_data, m1, s_vec, u, m1, vt,
174  nrow_vt1, work.fortran_vec (), lwork,
175  rwork.fortran_vec (), info
176  F77_CHAR_ARG_LEN (1)
177  F77_CHAR_ARG_LEN (1)));
178  }
179  else if (svd_driver == SVD::GESDD)
180  {
181  assert (jobu == jobv);
182  char jobz = jobu;
183 
184  octave_idx_type lrwork;
185  if (jobz == 'N')
186  lrwork = 7*min_mn;
187  else
188  lrwork = 5*min_mn*min_mn + 5*min_mn;
189  Array<double> rwork (dim_vector (lrwork, 1));
190 
191  OCTAVE_LOCAL_BUFFER (octave_idx_type, iwork, 8*min_mn);
192 
193  F77_XFCN (zgesdd, ZGESDD, (F77_CONST_CHAR_ARG2 (&jobz, 1),
194  m, n, tmp_data, m1, s_vec, u, m1, vt,
195  nrow_vt1, work.fortran_vec (), lwork,
196  rwork.fortran_vec (), iwork, info
197  F77_CHAR_ARG_LEN (1)));
198 
199  lwork = static_cast<octave_idx_type> (work(0).real ());
200  work.resize (dim_vector (lwork, 1));
201 
202  F77_XFCN (zgesdd, ZGESDD, (F77_CONST_CHAR_ARG2 (&jobz, 1),
203  m, n, tmp_data, m1, s_vec, u, m1, vt,
204  nrow_vt1, work.fortran_vec (), lwork,
205  rwork.fortran_vec (), iwork, info
206  F77_CHAR_ARG_LEN (1)));
207  }
208  else
209  assert (0); // impossible
210 
211  if (! (jobv == 'N' || jobv == 'O'))
213 
214  return info;
215 }