Mercurial > hg > octave-lyh

comparison src/DLD-FUNCTIONS/svd.cc @ 7789:82be108cc558

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
First attempt at single precision tyeps * * * corrections to qrupdate single precision routines * * * prefer demotion to single over promotion to double * * * Add single precision support to log2 function * * * Trivial PROJECT file update * * * Cache optimized hermitian/transpose methods * * * Add tests for tranpose/hermitian and ChangeLog entry for new transpose code
author David Bateman <dbateman@free.fr>
date Sun, 27 Apr 2008 22:34:17 +0200
parents a1dbe9d80eee
children 87865ed7405f
comparison
equal deleted inserted replaced
7788:45f5faba05a2 7789:82be108cc558
25 #include <config.h> 25 #include <config.h>
26 #endif 26 #endif
27 27
28 #include "CmplxSVD.h" 28 #include "CmplxSVD.h"
29 #include "dbleSVD.h" 29 #include "dbleSVD.h"
30 #include "fCmplxSVD.h"
31 #include "floatSVD.h"
30 32
31 #include "defun-dld.h" 33 #include "defun-dld.h"
32 #include "error.h" 34 #include "error.h"
33 #include "gripes.h" 35 #include "gripes.h"
34 #include "oct-obj.h" 36 #include "oct-obj.h"
130 octave_value arg = args(0); 132 octave_value arg = args(0);
131 133
132 octave_idx_type nr = arg.rows (); 134 octave_idx_type nr = arg.rows ();
133 octave_idx_type nc = arg.columns (); 135 octave_idx_type nc = arg.columns ();
134 136
137 bool isfloat = arg.is_single_type ();
138
135 if (nr == 0 || nc == 0) 139 if (nr == 0 || nc == 0)
136 { 140 {
137 if (nargout == 3) 141 if (isfloat)
138 { 142 {
139 retval(3) = identity_matrix (nr, nr); 143 if (nargout == 3)
140 retval(2) = Matrix (nr, nc); 144 {
141 retval(1) = identity_matrix (nc, nc); 145 retval(3) = float_identity_matrix (nr, nr);
146 retval(2) = FloatMatrix (nr, nc);
147 retval(1) = float_identity_matrix (nc, nc);
148 }
149 else
150 retval(0) = FloatMatrix (0, 1);
142 } 151 }
143 else 152 else
144 retval(0) = Matrix (0, 1); 153 {
154 if (nargout == 3)
155 {
156 retval(3) = identity_matrix (nr, nr);
157 retval(2) = Matrix (nr, nc);
158 retval(1) = identity_matrix (nc, nc);
159 }
160 else
161 retval(0) = Matrix (0, 1);
162 }
145 } 163 }
146 else 164 else
147 { 165 {
148 SVD::type type = ((nargout == 0 || nargout == 1) 166 SVD::type type = ((nargout == 0 || nargout == 1)
149 ? SVD::sigma_only 167 ? SVD::sigma_only
150 : (nargin == 2) ? SVD::economy : SVD::std); 168 : (nargin == 2) ? SVD::economy : SVD::std);
151 169
152 if (arg.is_real_type ()) 170 if (isfloat)
153 { 171 {
154 Matrix tmp = arg.matrix_value (); 172 if (arg.is_real_type ())
155 173 {
156 if (! error_state) 174 FloatMatrix tmp = arg.float_matrix_value ();
157 { 175
158 if (tmp.any_element_is_inf_or_nan ()) 176 if (! error_state)
159 { 177 {
160 error ("svd: cannot take SVD of matrix containing Inf or NaN values"); 178 if (tmp.any_element_is_inf_or_nan ())
161 return retval; 179 {
180 error ("svd: cannot take SVD of matrix containing Inf or NaN values");
181 return retval;
182 }
183
184 FloatSVD result (tmp, type);
185
186 FloatDiagMatrix sigma = result.singular_values ();
187
188 if (nargout == 0 || nargout == 1)
189 {
190 retval(0) = sigma.diag ();
191 }
192 else
193 {
194 retval(2) = result.right_singular_matrix ();
195 retval(1) = sigma;
196 retval(0) = result.left_singular_matrix ();
197 }
162 } 198 }
163 199 }
164 SVD result (tmp, type); 200 else if (arg.is_complex_type ())
165 201 {
166 DiagMatrix sigma = result.singular_values (); 202 FloatComplexMatrix ctmp = arg.float_complex_matrix_value ();
167 203
168 if (nargout == 0 || nargout == 1) 204 if (! error_state)
169 { 205 {
170 retval(0) = sigma.diag (); 206 if (ctmp.any_element_is_inf_or_nan ())
171 } 207 {
172 else 208 error ("svd: cannot take SVD of matrix containing Inf or NaN values");
173 { 209 return retval;
174 retval(2) = result.right_singular_matrix (); 210 }
175 retval(1) = sigma; 211
176 retval(0) = result.left_singular_matrix (); 212 FloatComplexSVD result (ctmp, type);
177 } 213
178 } 214 FloatDiagMatrix sigma = result.singular_values ();
179 } 215
180 else if (arg.is_complex_type ()) 216 if (nargout == 0 || nargout == 1)
181 { 217 {
182 ComplexMatrix ctmp = arg.complex_matrix_value (); 218 retval(0) = sigma.diag ();
183 219 }
184 if (! error_state) 220 else
185 { 221 {
186 if (ctmp.any_element_is_inf_or_nan ()) 222 retval(2) = result.right_singular_matrix ();
187 { 223 retval(1) = sigma;
188 error ("svd: cannot take SVD of matrix containing Inf or NaN values"); 224 retval(0) = result.left_singular_matrix ();
189 return retval; 225 }
190 }
191
192 ComplexSVD result (ctmp, type);
193
194 DiagMatrix sigma = result.singular_values ();
195
196 if (nargout == 0 || nargout == 1)
197 {
198 retval(0) = sigma.diag ();
199 }
200 else
201 {
202 retval(2) = result.right_singular_matrix ();
203 retval(1) = sigma;
204 retval(0) = result.left_singular_matrix ();
205 } 226 }
206 } 227 }
207 } 228 }
208 else 229 else
209 { 230 {
210 gripe_wrong_type_arg ("svd", arg); 231 if (arg.is_real_type ())
211 return retval; 232 {
233 Matrix tmp = arg.matrix_value ();
234
235 if (! error_state)
236 {
237 if (tmp.any_element_is_inf_or_nan ())
238 {
239 error ("svd: cannot take SVD of matrix containing Inf or NaN values");
240 return retval;
241 }
242
243 SVD result (tmp, type);
244
245 DiagMatrix sigma = result.singular_values ();
246
247 if (nargout == 0 || nargout == 1)
248 {
249 retval(0) = sigma.diag ();
250 }
251 else
252 {
253 retval(2) = result.right_singular_matrix ();
254 retval(1) = sigma;
255 retval(0) = result.left_singular_matrix ();
256 }
257 }
258 }
259 else if (arg.is_complex_type ())
260 {
261 ComplexMatrix ctmp = arg.complex_matrix_value ();
262
263 if (! error_state)
264 {
265 if (ctmp.any_element_is_inf_or_nan ())
266 {
267 error ("svd: cannot take SVD of matrix containing Inf or NaN values");
268 return retval;
269 }
270
271 ComplexSVD result (ctmp, type);
272
273 DiagMatrix sigma = result.singular_values ();
274
275 if (nargout == 0 || nargout == 1)
276 {
277 retval(0) = sigma.diag ();
278 }
279 else
280 {
281 retval(2) = result.right_singular_matrix ();
282 retval(1) = sigma;
283 retval(0) = result.left_singular_matrix ();
284 }
285 }
286 }
287 else
288 {
289 gripe_wrong_type_arg ("svd", arg);
290 return retval;
291 }
212 } 292 }
213 } 293 }
214 294
215 return retval; 295 return retval;
216 } 296 }