Mercurial > octave

comparison liboctave/numeric/svd.cc @ 31607:aac27ad79be6 stable

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
maint: Re-indent code after switch to using namespace macros. * build-env.h, build-env.in.cc, Cell.h, __betainc__.cc, __eigs__.cc, __ftp__.cc, __ichol__.cc, __ilu__.cc, __isprimelarge__.cc, __magick_read__.cc, __pchip_deriv__.cc, amd.cc, base-text-renderer.cc, base-text-renderer.h, besselj.cc, bitfcns.cc, bsxfun.cc, c-file-ptr-stream.h, call-stack.cc, call-stack.h, ccolamd.cc, cellfun.cc, chol.cc, colamd.cc, dasrt.cc, data.cc, debug.cc, defaults.cc, defaults.h, det.cc, display.cc, display.h, dlmread.cc, dynamic-ld.cc, dynamic-ld.h, ellipj.cc, environment.cc, environment.h, error.cc, error.h, errwarn.h, event-manager.cc, event-manager.h, event-queue.cc, event-queue.h, fcn-info.cc, fcn-info.h, fft.cc, fft2.cc, file-io.cc, filter.cc, find.cc, ft-text-renderer.cc, ft-text-renderer.h, gcd.cc, gl-render.cc, gl-render.h, gl2ps-print.cc, gl2ps-print.h, graphics-toolkit.cc, graphics-toolkit.h, graphics.cc, gsvd.cc, gtk-manager.cc, gtk-manager.h, help.cc, help.h, hook-fcn.cc, hook-fcn.h, input.cc, input.h, interpreter-private.cc, interpreter-private.h, interpreter.cc, interpreter.h, inv.cc, jsondecode.cc, jsonencode.cc, latex-text-renderer.cc, latex-text-renderer.h, load-path.cc, load-path.h, load-save.cc, load-save.h, lookup.cc, ls-hdf5.cc, ls-mat4.cc, ls-mat5.cc, lsode.cc, lu.cc, mappers.cc, matrix_type.cc, max.cc, mex.cc, mexproto.h, mxarray.h, mxtypes.in.h, oct-errno.in.cc, oct-hdf5-types.cc, oct-hist.cc, oct-hist.h, oct-map.cc, oct-map.h, oct-opengl.h, oct-prcstrm.h, oct-process.cc, oct-process.h, oct-stdstrm.h, oct-stream.cc, oct-stream.h, oct-strstrm.h, octave-default-image.h, ordqz.cc, ordschur.cc, pager.cc, pager.h, pinv.cc, pow2.cc, pr-output.cc, psi.cc, qr.cc, quadcc.cc, rand.cc, regexp.cc, settings.cc, settings.h, sighandlers.cc, sighandlers.h, sparse-xpow.cc, sqrtm.cc, stack-frame.cc, stack-frame.h, stream-euler.cc, strfns.cc, svd.cc, syminfo.cc, syminfo.h, symrcm.cc, symrec.cc, symrec.h, symscope.cc, symscope.h, symtab.cc, symtab.h, sysdep.cc, sysdep.h, text-engine.cc, text-engine.h, text-renderer.cc, text-renderer.h, time.cc, toplev.cc, typecast.cc, url-handle-manager.cc, url-handle-manager.h, urlwrite.cc, utils.cc, utils.h, variables.cc, variables.h, xdiv.cc, __delaunayn__.cc, __init_fltk__.cc, __init_gnuplot__.cc, __ode15__.cc, __voronoi__.cc, audioread.cc, convhulln.cc, gzip.cc, cdef-class.cc, cdef-class.h, cdef-fwd.h, cdef-manager.cc, cdef-manager.h, cdef-method.cc, cdef-method.h, cdef-object.cc, cdef-object.h, cdef-package.cc, cdef-package.h, cdef-property.cc, cdef-property.h, cdef-utils.cc, cdef-utils.h, ov-base-diag.cc, ov-base-int.cc, ov-base-mat.cc, ov-base-mat.h, ov-base-scalar.cc, ov-base.cc, ov-base.h, ov-bool-mat.cc, ov-bool-mat.h, ov-bool-sparse.cc, ov-bool.cc, ov-builtin.h, ov-cell.cc, ov-ch-mat.cc, ov-class.cc, ov-class.h, ov-classdef.cc, ov-classdef.h, ov-complex.cc, ov-cx-diag.cc, ov-cx-mat.cc, ov-cx-sparse.cc, ov-dld-fcn.cc, ov-dld-fcn.h, ov-fcn-handle.cc, ov-fcn-handle.h, ov-fcn.h, ov-float.cc, ov-flt-complex.cc, ov-flt-cx-diag.cc, ov-flt-cx-mat.cc, ov-flt-re-diag.cc, ov-flt-re-mat.cc, ov-flt-re-mat.h, ov-intx.h, ov-java.cc, ov-lazy-idx.cc, ov-legacy-range.cc, ov-magic-int.cc, ov-mex-fcn.cc, ov-mex-fcn.h, ov-null-mat.cc, ov-perm.cc, ov-range.cc, ov-re-diag.cc, ov-re-mat.cc, ov-re-mat.h, ov-re-sparse.cc, ov-scalar.cc, ov-str-mat.cc, ov-struct.cc, ov-typeinfo.cc, ov-typeinfo.h, ov-usr-fcn.cc, ov-usr-fcn.h, ov.cc, ov.h, ovl.h, octave.cc, octave.h, op-b-sbm.cc, op-bm-sbm.cc, op-cs-scm.cc, op-fm-fcm.cc, op-fs-fcm.cc, op-s-scm.cc, op-scm-cs.cc, op-scm-s.cc, op-sm-cs.cc, ops.h, anon-fcn-validator.cc, anon-fcn-validator.h, bp-table.cc, bp-table.h, comment-list.cc, comment-list.h, filepos.h, lex.h, oct-lvalue.cc, oct-lvalue.h, parse.h, profiler.cc, profiler.h, pt-anon-scopes.cc, pt-anon-scopes.h, pt-arg-list.cc, pt-arg-list.h, pt-args-block.cc, pt-args-block.h, pt-array-list.cc, pt-array-list.h, pt-assign.cc, pt-assign.h, pt-binop.cc, pt-binop.h, pt-bp.cc, pt-bp.h, pt-cbinop.cc, pt-cbinop.h, pt-cell.cc, pt-cell.h, pt-check.cc, pt-check.h, pt-classdef.cc, pt-classdef.h, pt-cmd.h, pt-colon.cc, pt-colon.h, pt-const.cc, pt-const.h, pt-decl.cc, pt-decl.h, pt-eval.cc, pt-eval.h, pt-except.cc, pt-except.h, pt-exp.cc, pt-exp.h, pt-fcn-handle.cc, pt-fcn-handle.h, pt-id.cc, pt-id.h, pt-idx.cc, pt-idx.h, pt-jump.h, pt-loop.cc, pt-loop.h, pt-mat.cc, pt-mat.h, pt-misc.cc, pt-misc.h, pt-pr-code.cc, pt-pr-code.h, pt-select.cc, pt-select.h, pt-spmd.cc, pt-spmd.h, pt-stmt.cc, pt-stmt.h, pt-tm-const.cc, pt-tm-const.h, pt-unop.cc, pt-unop.h, pt-walk.cc, pt-walk.h, pt.cc, pt.h, token.cc, token.h, Range.cc, Range.h, idx-vector.cc, idx-vector.h, range-fwd.h, CollocWt.cc, CollocWt.h, aepbalance.cc, aepbalance.h, chol.cc, chol.h, gepbalance.cc, gepbalance.h, gsvd.cc, gsvd.h, hess.cc, hess.h, lo-mappers.cc, lo-mappers.h, lo-specfun.cc, lo-specfun.h, lu.cc, lu.h, oct-convn.cc, oct-convn.h, oct-fftw.cc, oct-fftw.h, oct-norm.cc, oct-norm.h, oct-rand.cc, oct-rand.h, oct-spparms.cc, oct-spparms.h, qr.cc, qr.h, qrp.cc, qrp.h, randgamma.cc, randgamma.h, randmtzig.cc, randmtzig.h, randpoisson.cc, randpoisson.h, schur.cc, schur.h, sparse-chol.cc, sparse-chol.h, sparse-lu.cc, sparse-lu.h, sparse-qr.cc, sparse-qr.h, svd.cc, svd.h, child-list.cc, child-list.h, dir-ops.cc, dir-ops.h, file-ops.cc, file-ops.h, file-stat.cc, file-stat.h, lo-sysdep.cc, lo-sysdep.h, lo-sysinfo.cc, lo-sysinfo.h, mach-info.cc, mach-info.h, oct-env.cc, oct-env.h, oct-group.cc, oct-group.h, oct-password.cc, oct-password.h, oct-syscalls.cc, oct-syscalls.h, oct-time.cc, oct-time.h, oct-uname.cc, oct-uname.h, action-container.cc, action-container.h, base-list.h, cmd-edit.cc, cmd-edit.h, cmd-hist.cc, cmd-hist.h, f77-fcn.h, file-info.cc, file-info.h, lo-array-errwarn.cc, lo-array-errwarn.h, lo-hash.cc, lo-hash.h, lo-ieee.h, lo-regexp.cc, lo-regexp.h, lo-utils.cc, lo-utils.h, oct-base64.cc, oct-base64.h, oct-glob.cc, oct-glob.h, oct-inttypes.h, oct-mutex.cc, oct-mutex.h, oct-refcount.h, oct-shlib.cc, oct-shlib.h, oct-sparse.cc, oct-sparse.h, oct-string.h, octave-preserve-stream-state.h, pathsearch.cc, pathsearch.h, quit.cc, quit.h, unwind-prot.cc, unwind-prot.h, url-transfer.cc, url-transfer.h: Re-indent code after switch to using namespace macros.
author Rik <rik@octave.org>
date Thu, 01 Dec 2022 18:02:15 -0800
parents e88a07dec498
children 597f3ee61a48
comparison
equal deleted inserted replaced
31605:e88a07dec498 31607:aac27ad79be6
303 303
304 OCTAVE_BEGIN_NAMESPACE(octave) 304 OCTAVE_BEGIN_NAMESPACE(octave)
305 305
306 OCTAVE_BEGIN_NAMESPACE(math) 306 OCTAVE_BEGIN_NAMESPACE(math)
307 307
308 template <typename T> 308 template <typename T>
309 T 309 T
310 svd<T>::left_singular_matrix (void) const 310 svd<T>::left_singular_matrix (void) const
311 { 311 {
312 if (m_type == svd::Type::sigma_only) 312 if (m_type == svd::Type::sigma_only)
313 (*current_liboctave_error_handler) 313 (*current_liboctave_error_handler)
314 ("svd: U not computed because type == svd::sigma_only"); 314 ("svd: U not computed because type == svd::sigma_only");
315 315
316 return m_left_sm; 316 return m_left_sm;
317 } 317 }
318 318
319 template <typename T> 319 template <typename T>
320 T 320 T
321 svd<T>::right_singular_matrix (void) const 321 svd<T>::right_singular_matrix (void) const
322 { 322 {
323 if (m_type == svd::Type::sigma_only) 323 if (m_type == svd::Type::sigma_only)
324 (*current_liboctave_error_handler) 324 (*current_liboctave_error_handler)
325 ("svd: V not computed because type == svd::sigma_only"); 325 ("svd: V not computed because type == svd::sigma_only");
326 326
327 return m_right_sm; 327 return m_right_sm;
328 } 328 }
329 329
330 // GESVD specializations 330 // GESVD specializations
331 331
332 #define GESVD_REAL_STEP(f, F) \ 332 #define GESVD_REAL_STEP(f, F) \
333 F77_XFCN (f, F, (F77_CONST_CHAR_ARG2 (&jobu, 1), \ 333 F77_XFCN (f, F, (F77_CONST_CHAR_ARG2 (&jobu, 1), \
334 F77_CONST_CHAR_ARG2 (&jobv, 1), \ 334 F77_CONST_CHAR_ARG2 (&jobv, 1), \
335 m, n, tmp_data, m1, s_vec, u, m1, vt, \ 335 m, n, tmp_data, m1, s_vec, u, m1, vt, \
346 CMPLX_ARG (work.data ()), \ 346 CMPLX_ARG (work.data ()), \
347 lwork, rwork.data (), info \ 347 lwork, rwork.data (), info \
348 F77_CHAR_ARG_LEN (1) \ 348 F77_CHAR_ARG_LEN (1) \
349 F77_CHAR_ARG_LEN (1))) 349 F77_CHAR_ARG_LEN (1)))
350 350
351 // DGESVD 351 // DGESVD
352 template<> 352 template<>
353 OCTAVE_API void 353 OCTAVE_API void
354 svd<Matrix>::gesvd (char& jobu, char& jobv, F77_INT m, F77_INT n, 354 svd<Matrix>::gesvd (char& jobu, char& jobv, F77_INT m, F77_INT n,
355 double *tmp_data, F77_INT m1, double *s_vec, 355 double *tmp_data, F77_INT m1, double *s_vec,
356 double *u, double *vt, F77_INT nrow_vt1, 356 double *u, double *vt, F77_INT nrow_vt1,
357 std::vector<double>& work, F77_INT& lwork, 357 std::vector<double>& work, F77_INT& lwork,
358 F77_INT& info) 358 F77_INT& info)
359 { 359 {
360 GESVD_REAL_STEP (dgesvd, DGESVD); 360 GESVD_REAL_STEP (dgesvd, DGESVD);
361 361
362 lwork = static_cast<F77_INT> (work[0]); 362 lwork = static_cast<F77_INT> (work[0]);
363 work.reserve (lwork); 363 work.reserve (lwork);
364 364
365 GESVD_REAL_STEP (dgesvd, DGESVD); 365 GESVD_REAL_STEP (dgesvd, DGESVD);
366 } 366 }
367 367
368 // SGESVD 368 // SGESVD
369 template<> 369 template<>
370 OCTAVE_API void 370 OCTAVE_API void
371 svd<FloatMatrix>::gesvd (char& jobu, char& jobv, F77_INT m, F77_INT n, 371 svd<FloatMatrix>::gesvd (char& jobu, char& jobv, F77_INT m, F77_INT n,
372 float *tmp_data, F77_INT m1, float *s_vec, 372 float *tmp_data, F77_INT m1, float *s_vec,
373 float *u, float *vt, F77_INT nrow_vt1, 373 float *u, float *vt, F77_INT nrow_vt1,
374 std::vector<float>& work, F77_INT& lwork, 374 std::vector<float>& work, F77_INT& lwork,
375 F77_INT& info) 375 F77_INT& info)
376 { 376 {
377 GESVD_REAL_STEP (sgesvd, SGESVD); 377 GESVD_REAL_STEP (sgesvd, SGESVD);
378 378
379 lwork = static_cast<F77_INT> (work[0]); 379 lwork = static_cast<F77_INT> (work[0]);
380 work.reserve (lwork); 380 work.reserve (lwork);
381 381
382 GESVD_REAL_STEP (sgesvd, SGESVD); 382 GESVD_REAL_STEP (sgesvd, SGESVD);
383 } 383 }
384 384
385 // ZGESVD 385 // ZGESVD
386 template<> 386 template<>
387 OCTAVE_API void 387 OCTAVE_API void
388 svd<ComplexMatrix>::gesvd (char& jobu, char& jobv, F77_INT m, F77_INT n, 388 svd<ComplexMatrix>::gesvd (char& jobu, char& jobv, F77_INT m, F77_INT n,
389 Complex *tmp_data, F77_INT m1, double *s_vec, 389 Complex *tmp_data, F77_INT m1, double *s_vec,
390 Complex *u, Complex *vt, F77_INT nrow_vt1, 390 Complex *u, Complex *vt, F77_INT nrow_vt1,
391 std::vector<Complex>& work, F77_INT& lwork, 391 std::vector<Complex>& work, F77_INT& lwork,
392 F77_INT& info) 392 F77_INT& info)
393 { 393 {
394 std::vector<double> rwork (5 * std::max (m, n)); 394 std::vector<double> rwork (5 * std::max (m, n));
395 395
396 GESVD_COMPLEX_STEP (zgesvd, ZGESVD, F77_DBLE_CMPLX_ARG); 396 GESVD_COMPLEX_STEP (zgesvd, ZGESVD, F77_DBLE_CMPLX_ARG);
397 397
398 lwork = static_cast<F77_INT> (work[0].real ()); 398 lwork = static_cast<F77_INT> (work[0].real ());
399 work.reserve (lwork); 399 work.reserve (lwork);
400 400
401 GESVD_COMPLEX_STEP (zgesvd, ZGESVD, F77_DBLE_CMPLX_ARG); 401 GESVD_COMPLEX_STEP (zgesvd, ZGESVD, F77_DBLE_CMPLX_ARG);
402 } 402 }
403 403
404 // CGESVD 404 // CGESVD
405 template<> 405 template<>
406 OCTAVE_API void 406 OCTAVE_API void
407 svd<FloatComplexMatrix>::gesvd (char& jobu, char& jobv, F77_INT m, 407 svd<FloatComplexMatrix>::gesvd (char& jobu, char& jobv, F77_INT m,
408 F77_INT n, FloatComplex *tmp_data, 408 F77_INT n, FloatComplex *tmp_data,
409 F77_INT m1, float *s_vec, FloatComplex *u, 409 F77_INT m1, float *s_vec, FloatComplex *u,
410 FloatComplex *vt, F77_INT nrow_vt1, 410 FloatComplex *vt, F77_INT nrow_vt1,
411 std::vector<FloatComplex>& work, 411 std::vector<FloatComplex>& work,
412 F77_INT& lwork, F77_INT& info) 412 F77_INT& lwork, F77_INT& info)
413 { 413 {
414 std::vector<float> rwork (5 * std::max (m, n)); 414 std::vector<float> rwork (5 * std::max (m, n));
415 415
416 GESVD_COMPLEX_STEP (cgesvd, CGESVD, F77_CMPLX_ARG); 416 GESVD_COMPLEX_STEP (cgesvd, CGESVD, F77_CMPLX_ARG);
417 417
418 lwork = static_cast<F77_INT> (work[0].real ()); 418 lwork = static_cast<F77_INT> (work[0].real ());
419 work.reserve (lwork); 419 work.reserve (lwork);
420 420
421 GESVD_COMPLEX_STEP (cgesvd, CGESVD, F77_CMPLX_ARG); 421 GESVD_COMPLEX_STEP (cgesvd, CGESVD, F77_CMPLX_ARG);
422 } 422 }
423 423
424 #undef GESVD_REAL_STEP 424 #undef GESVD_REAL_STEP
425 #undef GESVD_COMPLEX_STEP 425 #undef GESVD_COMPLEX_STEP
426 426
427 // GESDD specializations 427 // GESDD specializations
428 428
429 #define GESDD_REAL_STEP(f, F) \ 429 #define GESDD_REAL_STEP(f, F) \
430 F77_XFCN (f, F, (F77_CONST_CHAR_ARG2 (&jobz, 1), \ 430 F77_XFCN (f, F, (F77_CONST_CHAR_ARG2 (&jobz, 1), \
431 m, n, tmp_data, m1, s_vec, u, m1, vt, nrow_vt1, \ 431 m, n, tmp_data, m1, s_vec, u, m1, vt, nrow_vt1, \
432 work.data (), lwork, iwork, info \ 432 work.data (), lwork, iwork, info \
439 CMPLX_ARG (vt), nrow_vt1, \ 439 CMPLX_ARG (vt), nrow_vt1, \
440 CMPLX_ARG (work.data ()), lwork, \ 440 CMPLX_ARG (work.data ()), lwork, \
441 rwork.data (), iwork, info \ 441 rwork.data (), iwork, info \
442 F77_CHAR_ARG_LEN (1))) 442 F77_CHAR_ARG_LEN (1)))
443 443
444 // DGESDD 444 // DGESDD
445 template<> 445 template<>
446 OCTAVE_API void 446 OCTAVE_API void
447 svd<Matrix>::gesdd (char& jobz, F77_INT m, F77_INT n, double *tmp_data, 447 svd<Matrix>::gesdd (char& jobz, F77_INT m, F77_INT n, double *tmp_data,
448 F77_INT m1, double *s_vec, double *u, double *vt, 448 F77_INT m1, double *s_vec, double *u, double *vt,
449 F77_INT nrow_vt1, std::vector<double>& work, 449 F77_INT nrow_vt1, std::vector<double>& work,
450 F77_INT& lwork, F77_INT *iwork, F77_INT& info) 450 F77_INT& lwork, F77_INT *iwork, F77_INT& info)
451 { 451 {
452 GESDD_REAL_STEP (dgesdd, DGESDD); 452 GESDD_REAL_STEP (dgesdd, DGESDD);
453 453
454 lwork = static_cast<F77_INT> (work[0]); 454 lwork = static_cast<F77_INT> (work[0]);
455 work.reserve (lwork); 455 work.reserve (lwork);
456 456
457 GESDD_REAL_STEP (dgesdd, DGESDD); 457 GESDD_REAL_STEP (dgesdd, DGESDD);
458 } 458 }
459 459
460 // SGESDD 460 // SGESDD
461 template<> 461 template<>
462 OCTAVE_API void 462 OCTAVE_API void
463 svd<FloatMatrix>::gesdd (char& jobz, F77_INT m, F77_INT n, float *tmp_data, 463 svd<FloatMatrix>::gesdd (char& jobz, F77_INT m, F77_INT n, float *tmp_data,
464 F77_INT m1, float *s_vec, float *u, float *vt, 464 F77_INT m1, float *s_vec, float *u, float *vt,
465 F77_INT nrow_vt1, std::vector<float>& work, 465 F77_INT nrow_vt1, std::vector<float>& work,
466 F77_INT& lwork, F77_INT *iwork, F77_INT& info) 466 F77_INT& lwork, F77_INT *iwork, F77_INT& info)
467 { 467 {
468 GESDD_REAL_STEP (sgesdd, SGESDD); 468 GESDD_REAL_STEP (sgesdd, SGESDD);
469 469
470 lwork = static_cast<F77_INT> (work[0]); 470 lwork = static_cast<F77_INT> (work[0]);
471 work.reserve (lwork); 471 work.reserve (lwork);
472 472
473 GESDD_REAL_STEP (sgesdd, SGESDD); 473 GESDD_REAL_STEP (sgesdd, SGESDD);
474 } 474 }
475 475
476 // ZGESDD 476 // ZGESDD
477 template<> 477 template<>
478 OCTAVE_API void 478 OCTAVE_API void
479 svd<ComplexMatrix>::gesdd (char& jobz, F77_INT m, F77_INT n, 479 svd<ComplexMatrix>::gesdd (char& jobz, F77_INT m, F77_INT n,
480 Complex *tmp_data, F77_INT m1, double *s_vec, 480 Complex *tmp_data, F77_INT m1, double *s_vec,
481 Complex *u, Complex *vt, F77_INT nrow_vt1, 481 Complex *u, Complex *vt, F77_INT nrow_vt1,
482 std::vector<Complex>& work, F77_INT& lwork, 482 std::vector<Complex>& work, F77_INT& lwork,
483 F77_INT *iwork, F77_INT& info) 483 F77_INT *iwork, F77_INT& info)
484 { 484 {
485 485
486 F77_INT min_mn = std::min (m, n); 486 F77_INT min_mn = std::min (m, n);
487 F77_INT max_mn = std::max (m, n); 487 F77_INT max_mn = std::max (m, n);
488 488
489 F77_INT lrwork; 489 F77_INT lrwork;
490 if (jobz == 'N') 490 if (jobz == 'N')
491 lrwork = 7*min_mn; 491 lrwork = 7*min_mn;
492 else 492 else
493 lrwork = min_mn * std::max (5*min_mn+5, 2*max_mn+2*min_mn+1); 493 lrwork = min_mn * std::max (5*min_mn+5, 2*max_mn+2*min_mn+1);
494 494
495 std::vector<double> rwork (lrwork); 495 std::vector<double> rwork (lrwork);
496 496
497 GESDD_COMPLEX_STEP (zgesdd, ZGESDD, F77_DBLE_CMPLX_ARG); 497 GESDD_COMPLEX_STEP (zgesdd, ZGESDD, F77_DBLE_CMPLX_ARG);
498 498
499 lwork = static_cast<F77_INT> (work[0].real ()); 499 lwork = static_cast<F77_INT> (work[0].real ());
500 work.reserve (lwork); 500 work.reserve (lwork);
501 501
502 GESDD_COMPLEX_STEP (zgesdd, ZGESDD, F77_DBLE_CMPLX_ARG); 502 GESDD_COMPLEX_STEP (zgesdd, ZGESDD, F77_DBLE_CMPLX_ARG);
503 } 503 }
504 504
505 // CGESDD 505 // CGESDD
506 template<> 506 template<>
507 OCTAVE_API void 507 OCTAVE_API void
508 svd<FloatComplexMatrix>::gesdd (char& jobz, F77_INT m, F77_INT n, 508 svd<FloatComplexMatrix>::gesdd (char& jobz, F77_INT m, F77_INT n,
509 FloatComplex *tmp_data, F77_INT m1, 509 FloatComplex *tmp_data, F77_INT m1,
510 float *s_vec, FloatComplex *u, 510 float *s_vec, FloatComplex *u,
511 FloatComplex *vt, F77_INT nrow_vt1, 511 FloatComplex *vt, F77_INT nrow_vt1,
512 std::vector<FloatComplex>& work, 512 std::vector<FloatComplex>& work,
513 F77_INT& lwork, F77_INT *iwork, 513 F77_INT& lwork, F77_INT *iwork,
514 F77_INT& info) 514 F77_INT& info)
515 { 515 {
516 F77_INT min_mn = std::min (m, n); 516 F77_INT min_mn = std::min (m, n);
517 F77_INT max_mn = std::max (m, n); 517 F77_INT max_mn = std::max (m, n);
518 518
519 F77_INT lrwork; 519 F77_INT lrwork;
520 if (jobz == 'N') 520 if (jobz == 'N')
521 lrwork = 7*min_mn; 521 lrwork = 7*min_mn;
522 else 522 else
523 lrwork = min_mn * std::max (5*min_mn+5, 2*max_mn+2*min_mn+1); 523 lrwork = min_mn * std::max (5*min_mn+5, 2*max_mn+2*min_mn+1);
524 std::vector<float> rwork (lrwork); 524 std::vector<float> rwork (lrwork);
525 525
526 GESDD_COMPLEX_STEP (cgesdd, CGESDD, F77_CMPLX_ARG); 526 GESDD_COMPLEX_STEP (cgesdd, CGESDD, F77_CMPLX_ARG);
527 527
528 lwork = static_cast<F77_INT> (work[0].real ()); 528 lwork = static_cast<F77_INT> (work[0].real ());
529 work.reserve (lwork); 529 work.reserve (lwork);
530 530
531 GESDD_COMPLEX_STEP (cgesdd, CGESDD, F77_CMPLX_ARG); 531 GESDD_COMPLEX_STEP (cgesdd, CGESDD, F77_CMPLX_ARG);
532 } 532 }
533 533
534 #undef GESDD_REAL_STEP 534 #undef GESDD_REAL_STEP
535 #undef GESDD_COMPLEX_STEP 535 #undef GESDD_COMPLEX_STEP
536 536
537 // GEJSV specializations 537 // GEJSV specializations
538 538
539 #define GEJSV_REAL_STEP(f, F) \ 539 #define GEJSV_REAL_STEP(f, F) \
540 F77_XFCN (f, F, (F77_CONST_CHAR_ARG2 (&joba, 1), \ 540 F77_XFCN (f, F, (F77_CONST_CHAR_ARG2 (&joba, 1), \
541 F77_CONST_CHAR_ARG2 (&jobu, 1), \ 541 F77_CONST_CHAR_ARG2 (&jobu, 1), \
542 F77_CONST_CHAR_ARG2 (&jobv, 1), \ 542 F77_CONST_CHAR_ARG2 (&jobv, 1), \
569 F77_CHAR_ARG_LEN (1) \ 569 F77_CHAR_ARG_LEN (1) \
570 F77_CHAR_ARG_LEN (1) \ 570 F77_CHAR_ARG_LEN (1) \
571 F77_CHAR_ARG_LEN (1) \ 571 F77_CHAR_ARG_LEN (1) \
572 F77_CHAR_ARG_LEN (1))) 572 F77_CHAR_ARG_LEN (1)))
573 573
574 // DGEJSV 574 // DGEJSV
575 template<> 575 template<>
576 void 576 void
577 svd<Matrix>::gejsv (char& joba, char& jobu, char& jobv, 577 svd<Matrix>::gejsv (char& joba, char& jobu, char& jobv,
578 char& jobr, char& jobt, char& jobp, 578 char& jobr, char& jobt, char& jobp,
579 F77_INT m, F77_INT n, 579 F77_INT m, F77_INT n,
580 P *tmp_data, F77_INT m1, DM_P *s_vec, P *u, 580 P *tmp_data, F77_INT m1, DM_P *s_vec, P *u,
581 P *v, F77_INT nrow_v1, std::vector<P>& work, 581 P *v, F77_INT nrow_v1, std::vector<P>& work,
582 F77_INT& lwork, std::vector<F77_INT>& iwork, 582 F77_INT& lwork, std::vector<F77_INT>& iwork,
583 F77_INT& info) 583 F77_INT& info)
584 { 584 {
585 lwork = gejsv_lwork<Matrix>::optimal (joba, jobu, jobv, m, n); 585 lwork = gejsv_lwork<Matrix>::optimal (joba, jobu, jobv, m, n);
586 work.reserve (lwork); 586 work.reserve (lwork);
587 587
588 GEJSV_REAL_STEP (dgejsv, DGEJSV); 588 GEJSV_REAL_STEP (dgejsv, DGEJSV);
589 } 589 }
590 590
591 // SGEJSV 591 // SGEJSV
592 template<> 592 template<>
593 void 593 void
594 svd<FloatMatrix>::gejsv (char& joba, char& jobu, char& jobv, 594 svd<FloatMatrix>::gejsv (char& joba, char& jobu, char& jobv,
595 char& jobr, char& jobt, char& jobp, 595 char& jobr, char& jobt, char& jobp,
596 F77_INT m, F77_INT n, 596 F77_INT m, F77_INT n,
597 P *tmp_data, F77_INT m1, DM_P *s_vec, P *u, 597 P *tmp_data, F77_INT m1, DM_P *s_vec, P *u,
598 P *v, F77_INT nrow_v1, std::vector<P>& work, 598 P *v, F77_INT nrow_v1, std::vector<P>& work,
599 F77_INT& lwork, std::vector<F77_INT>& iwork, 599 F77_INT& lwork, std::vector<F77_INT>& iwork,
600 F77_INT& info) 600 F77_INT& info)
601 { 601 {
602 lwork = gejsv_lwork<FloatMatrix>::optimal (joba, jobu, jobv, m, n); 602 lwork = gejsv_lwork<FloatMatrix>::optimal (joba, jobu, jobv, m, n);
603 work.reserve (lwork); 603 work.reserve (lwork);
604 604
605 GEJSV_REAL_STEP (sgejsv, SGEJSV); 605 GEJSV_REAL_STEP (sgejsv, SGEJSV);
606 } 606 }
607 607
608 // ZGEJSV 608 // ZGEJSV
609 template<> 609 template<>
610 void 610 void
611 svd<ComplexMatrix>::gejsv (char& joba, char& jobu, char& jobv, 611 svd<ComplexMatrix>::gejsv (char& joba, char& jobu, char& jobv,
612 char& jobr, char& jobt, char& jobp, 612 char& jobr, char& jobt, char& jobp,
613 F77_INT m, F77_INT n, 613 F77_INT m, F77_INT n,
614 P *tmp_data, F77_INT m1, DM_P *s_vec, P *u, 614 P *tmp_data, F77_INT m1, DM_P *s_vec, P *u,
615 P *v, F77_INT nrow_v1, std::vector<P>& work, 615 P *v, F77_INT nrow_v1, std::vector<P>& work,
616 F77_INT& lwork, std::vector<F77_INT>& iwork, 616 F77_INT& lwork, std::vector<F77_INT>& iwork,
617 F77_INT& info) 617 F77_INT& info)
618 { 618 {
619 F77_INT lrwork = -1; // work space size query 619 F77_INT lrwork = -1; // work space size query
620 std::vector<double> rwork (1); 620 std::vector<double> rwork (1);
621 work.reserve (2); 621 work.reserve (2);
622 622
623 GEJSV_COMPLEX_STEP (zgejsv, ZGEJSV, F77_DBLE_CMPLX_ARG); 623 GEJSV_COMPLEX_STEP (zgejsv, ZGEJSV, F77_DBLE_CMPLX_ARG);
624 624
625 lwork = static_cast<F77_INT> (work[0].real ()); 625 lwork = static_cast<F77_INT> (work[0].real ());
626 work.reserve (lwork); 626 work.reserve (lwork);
627 627
628 lrwork = static_cast<F77_INT> (rwork[0]); 628 lrwork = static_cast<F77_INT> (rwork[0]);
629 rwork.reserve (lrwork); 629 rwork.reserve (lrwork);
630 630
631 F77_INT liwork = static_cast<F77_INT> (iwork[0]); 631 F77_INT liwork = static_cast<F77_INT> (iwork[0]);
632 iwork.reserve (liwork); 632 iwork.reserve (liwork);
633 633
634 GEJSV_COMPLEX_STEP (zgejsv, ZGEJSV, F77_DBLE_CMPLX_ARG); 634 GEJSV_COMPLEX_STEP (zgejsv, ZGEJSV, F77_DBLE_CMPLX_ARG);
635 } 635 }
636 636
637 // CGEJSV 637 // CGEJSV
638 template<> 638 template<>
639 void 639 void
640 svd<FloatComplexMatrix>::gejsv (char& joba, char& jobu, char& jobv, 640 svd<FloatComplexMatrix>::gejsv (char& joba, char& jobu, char& jobv,
641 char& jobr, char& jobt, char& jobp, 641 char& jobr, char& jobt, char& jobp,
642 F77_INT m, F77_INT n, P *tmp_data, 642 F77_INT m, F77_INT n, P *tmp_data,
643 F77_INT m1, DM_P *s_vec, P *u, P *v, 643 F77_INT m1, DM_P *s_vec, P *u, P *v,
644 F77_INT nrow_v1, std::vector<P>& work, 644 F77_INT nrow_v1, std::vector<P>& work,
645 F77_INT& lwork, 645 F77_INT& lwork,
646 std::vector<F77_INT>& iwork, F77_INT& info) 646 std::vector<F77_INT>& iwork, F77_INT& info)
647 { 647 {
648 F77_INT lrwork = -1; // work space size query 648 F77_INT lrwork = -1; // work space size query
649 std::vector<float> rwork (1); 649 std::vector<float> rwork (1);
650 work.reserve (2); 650 work.reserve (2);
651 651
652 GEJSV_COMPLEX_STEP (cgejsv, CGEJSV, F77_CMPLX_ARG); 652 GEJSV_COMPLEX_STEP (cgejsv, CGEJSV, F77_CMPLX_ARG);
653 653
654 lwork = static_cast<F77_INT> (work[0].real ()); 654 lwork = static_cast<F77_INT> (work[0].real ());
655 work.reserve (lwork); 655 work.reserve (lwork);
656 656
657 lrwork = static_cast<F77_INT> (rwork[0]); 657 lrwork = static_cast<F77_INT> (rwork[0]);
658 rwork.reserve (lrwork); 658 rwork.reserve (lrwork);
659 659
660 F77_INT liwork = static_cast<F77_INT> (iwork[0]); 660 F77_INT liwork = static_cast<F77_INT> (iwork[0]);
661 iwork.reserve (liwork); 661 iwork.reserve (liwork);
662 662
663 GEJSV_COMPLEX_STEP (cgejsv, CGEJSV, F77_CMPLX_ARG); 663 GEJSV_COMPLEX_STEP (cgejsv, CGEJSV, F77_CMPLX_ARG);
664 } 664 }
665 665
666 #undef GEJSV_REAL_STEP 666 #undef GEJSV_REAL_STEP
667 #undef GEJSV_COMPLEX_STEP 667 #undef GEJSV_COMPLEX_STEP
668 668
669 template<typename T> 669 template<typename T>
670 svd<T>::svd (const T& a, svd::Type type, svd::Driver driver) 670 svd<T>::svd (const T& a, svd::Type type, svd::Driver driver)
671 : m_type (type), m_driver (driver), m_left_sm (), m_sigma (), 671 : m_type (type), m_driver (driver), m_left_sm (), m_sigma (),
672 m_right_sm () 672 m_right_sm ()
673 { 673 {
674 F77_INT info; 674 F77_INT info;
675 675
676 F77_INT m = to_f77_int (a.rows ()); 676 F77_INT m = to_f77_int (a.rows ());
677 F77_INT n = to_f77_int (a.cols ()); 677 F77_INT n = to_f77_int (a.cols ());
678 678
679 if (m == 0 || n == 0) 679 if (m == 0 || n == 0)
680 { 680 {
681 switch (m_type)
682 {
683 case svd::Type::std:
684 m_left_sm = T (m, m, 0);
685 for (F77_INT i = 0; i < m; i++)
686 m_left_sm.xelem (i, i) = 1;
687 m_sigma = DM_T (m, n);
688 m_right_sm = T (n, n, 0);
689 for (F77_INT i = 0; i < n; i++)
690 m_right_sm.xelem (i, i) = 1;
691 break;
692
693 case svd::Type::economy:
694 m_left_sm = T (m, 0, 0);
695 m_sigma = DM_T (0, 0);
696 m_right_sm = T (n, 0, 0);
697 break;
698
699 case svd::Type::sigma_only:
700 default:
701 m_sigma = DM_T (0, 1);
702 break;
703 }
704 return;
705 }
706
707 T atmp = a;
708 P *tmp_data = atmp.fortran_vec ();
709
710 F77_INT min_mn = (m < n ? m : n);
711
712 char jobu = 'A';
713 char jobv = 'A';
714
715 F77_INT ncol_u = m;
716 F77_INT nrow_vt = n;
717 F77_INT nrow_s = m;
718 F77_INT ncol_s = n;
719
720 switch (m_type) 681 switch (m_type)
721 { 682 {
683 case svd::Type::std:
684 m_left_sm = T (m, m, 0);
685 for (F77_INT i = 0; i < m; i++)
686 m_left_sm.xelem (i, i) = 1;
687 m_sigma = DM_T (m, n);
688 m_right_sm = T (n, n, 0);
689 for (F77_INT i = 0; i < n; i++)
690 m_right_sm.xelem (i, i) = 1;
691 break;
692
722 case svd::Type::economy: 693 case svd::Type::economy:
723 jobu = jobv = 'S'; 694 m_left_sm = T (m, 0, 0);
724 ncol_u = nrow_vt = nrow_s = ncol_s = min_mn; 695 m_sigma = DM_T (0, 0);
696 m_right_sm = T (n, 0, 0);
725 break; 697 break;
726 698
727 case svd::Type::sigma_only: 699 case svd::Type::sigma_only:
728
729 // Note: for this case, both jobu and jobv should be 'N', but there
730 // seems to be a bug in dgesvd from Lapack V2.0. To demonstrate the
731 // bug, set both jobu and jobv to 'N' and find the singular values of
732 // [eye(3), eye(3)]. The result is [-sqrt(2), -sqrt(2), -sqrt(2)].
733 //
734 // For Lapack 3.0, this problem seems to be fixed.
735
736 jobu = jobv = 'N';
737 ncol_u = nrow_vt = 1;
738 break;
739
740 default: 700 default:
701 m_sigma = DM_T (0, 1);
741 break; 702 break;
742 } 703 }
743 704 return;
744 if (! (jobu == 'N' || jobu == 'O')) 705 }
745 m_left_sm.resize (m, ncol_u); 706
746 707 T atmp = a;
747 P *u = m_left_sm.fortran_vec (); 708 P *tmp_data = atmp.fortran_vec ();
748 709
749 m_sigma.resize (nrow_s, ncol_s); 710 F77_INT min_mn = (m < n ? m : n);
750 DM_P *s_vec = m_sigma.fortran_vec (); 711
751 712 char jobu = 'A';
752 if (! (jobv == 'N' || jobv == 'O')) 713 char jobv = 'A';
714
715 F77_INT ncol_u = m;
716 F77_INT nrow_vt = n;
717 F77_INT nrow_s = m;
718 F77_INT ncol_s = n;
719
720 switch (m_type)
721 {
722 case svd::Type::economy:
723 jobu = jobv = 'S';
724 ncol_u = nrow_vt = nrow_s = ncol_s = min_mn;
725 break;
726
727 case svd::Type::sigma_only:
728
729 // Note: for this case, both jobu and jobv should be 'N', but there
730 // seems to be a bug in dgesvd from Lapack V2.0. To demonstrate the
731 // bug, set both jobu and jobv to 'N' and find the singular values of
732 // [eye(3), eye(3)]. The result is [-sqrt(2), -sqrt(2), -sqrt(2)].
733 //
734 // For Lapack 3.0, this problem seems to be fixed.
735
736 jobu = jobv = 'N';
737 ncol_u = nrow_vt = 1;
738 break;
739
740 default:
741 break;
742 }
743
744 if (! (jobu == 'N' || jobu == 'O'))
745 m_left_sm.resize (m, ncol_u);
746
747 P *u = m_left_sm.fortran_vec ();
748
749 m_sigma.resize (nrow_s, ncol_s);
750 DM_P *s_vec = m_sigma.fortran_vec ();
751
752 if (! (jobv == 'N' || jobv == 'O'))
753 {
754 if (m_driver == svd::Driver::GEJSV)
755 m_right_sm.resize (n, nrow_vt);
756 else
757 m_right_sm.resize (nrow_vt, n);
758 }
759
760 P *vt = m_right_sm.fortran_vec ();
761
762 // Query _GESVD for the correct dimension of WORK.
763
764 F77_INT lwork = -1;
765
766 std::vector<P> work (1);
767
768 const F77_INT f77_int_one = static_cast<F77_INT> (1);
769 F77_INT m1 = std::max (m, f77_int_one);
770 F77_INT nrow_vt1 = std::max (nrow_vt, f77_int_one);
771
772 if (m_driver == svd::Driver::GESVD)
773 gesvd (jobu, jobv, m, n, tmp_data, m1, s_vec, u, vt, nrow_vt1,
774 work, lwork, info);
775 else if (m_driver == svd::Driver::GESDD)
776 {
777 assert (jobu == jobv);
778 char jobz = jobu;
779
780 std::vector<F77_INT> iwork (8 * std::min (m, n));
781
782 gesdd (jobz, m, n, tmp_data, m1, s_vec, u, vt, nrow_vt1,
783 work, lwork, iwork.data (), info);
784 }
785 else if (m_driver == svd::Driver::GEJSV)
786 {
787 bool transposed = false;
788 if (n > m)
753 { 789 {
754 if (m_driver == svd::Driver::GEJSV) 790 // GEJSV only accepts m >= n, thus we need to transpose here
755 m_right_sm.resize (n, nrow_vt); 791 transposed = true;
756 else 792
757 m_right_sm.resize (nrow_vt, n); 793 std::swap (m, n);
794 m1 = std::max (m, f77_int_one);
795 nrow_vt1 = std::max (n, f77_int_one); // we have m > n
796 if (m_type == svd::Type::sigma_only)
797 nrow_vt1 = 1;
798 std::swap (jobu, jobv);
799
800 atmp = atmp.hermitian ();
801 tmp_data = atmp.fortran_vec ();
802
803 // Swap pointers of U and V.
804 u = m_right_sm.fortran_vec ();
805 vt = m_left_sm.fortran_vec ();
758 } 806 }
759 807
760 P *vt = m_right_sm.fortran_vec (); 808 // translate jobu and jobv from gesvd to gejsv.
761 809 std::unordered_map<char, std::string> job_svd2jsv;
762 // Query _GESVD for the correct dimension of WORK. 810 job_svd2jsv['A'] = "FJ";
763 811 job_svd2jsv['S'] = "UV";
764 F77_INT lwork = -1; 812 job_svd2jsv['O'] = "WW";
765 813 job_svd2jsv['N'] = "NN";
766 std::vector<P> work (1); 814 jobu = job_svd2jsv[jobu][0];
767 815 jobv = job_svd2jsv[jobv][1];
768 const F77_INT f77_int_one = static_cast<F77_INT> (1); 816
769 F77_INT m1 = std::max (m, f77_int_one); 817 char joba = 'F'; // 'F': most conservative
770 F77_INT nrow_vt1 = std::max (nrow_vt, f77_int_one); 818 char jobr = 'R'; // 'R' is recommended.
771 819 char jobt = 'N'; // or 'T', but that requires U and V appear together
772 if (m_driver == svd::Driver::GESVD) 820 char jobp = 'N'; // use 'P' if denormal is poorly implemented.
773 gesvd (jobu, jobv, m, n, tmp_data, m1, s_vec, u, vt, nrow_vt1, 821
774 work, lwork, info); 822 std::vector<F77_INT> iwork (std::max<F77_INT> (m + 3*n, 1));
775 else if (m_driver == svd::Driver::GESDD) 823
776 { 824 gejsv (joba, jobu, jobv, jobr, jobt, jobp, m, n, tmp_data, m1,
777 assert (jobu == jobv); 825 s_vec, u, vt, nrow_vt1, work, lwork, iwork, info);
778 char jobz = jobu; 826
779 827 if (iwork[2] == 1)
780 std::vector<F77_INT> iwork (8 * std::min (m, n)); 828 (*current_liboctave_warning_with_id_handler)
781 829 ("Octave:convergence", "svd: (driver: GEJSV) "
782 gesdd (jobz, m, n, tmp_data, m1, s_vec, u, vt, nrow_vt1, 830 "Denormal occurred, possible loss of accuracy.");
783 work, lwork, iwork.data (), info); 831
784 } 832 if (info < 0)
785 else if (m_driver == svd::Driver::GEJSV) 833 (*current_liboctave_error_handler)
786 { 834 ("svd: (driver: GEJSV) Illegal argument at #%d",
787 bool transposed = false; 835 static_cast<int> (-info));
788 if (n > m) 836 else if (info > 0)
789 { 837 (*current_liboctave_warning_with_id_handler)
790 // GEJSV only accepts m >= n, thus we need to transpose here 838 ("Octave:convergence", "svd: (driver: GEJSV) "
791 transposed = true; 839 "Fail to converge within max sweeps, "
792 840 "possible inaccurate result.");
793 std::swap (m, n); 841
794 m1 = std::max (m, f77_int_one); 842 if (transposed) // put things that need to transpose back here
795 nrow_vt1 = std::max (n, f77_int_one); // we have m > n 843 std::swap (m, n);
796 if (m_type == svd::Type::sigma_only) 844 }
797 nrow_vt1 = 1; 845 else
798 std::swap (jobu, jobv); 846 (*current_liboctave_error_handler) ("svd: unknown driver");
799 847
800 atmp = atmp.hermitian (); 848 // LAPACK can return -0 which is a small problem (bug #55710).
801 tmp_data = atmp.fortran_vec (); 849 for (octave_idx_type i = 0; i < m_sigma.diag_length (); i++)
802 850 {
803 // Swap pointers of U and V. 851 if (! m_sigma.dgxelem (i))
804 u = m_right_sm.fortran_vec (); 852 m_sigma.dgxelem (i) = DM_P (0);
805 vt = m_left_sm.fortran_vec (); 853 }
806 } 854
807 855 // GESVD and GESDD return VT instead of V, GEJSV return V.
808 // translate jobu and jobv from gesvd to gejsv. 856 if (! (jobv == 'N' || jobv == 'O') && (m_driver != svd::Driver::GEJSV))
809 std::unordered_map<char, std::string> job_svd2jsv; 857 m_right_sm = m_right_sm.hermitian ();
810 job_svd2jsv['A'] = "FJ"; 858 }
811 job_svd2jsv['S'] = "UV"; 859
812 job_svd2jsv['O'] = "WW"; 860 // Instantiations we need.
813 job_svd2jsv['N'] = "NN"; 861
814 jobu = job_svd2jsv[jobu][0]; 862 template class svd<Matrix>;
815 jobv = job_svd2jsv[jobv][1]; 863
816 864 template class svd<FloatMatrix>;
817 char joba = 'F'; // 'F': most conservative 865
818 char jobr = 'R'; // 'R' is recommended. 866 template class svd<ComplexMatrix>;
819 char jobt = 'N'; // or 'T', but that requires U and V appear together 867
820 char jobp = 'N'; // use 'P' if denormal is poorly implemented. 868 template class svd<FloatComplexMatrix>;
821
822 std::vector<F77_INT> iwork (std::max<F77_INT> (m + 3*n, 1));
823
824 gejsv (joba, jobu, jobv, jobr, jobt, jobp, m, n, tmp_data, m1,
825 s_vec, u, vt, nrow_vt1, work, lwork, iwork, info);
826
827 if (iwork[2] == 1)
828 (*current_liboctave_warning_with_id_handler)
829 ("Octave:convergence", "svd: (driver: GEJSV) "
830 "Denormal occurred, possible loss of accuracy.");
831
832 if (info < 0)
833 (*current_liboctave_error_handler)
834 ("svd: (driver: GEJSV) Illegal argument at #%d",
835 static_cast<int> (-info));
836 else if (info > 0)
837 (*current_liboctave_warning_with_id_handler)
838 ("Octave:convergence", "svd: (driver: GEJSV) "
839 "Fail to converge within max sweeps, "
840 "possible inaccurate result.");
841
842 if (transposed) // put things that need to transpose back here
843 std::swap (m, n);
844 }
845 else
846 (*current_liboctave_error_handler) ("svd: unknown driver");
847
848 // LAPACK can return -0 which is a small problem (bug #55710).
849 for (octave_idx_type i = 0; i < m_sigma.diag_length (); i++)
850 {
851 if (! m_sigma.dgxelem (i))
852 m_sigma.dgxelem (i) = DM_P (0);
853 }
854
855 // GESVD and GESDD return VT instead of V, GEJSV return V.
856 if (! (jobv == 'N' || jobv == 'O') && (m_driver != svd::Driver::GEJSV))
857 m_right_sm = m_right_sm.hermitian ();
858 }
859
860 // Instantiations we need.
861
862 template class svd<Matrix>;
863
864 template class svd<FloatMatrix>;
865
866 template class svd<ComplexMatrix>;
867
868 template class svd<FloatComplexMatrix>;
869 869
870 OCTAVE_END_NAMESPACE(math) 870 OCTAVE_END_NAMESPACE(math)
871 OCTAVE_END_NAMESPACE(octave) 871 OCTAVE_END_NAMESPACE(octave)