octave: liboctave/numeric/svd.cc annotate

annotate liboctave/numeric/svd.cc @ 22135:407c66ae1e20

reduce warnings from GCC's link-time optimization feature (bug #48531) * f77-fcn.h (F77_DBLE_CMPLX, F77_CMPLX): Use C types instead of typedefs for std::complex<T> types. (F77_CMPLX_ARG, F77_CONST_CMPLX_ARG, F77_DBLE_CMPLX_ARG, F77_CONST_DBLE_CMPLX_ARG): New macros. * dot.cc, ordschur.cc, qz.cc, CColVector.cc, CMatrix.cc, CRowVector.cc, CSparse.cc, dSparse.cc, fCColVector.cc, fCMatrix.cc, fCRowVector.cc, f77-fcn.h, EIG.cc, aepbalance.cc, chol.cc, eigs-base.cc, fEIG.cc, gepbalance.cc, hess.cc, lo-specfun.cc, lu.cc, oct-convn.cc, qr.cc, qrp.cc, schur.cc, svd.cc: Use new macros for passing complex arguments to Fortran function. Always pass pointers to complex arguments.

author	John W. Eaton <jwe@octave.org>
date	Mon, 18 Jul 2016 09:38:57 -0400
parents	59cadee1c74b
children	469c817eb256

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,
22133 59cadee1c74b new macros for F77 data types John W. Eaton <jwe@octave.org> parents: 21724 diff changeset	45 const F77_INT&, const F77_INT&,
59cadee1c74b new macros for F77 data types John W. Eaton <jwe@octave.org> parents: 21724 diff changeset	46 F77_DBLE, const F77_INT&, F77_DBLE,
59cadee1c74b new macros for F77 data types John W. Eaton <jwe@octave.org> parents: 21724 diff changeset	47 F77_DBLE, const F77_INT&, F77_DBLE,
59cadee1c74b new macros for F77 data types John W. Eaton <jwe@octave.org> parents: 21724 diff changeset	48 const F77_INT&, F77_DBLE*,
59cadee1c74b new macros for F77 data types John W. Eaton <jwe@octave.org> parents: 21724 diff changeset	49 const F77_INT&, F77_INT&
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,
22133 59cadee1c74b new macros for F77 data types John W. Eaton <jwe@octave.org> parents: 21724 diff changeset	55 const F77_INT&, const F77_INT&,
59cadee1c74b new macros for F77 data types John W. Eaton <jwe@octave.org> parents: 21724 diff changeset	56 F77_DBLE, const F77_INT&, F77_DBLE,
59cadee1c74b new macros for F77 data types John W. Eaton <jwe@octave.org> parents: 21724 diff changeset	57 F77_DBLE, const F77_INT&, F77_DBLE,
59cadee1c74b new macros for F77 data types John W. Eaton <jwe@octave.org> parents: 21724 diff changeset	58 const F77_INT&, F77_DBLE*,
59cadee1c74b new macros for F77 data types John W. Eaton <jwe@octave.org> parents: 21724 diff changeset	59 const F77_INT&, F77_INT *,
59cadee1c74b new macros for F77 data types John W. Eaton <jwe@octave.org> parents: 21724 diff changeset	60 F77_INT&
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,
22133 59cadee1c74b new macros for F77 data types John W. Eaton <jwe@octave.org> parents: 21724 diff changeset	66 const F77_INT&, const F77_INT&,
59cadee1c74b new macros for F77 data types John W. Eaton <jwe@octave.org> parents: 21724 diff changeset	67 F77_REAL, const F77_INT&, F77_REAL,
59cadee1c74b new macros for F77 data types John W. Eaton <jwe@octave.org> parents: 21724 diff changeset	68 F77_REAL, const F77_INT&, F77_REAL,
59cadee1c74b new macros for F77 data types John W. Eaton <jwe@octave.org> parents: 21724 diff changeset	69 const F77_INT&, F77_REAL*,
59cadee1c74b new macros for F77 data types John W. Eaton <jwe@octave.org> parents: 21724 diff changeset	70 const F77_INT&, F77_INT&
21273 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,
22133 59cadee1c74b new macros for F77 data types John W. Eaton <jwe@octave.org> parents: 21724 diff changeset	76 const F77_INT&, const F77_INT&,
59cadee1c74b new macros for F77 data types John W. Eaton <jwe@octave.org> parents: 21724 diff changeset	77 F77_REAL, const F77_INT&, F77_REAL,
59cadee1c74b new macros for F77 data types John W. Eaton <jwe@octave.org> parents: 21724 diff changeset	78 F77_REAL, const F77_INT&, F77_REAL,
59cadee1c74b new macros for F77 data types John W. Eaton <jwe@octave.org> parents: 21724 diff changeset	79 const F77_INT&, F77_REAL*,
59cadee1c74b new macros for F77 data types John W. Eaton <jwe@octave.org> parents: 21724 diff changeset	80 const F77_INT&, F77_INT *,
59cadee1c74b new macros for F77 data types John W. Eaton <jwe@octave.org> parents: 21724 diff changeset	81 F77_INT&
21273 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,
22133 59cadee1c74b new macros for F77 data types John W. Eaton <jwe@octave.org> parents: 21724 diff changeset	87 const F77_INT&, const F77_INT&,
59cadee1c74b new macros for F77 data types John W. Eaton <jwe@octave.org> parents: 21724 diff changeset	88 F77_DBLE_CMPLX*, const F77_INT&,
59cadee1c74b new macros for F77 data types John W. Eaton <jwe@octave.org> parents: 21724 diff changeset	89 F77_DBLE, F77_DBLE_CMPLX, const F77_INT&,
59cadee1c74b new macros for F77 data types John W. Eaton <jwe@octave.org> parents: 21724 diff changeset	90 F77_DBLE_CMPLX, const F77_INT&, F77_DBLE_CMPLX,
59cadee1c74b new macros for F77 data types John W. Eaton <jwe@octave.org> parents: 21724 diff changeset	91 const F77_INT&, F77_DBLE*, F77_INT&
21273 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,
22133 59cadee1c74b new macros for F77 data types John W. Eaton <jwe@octave.org> parents: 21724 diff changeset	97 const F77_INT&, const F77_INT&,
59cadee1c74b new macros for F77 data types John W. Eaton <jwe@octave.org> parents: 21724 diff changeset	98 F77_DBLE_CMPLX*, const F77_INT&,
59cadee1c74b new macros for F77 data types John W. Eaton <jwe@octave.org> parents: 21724 diff changeset	99 F77_DBLE, F77_DBLE_CMPLX, const F77_INT&,
59cadee1c74b new macros for F77 data types John W. Eaton <jwe@octave.org> parents: 21724 diff changeset	100 F77_DBLE_CMPLX, const F77_INT&, F77_DBLE_CMPLX,
59cadee1c74b new macros for F77 data types John W. Eaton <jwe@octave.org> parents: 21724 diff changeset	101 const F77_INT&, F77_DBLE*,
59cadee1c74b new macros for F77 data types John W. Eaton <jwe@octave.org> parents: 21724 diff changeset	102 F77_INT *, F77_INT&
21273 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,
22133 59cadee1c74b new macros for F77 data types John W. Eaton <jwe@octave.org> parents: 21724 diff changeset	107 const F77_INT&, const F77_INT&,
59cadee1c74b new macros for F77 data types John W. Eaton <jwe@octave.org> parents: 21724 diff changeset	108 F77_CMPLX, const F77_INT&, F77_REAL,
59cadee1c74b new macros for F77 data types John W. Eaton <jwe@octave.org> parents: 21724 diff changeset	109 F77_CMPLX*, const F77_INT&,
59cadee1c74b new macros for F77 data types John W. Eaton <jwe@octave.org> parents: 21724 diff changeset	110 F77_CMPLX*, const F77_INT&,
59cadee1c74b new macros for F77 data types John W. Eaton <jwe@octave.org> parents: 21724 diff changeset	111 F77_CMPLX*, const F77_INT&,
59cadee1c74b new macros for F77 data types John W. Eaton <jwe@octave.org> parents: 21724 diff changeset	112 F77_REAL*, F77_INT&
21273 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,
22133 59cadee1c74b new macros for F77 data types John W. Eaton <jwe@octave.org> parents: 21724 diff changeset	118 const F77_INT&, const F77_INT&,
59cadee1c74b new macros for F77 data types John W. Eaton <jwe@octave.org> parents: 21724 diff changeset	119 F77_CMPLX, const F77_INT&, F77_REAL,
59cadee1c74b new macros for F77 data types John W. Eaton <jwe@octave.org> parents: 21724 diff changeset	120 F77_CMPLX*, const F77_INT&,
59cadee1c74b new macros for F77 data types John W. Eaton <jwe@octave.org> parents: 21724 diff changeset	121 F77_CMPLX*, const F77_INT&,
59cadee1c74b new macros for F77 data types John W. Eaton <jwe@octave.org> parents: 21724 diff changeset	122 F77_CMPLX*, const F77_INT&,
59cadee1c74b new macros for F77 data types John W. Eaton <jwe@octave.org> parents: 21724 diff changeset	123 F77_REAL, F77_INT , F77_INT&
21273 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),
22135 407c66ae1e20 reduce warnings from GCC's link-time optimization feature (bug #48531) John W. Eaton <jwe@octave.org> parents: 22133 diff changeset	520 m, n, F77_DBLE_CMPLX_ARG (tmp_data), m1, s_vec, F77_DBLE_CMPLX_ARG (u), m1, F77_DBLE_CMPLX_ARG (vt),
407c66ae1e20 reduce warnings from GCC's link-time optimization feature (bug #48531) John W. Eaton <jwe@octave.org> parents: 22133 diff changeset	521 nrow_vt1, F77_DBLE_CMPLX_ARG (work.fortran_vec ()), lwork,
21273 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),
22135 407c66ae1e20 reduce warnings from GCC's link-time optimization feature (bug #48531) John W. Eaton <jwe@octave.org> parents: 22133 diff changeset	531 m, n, F77_DBLE_CMPLX_ARG (tmp_data), m1, s_vec, F77_DBLE_CMPLX_ARG (u), m1, F77_DBLE_CMPLX_ARG (vt),
407c66ae1e20 reduce warnings from GCC's link-time optimization feature (bug #48531) John W. Eaton <jwe@octave.org> parents: 22133 diff changeset	532 nrow_vt1, F77_DBLE_CMPLX_ARG (work.fortran_vec ()), lwork,
21273 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 = 5min_mnmin_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),
22135 407c66ae1e20 reduce warnings from GCC's link-time optimization feature (bug #48531) John W. Eaton <jwe@octave.org> parents: 22133 diff changeset	552 m, n, F77_DBLE_CMPLX_ARG (tmp_data), m1, s_vec, F77_DBLE_CMPLX_ARG (u), m1, F77_DBLE_CMPLX_ARG (vt),
407c66ae1e20 reduce warnings from GCC's link-time optimization feature (bug #48531) John W. Eaton <jwe@octave.org> parents: 22133 diff changeset	553 nrow_vt1, F77_DBLE_CMPLX_ARG (work.fortran_vec ()), lwork,
21273 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),
22135 407c66ae1e20 reduce warnings from GCC's link-time optimization feature (bug #48531) John W. Eaton <jwe@octave.org> parents: 22133 diff changeset	561 m, n, F77_DBLE_CMPLX_ARG (tmp_data), m1, s_vec, F77_DBLE_CMPLX_ARG (u), m1, F77_DBLE_CMPLX_ARG (vt),
407c66ae1e20 reduce warnings from GCC's link-time optimization feature (bug #48531) John W. Eaton <jwe@octave.org> parents: 22133 diff changeset	562 nrow_vt1, F77_DBLE_CMPLX_ARG (work.fortran_vec ()), lwork,
21273 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),
22135 407c66ae1e20 reduce warnings from GCC's link-time optimization feature (bug #48531) John W. Eaton <jwe@octave.org> parents: 22133 diff changeset	659 m, n, F77_CMPLX_ARG (tmp_data), m1, s_vec, F77_CMPLX_ARG (u), m1, F77_CMPLX_ARG (vt),
407c66ae1e20 reduce warnings from GCC's link-time optimization feature (bug #48531) John W. Eaton <jwe@octave.org> parents: 22133 diff changeset	660 nrow_vt1, F77_CMPLX_ARG (work.fortran_vec ()), lwork,
21273 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),
22135 407c66ae1e20 reduce warnings from GCC's link-time optimization feature (bug #48531) John W. Eaton <jwe@octave.org> parents: 22133 diff changeset	670 m, n, F77_CMPLX_ARG (tmp_data), m1, s_vec, F77_CMPLX_ARG (u), m1, F77_CMPLX_ARG (vt),
407c66ae1e20 reduce warnings from GCC's link-time optimization feature (bug #48531) John W. Eaton <jwe@octave.org> parents: 22133 diff changeset	671 nrow_vt1, F77_CMPLX_ARG (work.fortran_vec ()), lwork,
21273 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 (5min_mn+7, 2max_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),
22135 407c66ae1e20 reduce warnings from GCC's link-time optimization feature (bug #48531) John W. Eaton <jwe@octave.org> parents: 22133 diff changeset	691 m, n, F77_CMPLX_ARG (tmp_data), m1, s_vec, F77_CMPLX_ARG (u), m1, F77_CMPLX_ARG (vt),
407c66ae1e20 reduce warnings from GCC's link-time optimization feature (bug #48531) John W. Eaton <jwe@octave.org> parents: 22133 diff changeset	692 nrow_vt1, F77_CMPLX_ARG (work.fortran_vec ()), lwork,
21273 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),
22135 407c66ae1e20 reduce warnings from GCC's link-time optimization feature (bug #48531) John W. Eaton <jwe@octave.org> parents: 22133 diff changeset	700 m, n, F77_CMPLX_ARG (tmp_data), m1, s_vec, F77_CMPLX_ARG (u), m1, F77_CMPLX_ARG (vt),
407c66ae1e20 reduce warnings from GCC's link-time optimization feature (bug #48531) John W. Eaton <jwe@octave.org> parents: 22133 diff changeset	701 nrow_vt1, F77_CMPLX_ARG (work.fortran_vec ()), lwork,
21273 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>;

Mercurial > octave

annotate liboctave/numeric/svd.cc @ 22135:407c66ae1e20