Mercurial > octave

annotate liboctave/numeric/svd.cc @ 21724:aba2e6293dd8

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
use "#if ..." consistently instead of "#ifdef" and "#ifndef" * configure.ac, Doxyfile.in, make_int.cc, Backend.cc, BaseControl.cc, ButtonControl.cc, Canvas.cc, CheckBoxControl.cc, Container.cc, ContextMenu.cc, EditControl.cc, Figure.cc, FigureWindow.cc, GLCanvas.cc, KeyMap.cc, ListBoxControl.cc, Logger.cc, Menu.cc, MouseModeActionGroup.cc, Object.cc, ObjectFactory.cc, ObjectProxy.cc, Panel.cc, PopupMenuControl.cc, PushButtonControl.cc, PushTool.cc, QtHandlesUtils.cc, RadioButtonControl.cc, SliderControl.cc, TextControl.cc, TextEdit.cc, ToggleButtonControl.cc, ToggleTool.cc, ToolBar.cc, __init_qt__.cc, annotation-dialog.cc, gl-select.cc, module.mk, color-picker.cc, dialog.cc, documentation-dock-widget.cc, files-dock-widget.cc, find-files-dialog.cc, find-files-model.cc, history-dock-widget.cc, liboctgui-build-info.in.cc, file-editor-tab.cc, file-editor.cc, find-dialog.cc, marker.cc, marker.h, octave-qscintilla.cc, octave-txt-lexer.cc, main-window.cc, octave-cmd.cc, octave-dock-widget.cc, octave-gui.cc, octave-interpreter.cc, octave-qt-link.cc, parser.cc, webinfo.cc, resource-manager.cc, settings-dialog.cc, shortcut-manager.cc, terminal-dock-widget.cc, thread-manager.cc, welcome-wizard.cc, workspace-model.cc, workspace-view.cc, build-env.in.cc, Cell.cc, __contourc__.cc, __dispatch__.cc, __dsearchn__.cc, __ichol__.cc, __ilu__.cc, __lin_interpn__.cc, __luinc__.cc, __magick_read__.cc, __pchip_deriv__.cc, __qp__.cc, balance.cc, besselj.cc, betainc.cc, bitfcns.cc, bsxfun.cc, c-file-ptr-stream.cc, c-file-ptr-stream.h, cdisplay.c, cdisplay.h, cellfun.cc, coct-hdf5-types.c, colloc.cc, comment-list.cc, conv2.cc, daspk.cc, dasrt.cc, dassl.cc, data.cc, debug.cc, defaults.cc, defaults.in.h, defun.cc, det.cc, dirfns.cc, display.cc, dlmread.cc, dot.cc, dynamic-ld.cc, eig.cc, ellipj.cc, error.cc, errwarn.cc, event-queue.cc, fft.cc, fft2.cc, fftn.cc, file-io.cc, filter.cc, find.cc, ft-text-renderer.cc, gammainc.cc, gcd.cc, getgrent.cc, getpwent.cc, getrusage.cc, givens.cc, gl-render.cc, gl2ps-print.cc, graphics.cc, gripes.cc, hash.cc, help.cc, hess.cc, hex2num.cc, input.cc, inv.cc, jit-ir.cc, jit-ir.h, jit-typeinfo.cc, jit-typeinfo.h, jit-util.cc, jit-util.h, kron.cc, load-path.cc, load-save.cc, lookup.cc, ls-ascii-helper.cc, ls-hdf5.cc, ls-mat-ascii.cc, ls-mat4.cc, ls-mat5.cc, ls-oct-binary.cc, ls-oct-text.cc, ls-oct-text.h, ls-utils.cc, lsode.cc, lu.cc, mappers.cc, matrix_type.cc, max.cc, mex.cc, mex.h, mexproto.h, mgorth.cc, nproc.cc, oct-errno.in.cc, oct-fstrm.cc, oct-hdf5-types.cc, oct-hdf5-types.h, oct-hist.cc, oct-iostrm.cc, oct-lvalue.cc, oct-map.cc, oct-prcstrm.cc, oct-procbuf.cc, oct-stdstrm.h, oct-stream.cc, oct-strstrm.cc, oct-tex-lexer.in.ll, oct-tex-parser.in.yy, octave-link.cc, ordschur.cc, pager.cc, pinv.cc, pr-output.cc, procstream.cc, profiler.cc, psi.cc, pt-jit.cc, pt-jit.h, quad.cc, quadcc.cc, qz.cc, rand.cc, rcond.cc, regexp.cc, schur.cc, sighandlers.cc, sighandlers.h, siglist.c, siglist.h, sparse-xdiv.cc, sparse-xpow.cc, sparse.cc, spparms.cc, sqrtm.cc, str2double.cc, strfind.cc, strfns.cc, sub2ind.cc, svd.cc, sylvester.cc, symtab.cc, syscalls.cc, sysdep.cc, text-renderer.cc, time.cc, toplev.cc, tril.cc, tsearch.cc, txt-eng.cc, typecast.cc, urlwrite.cc, utils.cc, variables.cc, xdiv.cc, xgl2ps.c, xnorm.cc, xpow.cc, zfstream.cc, zfstream.h, __delaunayn__.cc, __eigs__.cc, __fltk_uigetfile__.cc, __glpk__.cc, __init_fltk__.cc, __init_gnuplot__.cc, __osmesa_print__.cc, __voronoi__.cc, amd.cc, audiodevinfo.cc, audioread.cc, ccolamd.cc, chol.cc, colamd.cc, convhulln.cc, dmperm.cc, fftw.cc, qr.cc, symbfact.cc, symrcm.cc, liboctinterp-build-info.in.cc, mkbuiltins, mkops, ov-base.cc, ov-bool-mat.cc, ov-bool-sparse.cc, ov-bool.cc, ov-builtin.cc, ov-cell.cc, ov-ch-mat.cc, ov-class.cc, ov-classdef.cc, ov-colon.cc, ov-complex.cc, ov-cs-list.cc, ov-cx-diag.cc, ov-cx-mat.cc, ov-cx-sparse.cc, ov-dld-fcn.cc, ov-fcn-handle.cc, ov-fcn-inline.cc, ov-fcn.cc, 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-int16.cc, ov-int32.cc, ov-int64.cc, ov-int8.cc, ov-java.cc, ov-lazy-idx.cc, ov-mex-fcn.cc, ov-null-mat.cc, ov-oncleanup.cc, ov-perm.cc, ov-range.cc, ov-re-diag.cc, ov-re-mat.cc, ov-re-sparse.cc, ov-scalar.cc, ov-str-mat.cc, ov-struct.cc, ov-typeinfo.cc, ov-uint16.cc, ov-uint32.cc, ov-uint64.cc, ov-uint8.cc, ov-usr-fcn.cc, ov-usr-fcn.h, ov.cc, ovl.cc, octave.cc, octave.h, op-b-b.cc, op-b-bm.cc, op-b-sbm.cc, op-bm-b.cc, op-bm-bm.cc, op-bm-sbm.cc, op-cdm-cdm.cc, op-cdm-cm.cc, op-cdm-dm.cc, op-cdm-m.cc, op-cell.cc, op-chm.cc, op-class.cc, op-cm-cdm.cc, op-cm-cm.cc, op-cm-cs.cc, op-cm-dm.cc, op-cm-m.cc, op-cm-s.cc, op-cm-scm.cc, op-cm-sm.cc, op-cs-cm.cc, op-cs-cs.cc, op-cs-m.cc, op-cs-s.cc, op-cs-scm.cc, op-cs-sm.cc, op-dm-cdm.cc, op-dm-cm.cc, op-dm-dm.cc, op-dm-m.cc, op-dm-scm.cc, op-dm-sm.cc, op-dm-template.cc, op-dms-template.cc, op-double-conv.cc, op-fcdm-fcdm.cc, op-fcdm-fcm.cc, op-fcdm-fdm.cc, op-fcdm-fm.cc, op-fcm-fcdm.cc, op-fcm-fcm.cc, op-fcm-fcs.cc, op-fcm-fdm.cc, op-fcm-fm.cc, op-fcm-fs.cc, op-fcn.cc, op-fcs-fcm.cc, op-fcs-fcs.cc, op-fcs-fm.cc, op-fcs-fs.cc, op-fdm-fcdm.cc, op-fdm-fcm.cc, op-fdm-fdm.cc, op-fdm-fm.cc, op-float-conv.cc, op-fm-fcdm.cc, op-fm-fcm.cc, op-fm-fcs.cc, op-fm-fdm.cc, op-fm-fm.cc, op-fm-fs.cc, op-fs-fcm.cc, op-fs-fcs.cc, op-fs-fm.cc, op-fs-fs.cc, op-i16-i16.cc, op-i32-i32.cc, op-i64-i64.cc, op-i8-i8.cc, op-int-concat.cc, op-int-conv.cc, op-m-cdm.cc, op-m-cm.cc, op-m-cs.cc, op-m-dm.cc, op-m-m.cc, op-m-s.cc, op-m-scm.cc, op-m-sm.cc, op-pm-pm.cc, op-pm-scm.cc, op-pm-sm.cc, op-pm-template.cc, op-range.cc, op-s-cm.cc, op-s-cs.cc, op-s-m.cc, op-s-s.cc, op-s-scm.cc, op-s-sm.cc, op-sbm-b.cc, op-sbm-bm.cc, op-sbm-sbm.cc, op-scm-cm.cc, op-scm-cs.cc, op-scm-m.cc, op-scm-s.cc, op-scm-scm.cc, op-scm-sm.cc, op-sm-cm.cc, op-sm-cs.cc, op-sm-m.cc, op-sm-s.cc, op-sm-scm.cc, op-sm-sm.cc, op-str-m.cc, op-str-s.cc, op-str-str.cc, op-struct.cc, op-ui16-ui16.cc, op-ui32-ui32.cc, op-ui64-ui64.cc, op-ui8-ui8.cc, lex.ll, oct-parse.in.yy, pt-arg-list.cc, pt-array-list.cc, pt-assign.cc, pt-binop.cc, pt-bp.cc, pt-cbinop.cc, pt-cell.cc, pt-check.cc, pt-classdef.cc, pt-cmd.cc, pt-colon.cc, pt-const.cc, pt-decl.cc, pt-eval.cc, pt-except.cc, pt-exp.cc, pt-fcn-handle.cc, pt-funcall.cc, pt-id.cc, pt-idx.cc, pt-jump.cc, pt-loop.cc, pt-loop.h, pt-mat.cc, pt-misc.cc, pt-pr-code.cc, pt-select.cc, pt-stmt.cc, pt-unop.cc, pt.cc, token.cc, Array-jit.cc, Array-tc.cc, version.cc, Array-C.cc, Array-b.cc, Array-ch.cc, Array-d.cc, Array-f.cc, Array-fC.cc, Array-i.cc, Array-idx-vec.cc, Array-s.cc, Array-str.cc, Array-util.cc, Array-voidp.cc, Array.cc, CColVector.cc, CDiagMatrix.cc, CMatrix.cc, CNDArray.cc, CRowVector.cc, CSparse.cc, MArray-C.cc, MArray-d.cc, MArray-f.cc, MArray-fC.cc, MArray-i.cc, MArray-s.cc, MSparse-C.cc, MSparse-d.cc, MatrixType.cc, PermMatrix.cc, Range.cc, Sparse-C.cc, Sparse-b.cc, Sparse-d.cc, boolMatrix.cc, boolNDArray.cc, boolSparse.cc, chMatrix.cc, chNDArray.cc, dColVector.cc, dDiagMatrix.cc, dMatrix.cc, dNDArray.cc, dRowVector.cc, dSparse.cc, dim-vector.cc, dim-vector.h, fCColVector.cc, fCDiagMatrix.cc, fCMatrix.cc, fCNDArray.cc, fCRowVector.cc, fColVector.cc, fDiagMatrix.cc, fMatrix.cc, fNDArray.cc, fRowVector.cc, idx-vector.cc, int16NDArray.cc, int32NDArray.cc, int64NDArray.cc, int8NDArray.cc, uint16NDArray.cc, uint32NDArray.cc, uint64NDArray.cc, uint8NDArray.cc, Faddeeva.cc, blaswrap.c, cquit.c, f77-extern.cc, f77-fcn.c, f77-fcn.h, lo-error.c, lo-error.h, quit.cc, quit.h, liboctave-build-info.in.cc, CollocWt.cc, DASPK.cc, DASRT.cc, DASSL.cc, EIG.cc, LSODE.cc, ODES.cc, Quad.cc, aepbalance.cc, chol.cc, eigs-base.cc, fEIG.cc, gepbalance.cc, hess.cc, lo-mappers.cc, lo-specfun.cc, lu.cc, oct-convn.cc, oct-fftw.cc, oct-norm.cc, oct-rand.cc, oct-spparms.cc, qr.cc, qrp.cc, randmtzig.cc, randpoisson.cc, schur.cc, sparse-chol.cc, sparse-dmsolve.cc, sparse-lu.cc, sparse-qr.cc, svd.cc, mk-ops.awk, mx-defs.h, dir-ops.cc, file-ops.cc, file-stat.cc, lo-sysdep.cc, mach-info.cc, oct-env.cc, oct-group.cc, oct-passwd.cc, oct-syscalls.cc, oct-time.cc, oct-uname.cc, pathlen.h, syswait.h, cmd-edit.cc, cmd-hist.cc, data-conv.cc, f2c-main.c, glob-match.cc, kpse.cc, lo-array-errwarn.cc, lo-array-gripes.cc, lo-cutils.c, lo-cutils.h, lo-ieee.cc, lo-ieee.h, lo-regexp.cc, lo-utils.cc, oct-base64.cc, oct-glob.cc, oct-inttypes.cc, oct-inttypes.h, oct-locbuf.cc, oct-mutex.cc, oct-rl-edit.c, oct-rl-edit.h, oct-rl-hist.c, oct-rl-hist.h, oct-shlib.cc, oct-sort.cc, pathsearch.cc, singleton-cleanup.cc, sparse-sort.cc, sparse-util.cc, statdefs.h, str-vec.cc, unwind-prot.cc, url-transfer.cc, acinclude.m4, display-available.c, display-available.h, main-cli.cc, main-gui.cc, main.in.cc, mkoctfile.in.cc, octave-build-info.in.cc, octave-config.in.cc, shared-fcns.h: Use "#if ..." consistently instead of "#ifdef" and "#ifndef".
author John W. Eaton <jwe@octave.org>
date Tue, 17 May 2016 12:09:30 -0400
parents 40de9f8f23a6
children 59cadee1c74b
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
457
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
1 /*
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
2
19697
4197fc428c7d maint: Update copyright notices for 2015.
John W. Eaton <jwe@octave.org>
parents: 17769
diff changeset
3 Copyright (C) 1994-2015 John W. Eaton
457
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
4
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
5 This file is part of Octave.
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
6
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
7 Octave is free software; you can redistribute it and/or modify it
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
8 under the terms of the GNU General Public License as published by the
7016
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe]
jwe
parents: 5307
diff changeset
9 Free Software Foundation; either version 3 of the License, or (at your
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe]
jwe
parents: 5307
diff changeset
10 option) any later version.
457
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
11
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
12 Octave is distributed in the hope that it will be useful, but WITHOUT
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
13 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
14 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
15 for more details.
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
16
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
17 You should have received a copy of the GNU General Public License
7016
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe]
jwe
parents: 5307
diff changeset
18 along with Octave; see the file COPYING. If not, see
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe]
jwe
parents: 5307
diff changeset
19 <http://www.gnu.org/licenses/>.
457
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
20
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
21 */
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
22
21724
aba2e6293dd8 use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents: 21301
diff changeset
23 #if defined (HAVE_CONFIG_H)
21301
40de9f8f23a6 Use '#include "config.h"' rather than <config.h>.
Rik <rik@octave.org>
parents: 21273
diff changeset
24 # include "config.h"
457
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
25 #endif
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
26
21273
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
27 #include <cassert>
1631
b0b22b6ce22f [project @ 1995-11-14 20:24:15 by jwe]
jwe
parents: 1545
diff changeset
28
21273
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
29 #include "CMatrix.h"
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
30 #include "dDiagMatrix.h"
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
31 #include "fDiagMatrix.h"
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
32 #include "dMatrix.h"
1847
2ffe49eb95a5 [project @ 1996-02-03 12:47:55 by jwe]
jwe
parents: 1631
diff changeset
33 #include "f77-fcn.h"
21273
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
34 #include "fCMatrix.h"
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
35 #include "fMatrix.h"
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
36 #include "lo-error.h"
10601
3ce0c530a9c9 implement svd_driver
Jaroslav Hajek <highegg@gmail.com>
parents: 10350
diff changeset
37 #include "oct-locbuf.h"
21273
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
38 #include "svd.h"
457
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
39
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
40 extern "C"
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
41 {
4552
6f3382e08a52 [project @ 2003-10-27 20:38:02 by jwe]
jwe
parents: 4192
diff changeset
42 F77_RET_T
6f3382e08a52 [project @ 2003-10-27 20:38:02 by jwe]
jwe
parents: 4192
diff changeset
43 F77_FUNC (dgesvd, DGESVD) (F77_CONST_CHAR_ARG_DECL,
10314
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10258
diff changeset
44 F77_CONST_CHAR_ARG_DECL,
11495
8a5e980da6aa style fixes
John W. Eaton <jwe@octave.org>
parents: 10601
diff changeset
45 const octave_idx_type&, const octave_idx_type&,
8a5e980da6aa style fixes
John W. Eaton <jwe@octave.org>
parents: 10601
diff changeset
46 double*, const octave_idx_type&, double*,
8a5e980da6aa style fixes
John W. Eaton <jwe@octave.org>
parents: 10601
diff changeset
47 double*, const octave_idx_type&, double*,
8a5e980da6aa style fixes
John W. Eaton <jwe@octave.org>
parents: 10601
diff changeset
48 const octave_idx_type&, double*,
8a5e980da6aa style fixes
John W. Eaton <jwe@octave.org>
parents: 10601
diff changeset
49 const octave_idx_type&, octave_idx_type&
10314
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10258
diff changeset
50 F77_CHAR_ARG_LEN_DECL
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10258
diff changeset
51 F77_CHAR_ARG_LEN_DECL);
10601
3ce0c530a9c9 implement svd_driver
Jaroslav Hajek <highegg@gmail.com>
parents: 10350
diff changeset
52
3ce0c530a9c9 implement svd_driver
Jaroslav Hajek <highegg@gmail.com>
parents: 10350
diff changeset
53 F77_RET_T
3ce0c530a9c9 implement svd_driver
Jaroslav Hajek <highegg@gmail.com>
parents: 10350
diff changeset
54 F77_FUNC (dgesdd, DGESDD) (F77_CONST_CHAR_ARG_DECL,
11495
8a5e980da6aa style fixes
John W. Eaton <jwe@octave.org>
parents: 10601
diff changeset
55 const octave_idx_type&, const octave_idx_type&,
8a5e980da6aa style fixes
John W. Eaton <jwe@octave.org>
parents: 10601
diff changeset
56 double*, const octave_idx_type&, double*,
8a5e980da6aa style fixes
John W. Eaton <jwe@octave.org>
parents: 10601
diff changeset
57 double*, const octave_idx_type&, double*,
8a5e980da6aa style fixes
John W. Eaton <jwe@octave.org>
parents: 10601
diff changeset
58 const octave_idx_type&, double*,
8a5e980da6aa style fixes
John W. Eaton <jwe@octave.org>
parents: 10601
diff changeset
59 const octave_idx_type&, octave_idx_type *,
8a5e980da6aa style fixes
John W. Eaton <jwe@octave.org>
parents: 10601
diff changeset
60 octave_idx_type&
10601
3ce0c530a9c9 implement svd_driver
Jaroslav Hajek <highegg@gmail.com>
parents: 10350
diff changeset
61 F77_CHAR_ARG_LEN_DECL);
21273
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
62
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
63 F77_RET_T
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
64 F77_FUNC (sgesvd, SGESVD) (F77_CONST_CHAR_ARG_DECL,
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
65 F77_CONST_CHAR_ARG_DECL,
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
66 const octave_idx_type&, const octave_idx_type&,
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
67 float*, const octave_idx_type&, float*,
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
68 float*, const octave_idx_type&, float*,
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
69 const octave_idx_type&, float*,
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
70 const octave_idx_type&, octave_idx_type&
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
71 F77_CHAR_ARG_LEN_DECL
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
72 F77_CHAR_ARG_LEN_DECL);
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
73
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
74 F77_RET_T
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
75 F77_FUNC (sgesdd, SGESDD) (F77_CONST_CHAR_ARG_DECL,
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
76 const octave_idx_type&, const octave_idx_type&,
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
77 float*, const octave_idx_type&, float*,
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
78 float*, const octave_idx_type&, float*,
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
79 const octave_idx_type&, float*,
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
80 const octave_idx_type&, octave_idx_type *,
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
81 octave_idx_type&
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
82 F77_CHAR_ARG_LEN_DECL);
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
83
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
84 F77_RET_T
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
85 F77_FUNC (zgesvd, ZGESVD) (F77_CONST_CHAR_ARG_DECL,
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
86 F77_CONST_CHAR_ARG_DECL,
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
87 const octave_idx_type&, const octave_idx_type&,
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
88 Complex*, const octave_idx_type&,
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
89 double*, Complex*, const octave_idx_type&,
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
90 Complex*, const octave_idx_type&, Complex*,
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
91 const octave_idx_type&, double*, octave_idx_type&
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
92 F77_CHAR_ARG_LEN_DECL
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
93 F77_CHAR_ARG_LEN_DECL);
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
94
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
95 F77_RET_T
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
96 F77_FUNC (zgesdd, ZGESDD) (F77_CONST_CHAR_ARG_DECL,
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
97 const octave_idx_type&, const octave_idx_type&,
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
98 Complex*, const octave_idx_type&,
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
99 double*, Complex*, const octave_idx_type&,
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
100 Complex*, const octave_idx_type&, Complex*,
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
101 const octave_idx_type&, double*,
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
102 octave_idx_type *, octave_idx_type&
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
103 F77_CHAR_ARG_LEN_DECL);
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
104 F77_RET_T
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
105 F77_FUNC (cgesvd, CGESVD) (F77_CONST_CHAR_ARG_DECL,
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
106 F77_CONST_CHAR_ARG_DECL,
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
107 const octave_idx_type&, const octave_idx_type&,
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
108 FloatComplex*, const octave_idx_type&, float*,
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
109 FloatComplex*, const octave_idx_type&,
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
110 FloatComplex*, const octave_idx_type&,
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
111 FloatComplex*, const octave_idx_type&,
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
112 float*, octave_idx_type&
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
113 F77_CHAR_ARG_LEN_DECL
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
114 F77_CHAR_ARG_LEN_DECL);
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
115
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
116 F77_RET_T
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
117 F77_FUNC (cgesdd, CGESDD) (F77_CONST_CHAR_ARG_DECL,
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
118 const octave_idx_type&, const octave_idx_type&,
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
119 FloatComplex*, const octave_idx_type&, float*,
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
120 FloatComplex*, const octave_idx_type&,
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
121 FloatComplex*, const octave_idx_type&,
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
122 FloatComplex*, const octave_idx_type&,
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
123 float*, octave_idx_type *, octave_idx_type&
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
124 F77_CHAR_ARG_LEN_DECL);
457
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
125 }
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
126
21273
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
127 template <typename T>
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
128 T
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
129 svd<T>::left_singular_matrix (void) const
1543
d6e96e0bc681 [project @ 1995-10-06 05:49:21 by jwe]
jwe
parents: 1368
diff changeset
130 {
21273
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
131 if (type_computed == svd::sigma_only)
21136
7cac4e7458f2 maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents: 19697
diff changeset
132 (*current_liboctave_error_handler)
21273
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
133 ("svd: U not computed because type == svd::sigma_only");
21136
7cac4e7458f2 maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents: 19697
diff changeset
134
7cac4e7458f2 maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents: 19697
diff changeset
135 return left_sm;
1543
d6e96e0bc681 [project @ 1995-10-06 05:49:21 by jwe]
jwe
parents: 1368
diff changeset
136 }
d6e96e0bc681 [project @ 1995-10-06 05:49:21 by jwe]
jwe
parents: 1368
diff changeset
137
21273
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
138 template <typename T>
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
139 T
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
140 svd<T>::right_singular_matrix (void) const
1543
d6e96e0bc681 [project @ 1995-10-06 05:49:21 by jwe]
jwe
parents: 1368
diff changeset
141 {
21273
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
142 if (type_computed == svd::sigma_only)
21136
7cac4e7458f2 maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents: 19697
diff changeset
143 (*current_liboctave_error_handler)
21273
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
144 ("svd: V not computed because type == svd::sigma_only");
21136
7cac4e7458f2 maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents: 19697
diff changeset
145
7cac4e7458f2 maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents: 19697
diff changeset
146 return right_sm;
1543
d6e96e0bc681 [project @ 1995-10-06 05:49:21 by jwe]
jwe
parents: 1368
diff changeset
147 }
d6e96e0bc681 [project @ 1995-10-06 05:49:21 by jwe]
jwe
parents: 1368
diff changeset
148
21273
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
149 template <typename T>
5275
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5105
diff changeset
150 octave_idx_type
21273
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
151 svd<T>::empty_init (octave_idx_type nr, octave_idx_type nc, svd::type svd_type)
457
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
152 {
21273
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
153 assert (nr == 0 || nc == 0);
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
154
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
155 static typename T::element_type zero (0);
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
156 static typename T::element_type one (1);
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
157
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
158 switch (svd_type)
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
159 {
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
160 case svd::std:
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
161 left_sm = T (nr, nr, zero);
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
162 for (octave_idx_type i = 0; i < nr; i++)
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
163 left_sm.xelem (i, i) = one;
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
164 sigma = DM_T (nr, nc);
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
165 right_sm = T (nc, nc, zero);
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
166 for (octave_idx_type i = 0; i < nc; i++)
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
167 right_sm.xelem (i, i) = one;
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
168 break;
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
169
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
170 case svd::economy:
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
171 left_sm = T (nr, 0, zero);
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
172 sigma = DM_T (0, 0);
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
173 right_sm = T (0, nc, zero);
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
174 break;
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
175
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
176 case svd::sigma_only:
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
177 default:
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
178 sigma = DM_T (0, 1);
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
179 break;
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
180 }
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
181
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
182 return 0;
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
183 }
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
184
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
185 // Specializations.
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
186
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
187 template <>
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
188 octave_idx_type
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
189 svd<Matrix>::init (const Matrix& a, svd::type svd_type, svd::driver svd_driver)
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
190 {
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
191 octave_idx_type info = 0;
457
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
192
5275
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5105
diff changeset
193 octave_idx_type m = a.rows ();
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5105
diff changeset
194 octave_idx_type n = a.cols ();
457
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
195
21273
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
196 if (m == 0 || n == 0)
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
197 return empty_init (m, n, svd_type);
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
198
1930
d20ab06301e8 [project @ 1996-02-11 22:30:18 by jwe]
jwe
parents: 1882
diff changeset
199 Matrix atmp = a;
d20ab06301e8 [project @ 1996-02-11 22:30:18 by jwe]
jwe
parents: 1882
diff changeset
200 double *tmp_data = atmp.fortran_vec ();
457
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
201
5275
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5105
diff changeset
202 octave_idx_type min_mn = m < n ? m : n;
457
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
203
1930
d20ab06301e8 [project @ 1996-02-11 22:30:18 by jwe]
jwe
parents: 1882
diff changeset
204 char jobu = 'A';
d20ab06301e8 [project @ 1996-02-11 22:30:18 by jwe]
jwe
parents: 1882
diff changeset
205 char jobv = 'A';
537
4ecbfd3c3710 [project @ 1994-07-21 22:30:34 by jwe]
jwe
parents: 457
diff changeset
206
5275
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5105
diff changeset
207 octave_idx_type ncol_u = m;
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5105
diff changeset
208 octave_idx_type nrow_vt = n;
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5105
diff changeset
209 octave_idx_type nrow_s = m;
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5105
diff changeset
210 octave_idx_type ncol_s = n;
537
4ecbfd3c3710 [project @ 1994-07-21 22:30:34 by jwe]
jwe
parents: 457
diff changeset
211
1543
d6e96e0bc681 [project @ 1995-10-06 05:49:21 by jwe]
jwe
parents: 1368
diff changeset
212 switch (svd_type)
537
4ecbfd3c3710 [project @ 1994-07-21 22:30:34 by jwe]
jwe
parents: 457
diff changeset
213 {
21273
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
214 case svd::economy:
1930
d20ab06301e8 [project @ 1996-02-11 22:30:18 by jwe]
jwe
parents: 1882
diff changeset
215 jobu = jobv = 'S';
537
4ecbfd3c3710 [project @ 1994-07-21 22:30:34 by jwe]
jwe
parents: 457
diff changeset
216 ncol_u = nrow_vt = nrow_s = ncol_s = min_mn;
1543
d6e96e0bc681 [project @ 1995-10-06 05:49:21 by jwe]
jwe
parents: 1368
diff changeset
217 break;
d6e96e0bc681 [project @ 1995-10-06 05:49:21 by jwe]
jwe
parents: 1368
diff changeset
218
21273
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
219 case svd::sigma_only:
2621
337a09dd1c06 [project @ 1997-01-24 21:49:41 by jwe]
jwe
parents: 1993
diff changeset
220
337a09dd1c06 [project @ 1997-01-24 21:49:41 by jwe]
jwe
parents: 1993
diff changeset
221 // Note: for this case, both jobu and jobv should be 'N', but
337a09dd1c06 [project @ 1997-01-24 21:49:41 by jwe]
jwe
parents: 1993
diff changeset
222 // there seems to be a bug in dgesvd from Lapack V2.0. To
337a09dd1c06 [project @ 1997-01-24 21:49:41 by jwe]
jwe
parents: 1993
diff changeset
223 // demonstrate the bug, set both jobu and jobv to 'N' and find
337a09dd1c06 [project @ 1997-01-24 21:49:41 by jwe]
jwe
parents: 1993
diff changeset
224 // the singular values of [eye(3), eye(3)]. The result is
337a09dd1c06 [project @ 1997-01-24 21:49:41 by jwe]
jwe
parents: 1993
diff changeset
225 // [-sqrt(2), -sqrt(2), -sqrt(2)].
3335
f39b97e13cf2 [project @ 1999-11-03 21:05:01 by jwe]
jwe
parents: 2847
diff changeset
226 //
f39b97e13cf2 [project @ 1999-11-03 21:05:01 by jwe]
jwe
parents: 2847
diff changeset
227 // For Lapack 3.0, this problem seems to be fixed.
2621
337a09dd1c06 [project @ 1997-01-24 21:49:41 by jwe]
jwe
parents: 1993
diff changeset
228
15887
8ced82e96b48 Fix segfaults with gesdd driver for svd (bug #37998).
Rik <rik@octave.org>
parents: 14138
diff changeset
229 jobu = jobv = 'N';
1545
8bdfa6fe5d69 [project @ 1995-10-06 06:13:28 by jwe]
jwe
parents: 1544
diff changeset
230 ncol_u = nrow_vt = 1;
1543
d6e96e0bc681 [project @ 1995-10-06 05:49:21 by jwe]
jwe
parents: 1368
diff changeset
231 break;
d6e96e0bc681 [project @ 1995-10-06 05:49:21 by jwe]
jwe
parents: 1368
diff changeset
232
d6e96e0bc681 [project @ 1995-10-06 05:49:21 by jwe]
jwe
parents: 1368
diff changeset
233 default:
d6e96e0bc681 [project @ 1995-10-06 05:49:21 by jwe]
jwe
parents: 1368
diff changeset
234 break;
537
4ecbfd3c3710 [project @ 1994-07-21 22:30:34 by jwe]
jwe
parents: 457
diff changeset
235 }
4ecbfd3c3710 [project @ 1994-07-21 22:30:34 by jwe]
jwe
parents: 457
diff changeset
236
1544
f1931fc63ce9 [project @ 1995-10-06 05:57:56 by jwe]
jwe
parents: 1543
diff changeset
237 type_computed = svd_type;
f1931fc63ce9 [project @ 1995-10-06 05:57:56 by jwe]
jwe
parents: 1543
diff changeset
238
2621
337a09dd1c06 [project @ 1997-01-24 21:49:41 by jwe]
jwe
parents: 1993
diff changeset
239 if (! (jobu == 'N' || jobu == 'O'))
1931
c91f81d5f72c [project @ 1996-02-11 23:02:16 by jwe]
jwe
parents: 1930
diff changeset
240 left_sm.resize (m, ncol_u);
1930
d20ab06301e8 [project @ 1996-02-11 22:30:18 by jwe]
jwe
parents: 1882
diff changeset
241
d20ab06301e8 [project @ 1996-02-11 22:30:18 by jwe]
jwe
parents: 1882
diff changeset
242 double *u = left_sm.fortran_vec ();
d20ab06301e8 [project @ 1996-02-11 22:30:18 by jwe]
jwe
parents: 1882
diff changeset
243
d20ab06301e8 [project @ 1996-02-11 22:30:18 by jwe]
jwe
parents: 1882
diff changeset
244 sigma.resize (nrow_s, ncol_s);
15018
3d8ace26c5b4 maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents: 14138
diff changeset
245 double *s_vec = sigma.fortran_vec ();
1930
d20ab06301e8 [project @ 1996-02-11 22:30:18 by jwe]
jwe
parents: 1882
diff changeset
246
2621
337a09dd1c06 [project @ 1997-01-24 21:49:41 by jwe]
jwe
parents: 1993
diff changeset
247 if (! (jobv == 'N' || jobv == 'O'))
1930
d20ab06301e8 [project @ 1996-02-11 22:30:18 by jwe]
jwe
parents: 1882
diff changeset
248 right_sm.resize (nrow_vt, n);
d20ab06301e8 [project @ 1996-02-11 22:30:18 by jwe]
jwe
parents: 1882
diff changeset
249
d20ab06301e8 [project @ 1996-02-11 22:30:18 by jwe]
jwe
parents: 1882
diff changeset
250 double *vt = right_sm.fortran_vec ();
457
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
251
15887
8ced82e96b48 Fix segfaults with gesdd driver for svd (bug #37998).
Rik <rik@octave.org>
parents: 14138
diff changeset
252 // Query DGESVD for the correct dimension of WORK.
1930
d20ab06301e8 [project @ 1996-02-11 22:30:18 by jwe]
jwe
parents: 1882
diff changeset
253
5275
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5105
diff changeset
254 octave_idx_type lwork = -1;
3336
08ad797989f8 [project @ 1999-11-03 21:41:34 by jwe]
jwe
parents: 3335
diff changeset
255
11570
57632dea2446 attempt better backward compatibility for Array constructors
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
256 Array<double> work (dim_vector (1, 1));
457
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
257
10258
e317791645c4 64-bit fixes
John W. Eaton <jwe@octave.org>
parents: 10185
diff changeset
258 octave_idx_type one = 1;
15887
8ced82e96b48 Fix segfaults with gesdd driver for svd (bug #37998).
Rik <rik@octave.org>
parents: 14138
diff changeset
259 octave_idx_type m1 = std::max (m, one);
8ced82e96b48 Fix segfaults with gesdd driver for svd (bug #37998).
Rik <rik@octave.org>
parents: 14138
diff changeset
260 octave_idx_type nrow_vt1 = std::max (nrow_vt, one);
10185
455759a5fcbe fix norm and svd on empty matrices
Jaroslav Hajek <highegg@gmail.com>
parents: 10158
diff changeset
261
21273
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
262 if (svd_driver == svd::GESVD)
10601
3ce0c530a9c9 implement svd_driver
Jaroslav Hajek <highegg@gmail.com>
parents: 10350
diff changeset
263 {
3ce0c530a9c9 implement svd_driver
Jaroslav Hajek <highegg@gmail.com>
parents: 10350
diff changeset
264 F77_XFCN (dgesvd, DGESVD, (F77_CONST_CHAR_ARG2 (&jobu, 1),
3ce0c530a9c9 implement svd_driver
Jaroslav Hajek <highegg@gmail.com>
parents: 10350
diff changeset
265 F77_CONST_CHAR_ARG2 (&jobv, 1),
3ce0c530a9c9 implement svd_driver
Jaroslav Hajek <highegg@gmail.com>
parents: 10350
diff changeset
266 m, n, tmp_data, m1, s_vec, u, m1, vt,
3ce0c530a9c9 implement svd_driver
Jaroslav Hajek <highegg@gmail.com>
parents: 10350
diff changeset
267 nrow_vt1, work.fortran_vec (), lwork, info
3ce0c530a9c9 implement svd_driver
Jaroslav Hajek <highegg@gmail.com>
parents: 10350
diff changeset
268 F77_CHAR_ARG_LEN (1)
3ce0c530a9c9 implement svd_driver
Jaroslav Hajek <highegg@gmail.com>
parents: 10350
diff changeset
269 F77_CHAR_ARG_LEN (1)));
3ce0c530a9c9 implement svd_driver
Jaroslav Hajek <highegg@gmail.com>
parents: 10350
diff changeset
270
3ce0c530a9c9 implement svd_driver
Jaroslav Hajek <highegg@gmail.com>
parents: 10350
diff changeset
271 lwork = static_cast<octave_idx_type> (work(0));
11574
a83bad07f7e3 attempt better backward compatibility for Array resize functions
John W. Eaton <jwe@octave.org>
parents: 11570
diff changeset
272 work.resize (dim_vector (lwork, 1));
10601
3ce0c530a9c9 implement svd_driver
Jaroslav Hajek <highegg@gmail.com>
parents: 10350
diff changeset
273
3ce0c530a9c9 implement svd_driver
Jaroslav Hajek <highegg@gmail.com>
parents: 10350
diff changeset
274 F77_XFCN (dgesvd, DGESVD, (F77_CONST_CHAR_ARG2 (&jobu, 1),
3ce0c530a9c9 implement svd_driver
Jaroslav Hajek <highegg@gmail.com>
parents: 10350
diff changeset
275 F77_CONST_CHAR_ARG2 (&jobv, 1),
3ce0c530a9c9 implement svd_driver
Jaroslav Hajek <highegg@gmail.com>
parents: 10350
diff changeset
276 m, n, tmp_data, m1, s_vec, u, m1, vt,
3ce0c530a9c9 implement svd_driver
Jaroslav Hajek <highegg@gmail.com>
parents: 10350
diff changeset
277 nrow_vt1, work.fortran_vec (), lwork, info
3ce0c530a9c9 implement svd_driver
Jaroslav Hajek <highegg@gmail.com>
parents: 10350
diff changeset
278 F77_CHAR_ARG_LEN (1)
3ce0c530a9c9 implement svd_driver
Jaroslav Hajek <highegg@gmail.com>
parents: 10350
diff changeset
279 F77_CHAR_ARG_LEN (1)));
1543
d6e96e0bc681 [project @ 1995-10-06 05:49:21 by jwe]
jwe
parents: 1368
diff changeset
280
10601
3ce0c530a9c9 implement svd_driver
Jaroslav Hajek <highegg@gmail.com>
parents: 10350
diff changeset
281 }
21273
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
282 else if (svd_driver == svd::GESDD)
10601
3ce0c530a9c9 implement svd_driver
Jaroslav Hajek <highegg@gmail.com>
parents: 10350
diff changeset
283 {
3ce0c530a9c9 implement svd_driver
Jaroslav Hajek <highegg@gmail.com>
parents: 10350
diff changeset
284 assert (jobu == jobv);
3ce0c530a9c9 implement svd_driver
Jaroslav Hajek <highegg@gmail.com>
parents: 10350
diff changeset
285 char jobz = jobu;
3ce0c530a9c9 implement svd_driver
Jaroslav Hajek <highegg@gmail.com>
parents: 10350
diff changeset
286 OCTAVE_LOCAL_BUFFER (octave_idx_type, iwork, 8*min_mn);
3ce0c530a9c9 implement svd_driver
Jaroslav Hajek <highegg@gmail.com>
parents: 10350
diff changeset
287
3ce0c530a9c9 implement svd_driver
Jaroslav Hajek <highegg@gmail.com>
parents: 10350
diff changeset
288 F77_XFCN (dgesdd, DGESDD, (F77_CONST_CHAR_ARG2 (&jobz, 1),
17769
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
289 m, n, tmp_data, m1, s_vec, u, m1, vt, nrow_vt1,
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
290 work.fortran_vec (), lwork, iwork, info
10601
3ce0c530a9c9 implement svd_driver
Jaroslav Hajek <highegg@gmail.com>
parents: 10350
diff changeset
291 F77_CHAR_ARG_LEN (1)));
3336
08ad797989f8 [project @ 1999-11-03 21:41:34 by jwe]
jwe
parents: 3335
diff changeset
292
10601
3ce0c530a9c9 implement svd_driver
Jaroslav Hajek <highegg@gmail.com>
parents: 10350
diff changeset
293 lwork = static_cast<octave_idx_type> (work(0));
11574
a83bad07f7e3 attempt better backward compatibility for Array resize functions
John W. Eaton <jwe@octave.org>
parents: 11570
diff changeset
294 work.resize (dim_vector (lwork, 1));
10601
3ce0c530a9c9 implement svd_driver
Jaroslav Hajek <highegg@gmail.com>
parents: 10350
diff changeset
295
3ce0c530a9c9 implement svd_driver
Jaroslav Hajek <highegg@gmail.com>
parents: 10350
diff changeset
296 F77_XFCN (dgesdd, DGESDD, (F77_CONST_CHAR_ARG2 (&jobz, 1),
17769
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
297 m, n, tmp_data, m1, s_vec, u, m1, vt, nrow_vt1,
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
298 work.fortran_vec (), lwork, iwork, info
10601
3ce0c530a9c9 implement svd_driver
Jaroslav Hajek <highegg@gmail.com>
parents: 10350
diff changeset
299 F77_CHAR_ARG_LEN (1)));
3ce0c530a9c9 implement svd_driver
Jaroslav Hajek <highegg@gmail.com>
parents: 10350
diff changeset
300
3ce0c530a9c9 implement svd_driver
Jaroslav Hajek <highegg@gmail.com>
parents: 10350
diff changeset
301 }
3ce0c530a9c9 implement svd_driver
Jaroslav Hajek <highegg@gmail.com>
parents: 10350
diff changeset
302 else
21273
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
303 abort ();
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
304
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
305 if (! (jobv == 'N' || jobv == 'O'))
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
306 right_sm = right_sm.transpose ();
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
307
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
308 return info;
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
309 }
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
310
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
311 template <>
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
312 octave_idx_type
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
313 svd<FloatMatrix>::init (const FloatMatrix& a, svd::type svd_type,
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
314 svd::driver svd_driver)
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
315 {
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
316 octave_idx_type info;
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
317
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
318 octave_idx_type m = a.rows ();
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
319 octave_idx_type n = a.cols ();
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
320
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
321 if (m == 0 || n == 0)
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
322 return empty_init (m, n, svd_type);
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
323
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
324 FloatMatrix atmp = a;
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
325 float *tmp_data = atmp.fortran_vec ();
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
326
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
327 octave_idx_type min_mn = m < n ? m : n;
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
328
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
329 char jobu = 'A';
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
330 char jobv = 'A';
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
331
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
332 octave_idx_type ncol_u = m;
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
333 octave_idx_type nrow_vt = n;
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
334 octave_idx_type nrow_s = m;
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
335 octave_idx_type ncol_s = n;
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
336
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
337 switch (svd_type)
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
338 {
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
339 case svd::economy:
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
340 jobu = jobv = 'S';
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
341 ncol_u = nrow_vt = nrow_s = ncol_s = min_mn;
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
342 break;
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
343
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
344 case svd::sigma_only:
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
345
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
346 // Note: for this case, both jobu and jobv should be 'N', but
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
347 // there seems to be a bug in dgesvd from Lapack V2.0. To
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
348 // demonstrate the bug, set both jobu and jobv to 'N' and find
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
349 // the singular values of [eye(3), eye(3)]. The result is
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
350 // [-sqrt(2), -sqrt(2), -sqrt(2)].
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
351 //
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
352 // For Lapack 3.0, this problem seems to be fixed.
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
353
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
354 jobu = jobv = 'N';
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
355 ncol_u = nrow_vt = 1;
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
356 break;
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
357
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
358 default:
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
359 break;
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
360 }
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
361
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
362 type_computed = svd_type;
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
363
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
364 if (! (jobu == 'N' || jobu == 'O'))
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
365 left_sm.resize (m, ncol_u);
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
366
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
367 float *u = left_sm.fortran_vec ();
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
368
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
369 sigma.resize (nrow_s, ncol_s);
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
370 float *s_vec = sigma.fortran_vec ();
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
371
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
372 if (! (jobv == 'N' || jobv == 'O'))
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
373 right_sm.resize (nrow_vt, n);
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
374
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
375 float *vt = right_sm.fortran_vec ();
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
376
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
377 // Query SGESVD for the correct dimension of WORK.
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
378
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
379 octave_idx_type lwork = -1;
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
380
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
381 Array<float> work (dim_vector (1, 1));
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
382
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
383 octave_idx_type one = 1;
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
384 octave_idx_type m1 = std::max (m, one);
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
385 octave_idx_type nrow_vt1 = std::max (nrow_vt, one);
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
386
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
387 if (svd_driver == svd::GESVD)
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
388 {
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
389 F77_XFCN (sgesvd, SGESVD, (F77_CONST_CHAR_ARG2 (&jobu, 1),
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
390 F77_CONST_CHAR_ARG2 (&jobv, 1),
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
391 m, n, tmp_data, m1, s_vec, u, m1, vt,
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
392 nrow_vt1, work.fortran_vec (), lwork, info
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
393 F77_CHAR_ARG_LEN (1)
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
394 F77_CHAR_ARG_LEN (1)));
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
395
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
396 lwork = static_cast<octave_idx_type> (work(0));
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
397 work.resize (dim_vector (lwork, 1));
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
398
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
399 F77_XFCN (sgesvd, SGESVD, (F77_CONST_CHAR_ARG2 (&jobu, 1),
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
400 F77_CONST_CHAR_ARG2 (&jobv, 1),
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
401 m, n, tmp_data, m1, s_vec, u, m1, vt,
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
402 nrow_vt1, work.fortran_vec (), lwork, info
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
403 F77_CHAR_ARG_LEN (1)
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
404 F77_CHAR_ARG_LEN (1)));
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
405
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
406 }
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
407 else if (svd_driver == svd::GESDD)
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
408 {
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
409 assert (jobu == jobv);
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
410 char jobz = jobu;
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
411 OCTAVE_LOCAL_BUFFER (octave_idx_type, iwork, 8*min_mn);
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
412
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
413 F77_XFCN (sgesdd, SGESDD, (F77_CONST_CHAR_ARG2 (&jobz, 1),
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
414 m, n, tmp_data, m1, s_vec, u, m1, vt, nrow_vt1,
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
415 work.fortran_vec (), lwork, iwork, info
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
416 F77_CHAR_ARG_LEN (1)));
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
417
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
418 lwork = static_cast<octave_idx_type> (work(0));
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
419 work.resize (dim_vector (lwork, 1));
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
420
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
421 F77_XFCN (sgesdd, SGESDD, (F77_CONST_CHAR_ARG2 (&jobz, 1),
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
422 m, n, tmp_data, m1, s_vec, u, m1, vt, nrow_vt1,
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
423 work.fortran_vec (), lwork, iwork, info
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
424 F77_CHAR_ARG_LEN (1)));
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
425
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
426 }
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
427 else
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
428 abort ();
3336
08ad797989f8 [project @ 1999-11-03 21:41:34 by jwe]
jwe
parents: 3335
diff changeset
429
7482
29980c6b8604 don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents: 7017
diff changeset
430 if (! (jobv == 'N' || jobv == 'O'))
29980c6b8604 don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents: 7017
diff changeset
431 right_sm = right_sm.transpose ();
457
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
432
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
433 return info;
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
434 }
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
435
21273
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
436 template <>
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
437 octave_idx_type
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
438 svd<ComplexMatrix>::init (const ComplexMatrix& a, svd::type svd_type,
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
439 svd::driver svd_driver)
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
440 {
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
441 octave_idx_type info;
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
442
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
443 octave_idx_type m = a.rows ();
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
444 octave_idx_type n = a.cols ();
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
445
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
446 if (m == 0 || n == 0)
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
447 return empty_init (m, n, svd_type);
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
448
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
449 ComplexMatrix atmp = a;
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
450 Complex *tmp_data = atmp.fortran_vec ();
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
451
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
452 octave_idx_type min_mn = m < n ? m : n;
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
453 octave_idx_type max_mn = m > n ? m : n;
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
454
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
455 char jobu = 'A';
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
456 char jobv = 'A';
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
457
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
458 octave_idx_type ncol_u = m;
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
459 octave_idx_type nrow_vt = n;
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
460 octave_idx_type nrow_s = m;
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
461 octave_idx_type ncol_s = n;
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
462
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
463 switch (svd_type)
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
464 {
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
465 case svd::economy:
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
466 jobu = jobv = 'S';
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
467 ncol_u = nrow_vt = nrow_s = ncol_s = min_mn;
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
468 break;
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
469
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
470 case svd::sigma_only:
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
471
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
472 // Note: for this case, both jobu and jobv should be 'N', but
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
473 // there seems to be a bug in dgesvd from Lapack V2.0. To
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
474 // demonstrate the bug, set both jobu and jobv to 'N' and find
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
475 // the singular values of [eye(3), eye(3)]. The result is
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
476 // [-sqrt(2), -sqrt(2), -sqrt(2)].
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
477 //
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
478 // For Lapack 3.0, this problem seems to be fixed.
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
479
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
480 jobu = jobv = 'N';
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
481 ncol_u = nrow_vt = 1;
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
482 break;
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
483
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
484 default:
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
485 break;
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
486 }
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
487
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
488 type_computed = svd_type;
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
489
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
490 if (! (jobu == 'N' || jobu == 'O'))
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
491 left_sm.resize (m, ncol_u);
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
492
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
493 Complex *u = left_sm.fortran_vec ();
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
494
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
495 sigma.resize (nrow_s, ncol_s);
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
496 double *s_vec = sigma.fortran_vec ();
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
497
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
498 if (! (jobv == 'N' || jobv == 'O'))
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
499 right_sm.resize (nrow_vt, n);
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
500
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
501 Complex *vt = right_sm.fortran_vec ();
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
502
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
503 // Query ZGESVD for the correct dimension of WORK.
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
504
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
505 octave_idx_type lwork = -1;
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
506
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
507 Array<Complex> work (dim_vector (1, 1));
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
508
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
509 octave_idx_type one = 1;
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
510 octave_idx_type m1 = std::max (m, one);
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
511 octave_idx_type nrow_vt1 = std::max (nrow_vt, one);
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
512
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
513 if (svd_driver == svd::GESVD)
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
514 {
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
515 octave_idx_type lrwork = 5*max_mn;
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
516 Array<double> rwork (dim_vector (lrwork, 1));
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
517
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
518 F77_XFCN (zgesvd, ZGESVD, (F77_CONST_CHAR_ARG2 (&jobu, 1),
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
519 F77_CONST_CHAR_ARG2 (&jobv, 1),
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
520 m, n, tmp_data, m1, s_vec, u, m1, vt,
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
521 nrow_vt1, work.fortran_vec (), lwork,
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
522 rwork.fortran_vec (), info
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
523 F77_CHAR_ARG_LEN (1)
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
524 F77_CHAR_ARG_LEN (1)));
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
525
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
526 lwork = static_cast<octave_idx_type> (work(0).real ());
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
527 work.resize (dim_vector (lwork, 1));
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
528
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
529 F77_XFCN (zgesvd, ZGESVD, (F77_CONST_CHAR_ARG2 (&jobu, 1),
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
530 F77_CONST_CHAR_ARG2 (&jobv, 1),
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
531 m, n, tmp_data, m1, s_vec, u, m1, vt,
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
532 nrow_vt1, work.fortran_vec (), lwork,
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
533 rwork.fortran_vec (), info
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
534 F77_CHAR_ARG_LEN (1)
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
535 F77_CHAR_ARG_LEN (1)));
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
536 }
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
537 else if (svd_driver == svd::GESDD)
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
538 {
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
539 assert (jobu == jobv);
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
540 char jobz = jobu;
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
541
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
542 octave_idx_type lrwork;
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
543 if (jobz == 'N')
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
544 lrwork = 7*min_mn;
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
545 else
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
546 lrwork = 5*min_mn*min_mn + 5*min_mn;
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
547 Array<double> rwork (dim_vector (lrwork, 1));
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
548
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
549 OCTAVE_LOCAL_BUFFER (octave_idx_type, iwork, 8*min_mn);
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
550
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
551 F77_XFCN (zgesdd, ZGESDD, (F77_CONST_CHAR_ARG2 (&jobz, 1),
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
552 m, n, tmp_data, m1, s_vec, u, m1, vt,
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
553 nrow_vt1, work.fortran_vec (), lwork,
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
554 rwork.fortran_vec (), iwork, info
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
555 F77_CHAR_ARG_LEN (1)));
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
556
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
557 lwork = static_cast<octave_idx_type> (work(0).real ());
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
558 work.resize (dim_vector (lwork, 1));
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
559
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
560 F77_XFCN (zgesdd, ZGESDD, (F77_CONST_CHAR_ARG2 (&jobz, 1),
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
561 m, n, tmp_data, m1, s_vec, u, m1, vt,
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
562 nrow_vt1, work.fortran_vec (), lwork,
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
563 rwork.fortran_vec (), iwork, info
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
564 F77_CHAR_ARG_LEN (1)));
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
565 }
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
566 else
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
567 abort ();
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
568
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
569 if (! (jobv == 'N' || jobv == 'O'))
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
570 right_sm = right_sm.hermitian ();
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
571
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
572 return info;
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
573 }
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
574
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
575 template <>
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
576 octave_idx_type
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
577 svd<FloatComplexMatrix>::init (const FloatComplexMatrix& a, svd::type svd_type,
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
578 svd::driver svd_driver)
457
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
579 {
21273
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
580 octave_idx_type info;
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
581
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
582 octave_idx_type m = a.rows ();
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
583 octave_idx_type n = a.cols ();
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
584
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
585 if (m == 0 || n == 0)
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
586 return empty_init (m, n, svd_type);
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
587
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
588 FloatComplexMatrix atmp = a;
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
589 FloatComplex *tmp_data = atmp.fortran_vec ();
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
590
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
591 octave_idx_type min_mn = m < n ? m : n;
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
592 octave_idx_type max_mn = m > n ? m : n;
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
593
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
594 char jobu = 'A';
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
595 char jobv = 'A';
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
596
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
597 octave_idx_type ncol_u = m;
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
598 octave_idx_type nrow_vt = n;
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
599 octave_idx_type nrow_s = m;
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
600 octave_idx_type ncol_s = n;
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
601
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
602 switch (svd_type)
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
603 {
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
604 case svd::economy:
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
605 jobu = jobv = 'S';
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
606 ncol_u = nrow_vt = nrow_s = ncol_s = min_mn;
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
607 break;
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
608
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
609 case svd::sigma_only:
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
610
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
611 // Note: for this case, both jobu and jobv should be 'N', but
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
612 // there seems to be a bug in dgesvd from Lapack V2.0. To
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
613 // demonstrate the bug, set both jobu and jobv to 'N' and find
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
614 // the singular values of [eye(3), eye(3)]. The result is
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
615 // [-sqrt(2), -sqrt(2), -sqrt(2)].
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
616 //
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
617 // For Lapack 3.0, this problem seems to be fixed.
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
618
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
619 jobu = jobv = 'N';
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
620 ncol_u = nrow_vt = 1;
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
621 break;
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
622
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
623 default:
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
624 break;
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
625 }
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
626
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
627 type_computed = svd_type;
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
628
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
629 if (! (jobu == 'N' || jobu == 'O'))
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
630 left_sm.resize (m, ncol_u);
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
631
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
632 FloatComplex *u = left_sm.fortran_vec ();
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
633
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
634 sigma.resize (nrow_s, ncol_s);
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
635 float *s_vec = sigma.fortran_vec ();
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
636
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
637 if (! (jobv == 'N' || jobv == 'O'))
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
638 right_sm.resize (nrow_vt, n);
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
639
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
640 FloatComplex *vt = right_sm.fortran_vec ();
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
641
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
642 // Query CGESVD for the correct dimension of WORK.
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
643
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
644 octave_idx_type lwork = -1;
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
645
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
646 Array<FloatComplex> work (dim_vector (1, 1));
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
647
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
648 octave_idx_type one = 1;
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
649 octave_idx_type m1 = std::max (m, one);
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
650 octave_idx_type nrow_vt1 = std::max (nrow_vt, one);
457
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
651
21273
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
652 if (svd_driver == svd::GESVD)
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
653 {
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
654 octave_idx_type lrwork = 5*max_mn;
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
655 Array<float> rwork (dim_vector (lrwork, 1));
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
656
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
657 F77_XFCN (cgesvd, CGESVD, (F77_CONST_CHAR_ARG2 (&jobu, 1),
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
658 F77_CONST_CHAR_ARG2 (&jobv, 1),
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
659 m, n, tmp_data, m1, s_vec, u, m1, vt,
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
660 nrow_vt1, work.fortran_vec (), lwork,
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
661 rwork.fortran_vec (), info
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
662 F77_CHAR_ARG_LEN (1)
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
663 F77_CHAR_ARG_LEN (1)));
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
664
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
665 lwork = static_cast<octave_idx_type> (work(0).real ());
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
666 work.resize (dim_vector (lwork, 1));
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
667
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
668 F77_XFCN (cgesvd, CGESVD, (F77_CONST_CHAR_ARG2 (&jobu, 1),
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
669 F77_CONST_CHAR_ARG2 (&jobv, 1),
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
670 m, n, tmp_data, m1, s_vec, u, m1, vt,
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
671 nrow_vt1, work.fortran_vec (), lwork,
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
672 rwork.fortran_vec (), info
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
673 F77_CHAR_ARG_LEN (1)
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
674 F77_CHAR_ARG_LEN (1)));
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
675 }
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
676 else if (svd_driver == svd::GESDD)
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
677 {
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
678 assert (jobu == jobv);
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
679 char jobz = jobu;
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
680
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
681 octave_idx_type lrwork;
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
682 if (jobz == 'N')
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
683 lrwork = 5*min_mn;
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
684 else
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
685 lrwork = min_mn * std::max (5*min_mn+7, 2*max_mn+2*min_mn+1);
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
686 Array<float> rwork (dim_vector (lrwork, 1));
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
687
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
688 OCTAVE_LOCAL_BUFFER (octave_idx_type, iwork, 8*min_mn);
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
689
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
690 F77_XFCN (cgesdd, CGESDD, (F77_CONST_CHAR_ARG2 (&jobz, 1),
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
691 m, n, tmp_data, m1, s_vec, u, m1, vt,
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
692 nrow_vt1, work.fortran_vec (), lwork,
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
693 rwork.fortran_vec (), iwork, info
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
694 F77_CHAR_ARG_LEN (1)));
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
695
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
696 lwork = static_cast<octave_idx_type> (work(0).real ());
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
697 work.resize (dim_vector (lwork, 1));
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
698
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
699 F77_XFCN (cgesdd, CGESDD, (F77_CONST_CHAR_ARG2 (&jobz, 1),
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
700 m, n, tmp_data, m1, s_vec, u, m1, vt,
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
701 nrow_vt1, work.fortran_vec (), lwork,
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
702 rwork.fortran_vec (), iwork, info
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
703 F77_CHAR_ARG_LEN (1)));
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
704 }
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
705 else
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
706 abort ();
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
707
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
708 if (! (jobv == 'N' || jobv == 'O'))
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
709 right_sm = right_sm.hermitian ();
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
710
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
711 return info;
457
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
712 }
21273
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
713
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
714 // Instantiations we need.
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
715
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
716 template class svd<Matrix>;
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
717
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
718 template class svd<FloatMatrix>;
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
719
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
720 template class svd<ComplexMatrix>;
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
721
cbced1c09916 better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
722 template class svd<FloatComplexMatrix>;