Mercurial > octave-nkf

comparison liboctave/dbleSVD.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
27 27
28 #include <iostream> 28 #include <iostream>
29 29
30 #include "dbleSVD.h" 30 #include "dbleSVD.h"
31 #include "f77-fcn.h" 31 #include "f77-fcn.h"
32 #include "oct-locbuf.h"
32 33
33 extern "C" 34 extern "C"
34 { 35 {
35 F77_RET_T 36 F77_RET_T
36 F77_FUNC (dgesvd, DGESVD) (F77_CONST_CHAR_ARG_DECL, 37 F77_FUNC (dgesvd, DGESVD) (F77_CONST_CHAR_ARG_DECL,
39 const octave_idx_type&, double*, double*, 40 const octave_idx_type&, double*, double*,
40 const octave_idx_type&, double*, const octave_idx_type&, 41 const octave_idx_type&, double*, const octave_idx_type&,
41 double*, const octave_idx_type&, octave_idx_type& 42 double*, const octave_idx_type&, octave_idx_type&
42 F77_CHAR_ARG_LEN_DECL 43 F77_CHAR_ARG_LEN_DECL
43 F77_CHAR_ARG_LEN_DECL); 44 F77_CHAR_ARG_LEN_DECL);
45
46 F77_RET_T
47 F77_FUNC (dgesdd, DGESDD) (F77_CONST_CHAR_ARG_DECL,
48 const octave_idx_type&, const octave_idx_type&, double*,
49 const octave_idx_type&, double*, double*,
50 const octave_idx_type&, double*, const octave_idx_type&,
51 double*, const octave_idx_type&, octave_idx_type *, octave_idx_type&
52 F77_CHAR_ARG_LEN_DECL);
44 } 53 }
45 54
46 Matrix 55 Matrix
47 SVD::left_singular_matrix (void) const 56 SVD::left_singular_matrix (void) const
48 { 57 {
68 else 77 else
69 return right_sm; 78 return right_sm;
70 } 79 }
71 80
72 octave_idx_type 81 octave_idx_type
73 SVD::init (const Matrix& a, SVD::type svd_type) 82 SVD::init (const Matrix& a, SVD::type svd_type, SVD::driver svd_driver)
74 { 83 {
75 octave_idx_type info; 84 octave_idx_type info;
76 85
77 octave_idx_type m = a.rows (); 86 octave_idx_type m = a.rows ();
78 octave_idx_type n = a.cols (); 87 octave_idx_type n = a.cols ();
138 Array<double> work (1, 1); 147 Array<double> work (1, 1);
139 148
140 octave_idx_type one = 1; 149 octave_idx_type one = 1;
141 octave_idx_type m1 = std::max (m, one), nrow_vt1 = std::max (nrow_vt, one); 150 octave_idx_type m1 = std::max (m, one), nrow_vt1 = std::max (nrow_vt, one);
142 151
143 F77_XFCN (dgesvd, DGESVD, (F77_CONST_CHAR_ARG2 (&jobu, 1), 152 if (svd_driver == SVD::GESVD)
144 F77_CONST_CHAR_ARG2 (&jobv, 1), 153 {
145 m, n, tmp_data, m1, s_vec, u, m1, vt, 154 F77_XFCN (dgesvd, DGESVD, (F77_CONST_CHAR_ARG2 (&jobu, 1),
146 nrow_vt1, work.fortran_vec (), lwork, info 155 F77_CONST_CHAR_ARG2 (&jobv, 1),
147 F77_CHAR_ARG_LEN (1) 156 m, n, tmp_data, m1, s_vec, u, m1, vt,
148 F77_CHAR_ARG_LEN (1))); 157 nrow_vt1, work.fortran_vec (), lwork, info
149 158 F77_CHAR_ARG_LEN (1)
150 lwork = static_cast<octave_idx_type> (work(0)); 159 F77_CHAR_ARG_LEN (1)));
151 work.resize (lwork, 1); 160
152 161 lwork = static_cast<octave_idx_type> (work(0));
153 F77_XFCN (dgesvd, DGESVD, (F77_CONST_CHAR_ARG2 (&jobu, 1), 162 work.resize (lwork, 1);
154 F77_CONST_CHAR_ARG2 (&jobv, 1), 163
155 m, n, tmp_data, m1, s_vec, u, m1, vt, 164 F77_XFCN (dgesvd, DGESVD, (F77_CONST_CHAR_ARG2 (&jobu, 1),
156 nrow_vt1, work.fortran_vec (), lwork, info 165 F77_CONST_CHAR_ARG2 (&jobv, 1),
157 F77_CHAR_ARG_LEN (1) 166 m, n, tmp_data, m1, s_vec, u, m1, vt,
158 F77_CHAR_ARG_LEN (1))); 167 nrow_vt1, work.fortran_vec (), lwork, info
168 F77_CHAR_ARG_LEN (1)
169 F77_CHAR_ARG_LEN (1)));
170
171 }
172 else if (svd_driver == SVD::GESDD)
173 {
174 assert (jobu == jobv);
175 char jobz = jobu;
176 OCTAVE_LOCAL_BUFFER (octave_idx_type, iwork, 8*min_mn);
177
178 F77_XFCN (dgesdd, DGESDD, (F77_CONST_CHAR_ARG2 (&jobz, 1),
179 m, n, tmp_data, m1, s_vec, u, m1, vt,
180 nrow_vt1, work.fortran_vec (), lwork, iwork, info
181 F77_CHAR_ARG_LEN (1)));
182
183 lwork = static_cast<octave_idx_type> (work(0));
184 work.resize (lwork, 1);
185
186 F77_XFCN (dgesdd, DGESDD, (F77_CONST_CHAR_ARG2 (&jobz, 1),
187 m, n, tmp_data, m1, s_vec, u, m1, vt,
188 nrow_vt1, work.fortran_vec (), lwork, iwork, info
189 F77_CHAR_ARG_LEN (1)));
190
191 }
192 else
193 assert (0); // impossible
159 194
160 if (! (jobv == 'N' || jobv == 'O')) 195 if (! (jobv == 'N' || jobv == 'O'))
161 right_sm = right_sm.transpose (); 196 right_sm = right_sm.transpose ();
162 197
163 return info; 198 return info;