Mercurial > octave-nkf

comparison liboctave/fCmplxSVD.cc @ 10601:3ce0c530a9c9

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
implement svd_driver
author Jaroslav Hajek <highegg@gmail.com>
date Mon, 03 May 2010 13:21:35 +0200
parents 12884915a8e4
children 8a5e980da6aa
comparison
equal deleted inserted replaced
10600:036bdc2d0af0 10601:3ce0c530a9c9
26 #endif 26 #endif
27 27
28 #include "fCmplxSVD.h" 28 #include "fCmplxSVD.h"
29 #include "f77-fcn.h" 29 #include "f77-fcn.h"
30 #include "lo-error.h" 30 #include "lo-error.h"
31 #include "oct-locbuf.h"
31 32
32 extern "C" 33 extern "C"
33 { 34 {
34 F77_RET_T 35 F77_RET_T
35 F77_FUNC (cgesvd, CGESVD) (F77_CONST_CHAR_ARG_DECL, 36 F77_FUNC (cgesvd, CGESVD) (F77_CONST_CHAR_ARG_DECL,
38 const octave_idx_type&, float*, FloatComplex*, const octave_idx_type&, 39 const octave_idx_type&, float*, FloatComplex*, const octave_idx_type&,
39 FloatComplex*, const octave_idx_type&, FloatComplex*, const octave_idx_type&, 40 FloatComplex*, const octave_idx_type&, FloatComplex*, const octave_idx_type&,
40 float*, octave_idx_type& 41 float*, octave_idx_type&
41 F77_CHAR_ARG_LEN_DECL 42 F77_CHAR_ARG_LEN_DECL
42 F77_CHAR_ARG_LEN_DECL); 43 F77_CHAR_ARG_LEN_DECL);
44
45 F77_RET_T
46 F77_FUNC (cgesdd, CGESDD) (F77_CONST_CHAR_ARG_DECL,
47 const octave_idx_type&, const octave_idx_type&, FloatComplex*,
48 const octave_idx_type&, float*, FloatComplex*, const octave_idx_type&,
49 FloatComplex*, const octave_idx_type&, FloatComplex*, const octave_idx_type&,
50 float*, octave_idx_type *, octave_idx_type&
51 F77_CHAR_ARG_LEN_DECL);
43 } 52 }
44 53
45 FloatComplexMatrix 54 FloatComplexMatrix
46 FloatComplexSVD::left_singular_matrix (void) const 55 FloatComplexSVD::left_singular_matrix (void) const
47 { 56 {
67 else 76 else
68 return right_sm; 77 return right_sm;
69 } 78 }
70 79
71 octave_idx_type 80 octave_idx_type
72 FloatComplexSVD::init (const FloatComplexMatrix& a, SVD::type svd_type) 81 FloatComplexSVD::init (const FloatComplexMatrix& a, SVD::type svd_type, SVD::driver svd_driver)
73 { 82 {
74 octave_idx_type info; 83 octave_idx_type info;
75 84
76 octave_idx_type m = a.rows (); 85 octave_idx_type m = a.rows ();
77 octave_idx_type n = a.cols (); 86 octave_idx_type n = a.cols ();
142 Array<FloatComplex> work (1, 1); 151 Array<FloatComplex> work (1, 1);
143 152
144 octave_idx_type one = 1; 153 octave_idx_type one = 1;
145 octave_idx_type m1 = std::max (m, one), nrow_vt1 = std::max (nrow_vt, one); 154 octave_idx_type m1 = std::max (m, one), nrow_vt1 = std::max (nrow_vt, one);
146 155
147 F77_XFCN (cgesvd, CGESVD, (F77_CONST_CHAR_ARG2 (&jobu, 1), 156 if (svd_driver == SVD::GESVD)
148 F77_CONST_CHAR_ARG2 (&jobv, 1), 157 {
149 m, n, tmp_data, m1, s_vec, u, m1, vt, 158 F77_XFCN (cgesvd, CGESVD, (F77_CONST_CHAR_ARG2 (&jobu, 1),
150 nrow_vt1, work.fortran_vec (), lwork, 159 F77_CONST_CHAR_ARG2 (&jobv, 1),
151 rwork.fortran_vec (), info 160 m, n, tmp_data, m1, s_vec, u, m1, vt,
152 F77_CHAR_ARG_LEN (1) 161 nrow_vt1, work.fortran_vec (), lwork,
153 F77_CHAR_ARG_LEN (1))); 162 rwork.fortran_vec (), info
154 163 F77_CHAR_ARG_LEN (1)
155 lwork = static_cast<octave_idx_type> (work(0).real ()); 164 F77_CHAR_ARG_LEN (1)));
156 work.resize (lwork, 1); 165
157 166 lwork = static_cast<octave_idx_type> (work(0).real ());
158 F77_XFCN (cgesvd, CGESVD, (F77_CONST_CHAR_ARG2 (&jobu, 1), 167 work.resize (lwork, 1);
159 F77_CONST_CHAR_ARG2 (&jobv, 1), 168
160 m, n, tmp_data, m1, s_vec, u, m1, vt, 169 F77_XFCN (cgesvd, CGESVD, (F77_CONST_CHAR_ARG2 (&jobu, 1),
161 nrow_vt1, work.fortran_vec (), lwork, 170 F77_CONST_CHAR_ARG2 (&jobv, 1),
162 rwork.fortran_vec (), info 171 m, n, tmp_data, m1, s_vec, u, m1, vt,
163 F77_CHAR_ARG_LEN (1) 172 nrow_vt1, work.fortran_vec (), lwork,
164 F77_CHAR_ARG_LEN (1))); 173 rwork.fortran_vec (), info
174 F77_CHAR_ARG_LEN (1)
175 F77_CHAR_ARG_LEN (1)));
176 }
177 else if (svd_driver == SVD::GESDD)
178 {
179 assert (jobu == jobv);
180 char jobz = jobu;
181 OCTAVE_LOCAL_BUFFER (octave_idx_type, iwork, 8*min_mn);
182
183 F77_XFCN (cgesdd, CGESDD, (F77_CONST_CHAR_ARG2 (&jobz, 1),
184 m, n, tmp_data, m1, s_vec, u, m1, vt,
185 nrow_vt1, work.fortran_vec (), lwork,
186 rwork.fortran_vec (), iwork, info
187 F77_CHAR_ARG_LEN (1)));
188
189 lwork = static_cast<octave_idx_type> (work(0).real ());
190 work.resize (lwork, 1);
191
192 F77_XFCN (cgesdd, CGESDD, (F77_CONST_CHAR_ARG2 (&jobz, 1),
193 m, n, tmp_data, m1, s_vec, u, m1, vt,
194 nrow_vt1, work.fortran_vec (), lwork,
195 rwork.fortran_vec (), iwork, info
196 F77_CHAR_ARG_LEN (1)));
197 }
198 else
199 assert (0); // impossible
165 200
166 if (! (jobv == 'N' || jobv == 'O')) 201 if (! (jobv == 'N' || jobv == 'O'))
167 right_sm = right_sm.hermitian (); 202 right_sm = right_sm.hermitian ();
168 203
169 return info; 204 return info;