Mercurial > octave
changeset 22236:065a44375723
gsvd: reduce code duplication with templates.
* CmplxGSVD.cc, CmplxGSVD.h, dbleGSVD.cc, dbleGSVD.h: Remove files for
no longer existing classes. Replaced by gsvd template class. This
classes never existed in an Octave release, this was freshly imported
from Octave Forge so backwards compatibility is not an issue.
* liboctave/numeric/gsvd.h, liboctave/numeric/gsvd.cc: New files for gsvd
class template generated from CmplxGSVD.cc, CmplxGSVD.h, dbleGSVD.cc, and
dbleGSVD.h and converted to template. Removed unused << operator, unused
constructor with &info, and commented code. Only instantiated for Matrix
and ComplexMatrix, providing interface to DGGSVD and ZGGSVD.
* liboctave/numeric/module.mk: Update.
* mx-defs.h, mx-ext.h: Use new classes.
| author | Barbara Locsi <locsi.barbara@gmail.com> |
|---|---|
| date | Tue, 09 Aug 2016 18:02:11 +0200 |
| parents | 63b41167ef1e |
| children | 867b177af1fe |
| files | libinterp/corefcn/gsvd.cc liboctave/numeric/CmplxGSVD.cc liboctave/numeric/CmplxGSVD.h liboctave/numeric/dbleGSVD.cc liboctave/numeric/dbleGSVD.h liboctave/numeric/gsvd.cc liboctave/numeric/gsvd.h liboctave/numeric/module.mk liboctave/operators/mx-defs.h liboctave/operators/mx-ext.h |
| diffstat | 10 files changed, 486 insertions(+), 796 deletions(-) [+] |
line wrap: on
line diff
--- a/libinterp/corefcn/gsvd.cc Thu Aug 04 07:50:31 2016 +0200 +++ b/libinterp/corefcn/gsvd.cc Tue Aug 09 18:02:11 2016 +0200 @@ -1,5 +1,6 @@ // Copyright (C) 1996, 1997 John W. Eaton // Copyright (C) 2006, 2010 Pascal Dupuis <Pascal.Dupuis@uclouvain.be> +// Copyright (C) 2016 Barbara Lócsi // // This program is free software; you can redistribute it and/or modify it under // the terms of the GNU General Public License as published by the Free Software @@ -25,8 +26,17 @@ #include "utils.h" #include "ovl.h" -#include "CmplxGSVD.h" -#include "dbleGSVD.h" +#include "gsvd.h" + +template <typename T> +static typename gsvd<T>::Type +gsvd_type (int nargout) +{ + return ((nargout == 0 || nargout == 1) + ? gsvd<T>::Type::sigma_only + : (nargout > 5) ? gsvd<T>::Type::std : gsvd<T>::Type::economy); +} + DEFUN (gsvd, args, nargout, doc: /* -*- texinfo -*- @@ -86,7 +96,7 @@ @end ifinfo The common upper triangular right term. Other authors, like S. Van Huffel, define this transformation as the simulatenous diagonalisation of the -input matrices, this can be achieved by multiplying +input matrices, this can be achieved by multiplying @iftex @tex X @@ -185,7 +195,7 @@ // octave_idx_type nn = argB.rows (); octave_idx_type np = argB.columns (); - + if (nr == 0 || nc == 0) { if (nargout == 5) @@ -208,10 +218,6 @@ return retval; } - GSVD::type type = ((nargout == 0 || nargout == 1) - ? GSVD::sigma_only - : (nargout > 5) ? GSVD::std : GSVD::economy ); - if (argA.is_real_type () && argB.is_real_type ()) { Matrix tmpA = argA.matrix_value (); @@ -221,17 +227,17 @@ { if (tmpA.any_element_is_inf_or_nan ()) { - error ("gsvd: cannot take GSVD of matrix containing Inf or NaN values"); - return retval; - } - - if (tmpB.any_element_is_inf_or_nan ()) - { - error ("gsvd: cannot take GSVD of matrix containing Inf or NaN values"); + error ("gsvd: cannot take gsvd of matrix containing Inf or NaN values"); return retval; } - GSVD result (tmpA, tmpB, type); + if (tmpB.any_element_is_inf_or_nan ()) + { + error ("gsvd: cannot take gsvd of matrix containing Inf or NaN values"); + return retval; + } + + gsvd<Matrix> result (tmpA, tmpB, gsvd_type<Matrix> (nargout)); // DiagMatrix sigma = result.singular_values (); @@ -244,7 +250,7 @@ retval = ovl (sigA.diag()); } else - { + { if (nargout > 5) retval = ovl (result.left_singular_matrix_A (), result.left_singular_matrix_B (), @@ -270,16 +276,16 @@ { if (ctmpA.any_element_is_inf_or_nan ()) { - error ("gsvd: cannot take GSVD of matrix containing Inf or NaN values"); + error ("gsvd: cannot take gsvd of matrix containing Inf or NaN values"); return retval; } if (ctmpB.any_element_is_inf_or_nan ()) { - error ("gsvd: cannot take GSVD of matrix containing Inf or NaN values"); + error ("gsvd: cannot take gsvd of matrix containing Inf or NaN values"); return retval; } - ComplexGSVD result (ctmpA, ctmpB, type); + gsvd<ComplexMatrix> result (ctmpA, ctmpB, gsvd_type<ComplexMatrix> (nargout)); // DiagMatrix sigma = result.singular_values (); @@ -370,7 +376,7 @@ %! [U, V, C, S, X, R] = gsvd(A, B); %! D1 = [C zeros(3,2)]; -%! D2 = [S zeros(3,2); zeros(2, 3) eye(2)]; +%! D2 = [S zeros(3,2); zeros(2, 3) eye(2)]; %!assert(norm(diag(C).^2+diag(S).^2 - ones(3, 1)) <= 1e-6) %!assert(norm((U'*A*X)-D1*R) <= 1e-6) %!assert(norm((V'*B*X)-D2*R) <= 1e-6) @@ -453,13 +459,13 @@ %!assert(norm((V'*B*X)-D2*[zeros(2, 1) R]) <= 1e-6) %! # now, A is 3x5 -%! A = A0.'; B0=diag([1 2 4 8 16])+j*diag([-5 4 -3 2 -1]); +%! A = A0.'; B0=diag([1 2 4 8 16])+j*diag([-5 4 -3 2 -1]); %! B = B0; %! # disp('Complex: Full rank, 3x5 by 5x5 matrices'); %! # disp([rank(A) rank(B) rank([A' B'])]); %! [U, V, C, S, X, R] = gsvd(A, B); %! D1 = [C zeros(3,2)]; -%! D2 = [S zeros(3,2); zeros(2, 3) eye(2)]; +%! D2 = [S zeros(3,2); zeros(2, 3) eye(2)]; %!assert(norm(diag(C).^2+diag(S).^2 - ones(3, 1)) <= 1e-6) %!assert(norm((U'*A*X)-D1*R) <= 1e-6) %!assert(norm((V'*B*X)-D2*R) <= 1e-6) @@ -497,4 +503,4 @@ %!assert(norm((U'*A*X)-D1*[zeros(4, 1) R]) <= 1e-6) %!assert(norm((V'*B*X)-D2*[zeros(4, 1) R]) <= 1e-6) - */ +*/
--- a/liboctave/numeric/CmplxGSVD.cc Thu Aug 04 07:50:31 2016 +0200 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,287 +0,0 @@ -// Copyright (C) 1996, 1997 John W. Eaton -// Copyright (C) 2006 Pascal Dupuis <Pascal.Dupuis@uclouvain.be> -// -// This program is free software; you can redistribute it and/or modify it under -// the terms of the GNU General Public License as published by the Free Software -// Foundation; either version 3 of the License, or (at your option) any later -// version. -// -// This program is distributed in the hope that it will be useful, but WITHOUT -// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or -// FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more -// details. -// -// You should have received a copy of the GNU General Public License along with -// this program; if not, see <http://www.gnu.org/licenses/>. - -#ifdef HAVE_CONFIG_H -# include <config.h> -#endif - -#include "CmplxGSVD.h" -#include "f77-fcn.h" -#include "lo-error.h" - - -/* - uncomment those lines to monitor k and l - #include "oct-obj.h" - #include "pager.h" -*/ - -extern "C" -{ - F77_RET_T - F77_FUNC (zggsvd, ZGGSVD) - ( - F77_CONST_CHAR_ARG_DECL, // JOBU (input) CHARACTER*1 - F77_CONST_CHAR_ARG_DECL, // JOBV (input) CHARACTER*1 - F77_CONST_CHAR_ARG_DECL, // JOBQ (input) CHARACTER*1 - const octave_idx_type&, // M (input) INTEGER - const octave_idx_type&, // N (input) INTEGER - const octave_idx_type&, // P (input) INTEGER - octave_idx_type &, // K (output) INTEGER - octave_idx_type &, // L (output) INTEGER - Complex*, // A (input/output) COMPLEX*16 array, dimension (LDA,N) - const octave_idx_type&, // LDA (input) INTEGER - Complex*, // B (input/output) COMPLEX*16 array, dimension (LDB,N) - const octave_idx_type&, // LDB (input) INTEGER - double*, // ALPHA (output) DOUBLE PRECISION array, dimension (N) - double*, // BETA (output) DOUBLE PRECISION array, dimension (N) - Complex*, // U (output) COMPLEX*16 array, dimension (LDU,M) - const octave_idx_type&, // LDU (input) INTEGER - Complex*, // V (output) COMPLEX*16 array, dimension (LDV,P) - const octave_idx_type&, // LDV (input) INTEGER - Complex*, // Q (output) COMPLEX*16 array, dimension (LDQ,N) - const octave_idx_type&, // LDQ (input) INTEGER - Complex*, // WORK (workspace) COMPLEX*16 array - double*, // RWORK (workspace) DOUBLE PRECISION array - int*, // IWORK (workspace/output) INTEGER array, dimension (N) - octave_idx_type& // INFO (output)INTEGER - F77_CHAR_ARG_LEN_DECL - F77_CHAR_ARG_LEN_DECL - F77_CHAR_ARG_LEN_DECL - ); -} - -ComplexMatrix -ComplexGSVD::left_singular_matrix_A (void) const -{ - if (type_computed == GSVD::sigma_only) - { - (*current_liboctave_error_handler) - ("CmplxGSVD: U not computed because type == GSVD::sigma_only"); - return ComplexMatrix (); - } - else - return left_smA; -} - -ComplexMatrix -ComplexGSVD::left_singular_matrix_B (void) const -{ - if (type_computed == GSVD::sigma_only) - { - (*current_liboctave_error_handler) - ("CmplxGSVD: V not computed because type == GSVD::sigma_only"); - return ComplexMatrix (); - } - else - return left_smB; -} - -ComplexMatrix -ComplexGSVD::right_singular_matrix (void) const -{ - if (type_computed == GSVD::sigma_only) - { - (*current_liboctave_error_handler) - ("CmplxGSVD: X not computed because type == GSVD::sigma_only"); - return ComplexMatrix (); - } - else - return right_sm; -} - -ComplexMatrix -ComplexGSVD::R_matrix (void) const -{ - if (type_computed != GSVD::std) - { - (*current_liboctave_error_handler) - ("CmplxGSVD: R not computed because type != GSVD::std"); - return ComplexMatrix (); - } - else - return R; -} - -octave_idx_type -ComplexGSVD::init (const ComplexMatrix& a, const ComplexMatrix& b, - GSVD::type gsvd_type) -{ - octave_idx_type info; - - octave_idx_type m = a.rows (); - octave_idx_type n = a.cols (); - octave_idx_type p = b.rows (); - - ComplexMatrix atmp = a; - Complex *tmp_dataA = atmp.fortran_vec (); - - ComplexMatrix btmp = b; - Complex *tmp_dataB = btmp.fortran_vec (); - - // octave_idx_type min_mn = m < n ? m : n; - - char jobu = 'U'; - char jobv = 'V'; - char jobq = 'Q'; - - octave_idx_type nrow_u = m; - octave_idx_type nrow_v = p; - octave_idx_type nrow_q = n; - - octave_idx_type k, l; - - switch (gsvd_type) - { - - case GSVD::sigma_only: - - // Note: for this case, both jobu and jobv should be 'N', but - // there seems to be a bug in dgesvd from Lapack V2.0. To - // demonstrate the bug, set both jobu and jobv to 'N' and find - // the singular values of [eye(3), eye(3)]. The result is - // [-sqrt(2), -sqrt(2), -sqrt(2)]. - // - // For Lapack 3.0, this problem seems to be fixed. - - jobu = 'N'; - jobv = 'N'; - jobq = 'N'; - nrow_u = nrow_v = nrow_q = 1; - break; - - default: - break; - } - - type_computed = gsvd_type; - - if (! (jobu == 'N' || jobu == 'O')) { - left_smA.resize (nrow_u, m); - } - - Complex *u = left_smA.fortran_vec (); - - if (! (jobv == 'N' || jobv == 'O')) { - left_smB.resize (nrow_v, p); - } - - Complex *v = left_smB.fortran_vec (); - - if (! (jobq == 'N' || jobq == 'O')) { - right_sm.resize (nrow_q, n); - } - Complex *q = right_sm.fortran_vec (); - - octave_idx_type lwork = 3*n; - lwork = lwork > m ? lwork : m; - lwork = (lwork > p ? lwork : p) + n; - - Array<Complex> work (dim_vector (lwork, 1)); - Array<double> alpha (dim_vector (n, 1)); - Array<double> beta (dim_vector (n, 1)); - Array<double> rwork(dim_vector (2*n, 1)); - Array<int> iwork (dim_vector (n, 1)); - - F77_XFCN (zggsvd, ZGGSVD, (F77_CONST_CHAR_ARG2 (&jobu, 1), - F77_CONST_CHAR_ARG2 (&jobv, 1), - F77_CONST_CHAR_ARG2 (&jobq, 1), - m, n, p, k, l, tmp_dataA, m, - tmp_dataB, p, alpha.fortran_vec (), - beta.fortran_vec (), u, nrow_u, - v, nrow_v, q, nrow_q, work.fortran_vec (), - rwork.fortran_vec (), iwork.fortran_vec (), - info - F77_CHAR_ARG_LEN (1) - F77_CHAR_ARG_LEN (1) - F77_CHAR_ARG_LEN (1))); - - if (f77_exception_encountered) - (*current_liboctave_error_handler) ("unrecoverable error in zggsvd"); - - if (info < 0) { - (*current_liboctave_error_handler) ("zggsvd.f: argument %d illegal", -info); - } else { - if (info > 0) { - (*current_liboctave_error_handler) ("zggsvd.f: Jacobi-type procedure failed to converge."); - } else { - octave_idx_type i, j; - - if (GSVD::std == gsvd_type) { - R.resize(k+l, k+l); - int astart = n-k-l; - if (m - k - l >= 0) { - int astart = n-k-l; - /* - * R is stored in A(1:K+L,N-K-L+1:N) - */ - for (i = 0; i < k+l; i++) - for (j = 0; j < k+l; j++) - R.xelem(i, j) = atmp.xelem(i, astart + j); - } else { - /* - * (R11 R12 R13 ) is stored in A(1:M, N-K-L+1:N), - * ( 0 R22 R23 ) - */ - - for (i = 0; i < m; i++) - for (j = 0; j < k+l; j++) - R.xelem(i, j) = atmp.xelem(i, astart + j); - /* - * and R33 is stored in B(M-K+1:L,N+M-K-L+1:N) - */ - for (i = k+l-1; i >=m; i--) { - for (j = 0; j < m; j++) - R.xelem(i, j) = 0.0; - for (j = m; j < k+l; j++) - R.xelem(i, j) = btmp.xelem(i - k, astart + j); - } - } - } - /* - uncomment this to monitor k and l - octave_value tmp; - octave_stdout << "CmplxGSVD k: "; - tmp = k; - tmp.print(octave_stdout); - octave_stdout << "\n"; - octave_stdout << "CmplxGSVD l: "; - tmp = l; - tmp.print(octave_stdout); - octave_stdout << "\n"; - */ - - if (m-k-l >= 0) { - // Fills in C and S - sigmaA.resize (l, l); - sigmaB.resize (l, l); - for (i = 0; i < l; i++) { - sigmaA.dgxelem(i) = alpha.elem(k+i); - sigmaB.dgxelem(i) = beta.elem(k+i); - } - } else { - // Fills in C and S - sigmaA.resize (m-k, m-k); - sigmaB.resize (m-k, m-k); - for (i = 0; i < m-k; i++) { - sigmaA.dgxelem(i) = alpha.elem(k+i); - sigmaB.dgxelem(i) = beta.elem(k+i); - } - } - } - } - return info; -}
--- a/liboctave/numeric/CmplxGSVD.h Thu Aug 04 07:50:31 2016 +0200 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,94 +0,0 @@ -// Copyright (C) 1996, 1997 John W. Eaton -// Copyright (C) 2006 Pascal Dupuis <Pascal.Dupuis@uclouvain.be> -// -// This program is free software; you can redistribute it and/or modify it under -// the terms of the GNU General Public License as published by the Free Software -// Foundation; either version 3 of the License, or (at your option) any later -// version. -// -// This program is distributed in the hope that it will be useful, but WITHOUT -// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or -// FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more -// details. -// -// You should have received a copy of the GNU General Public License along with -// this program; if not, see <http://www.gnu.org/licenses/>. - -#if !defined (octave_ComplexGSVD_h) -#define octave_ComplexGSVD_h 1 - -#include "octave-config.h" - -#include <iosfwd> - -#include "dDiagMatrix.h" -#include "CMatrix.h" -#include "dbleGSVD.h" - -class -ComplexGSVD -{ -public: - - ComplexGSVD (void) { } - - ComplexGSVD (const ComplexMatrix& a, const ComplexMatrix& b, - GSVD::type gsvd_type = GSVD::economy) - { - init (a, b, gsvd_type); - } - - ComplexGSVD (const ComplexMatrix& a, const ComplexMatrix& b, - octave_idx_type& info, GSVD::type gsvd_type = GSVD::economy) - { - info = init (a, b, gsvd_type); - } - - ComplexGSVD (const ComplexGSVD& a) - : type_computed (a.type_computed), - sigmaA (a.sigmaA), sigmaB (a.sigmaB), - left_smA (a.left_smA), left_smB (a.left_smB), right_sm (a.right_sm), - R(a.R) { } - - ComplexGSVD& operator = (const ComplexGSVD& a) - { - if (this != &a) - { - type_computed = a.type_computed; - sigmaA = a.sigmaA; - sigmaB = a.sigmaB; - left_smA = a.left_smA; - left_smB = a.left_smB; - right_sm = a.right_sm; - R = a.R; - } - - return *this; - } - - ~ComplexGSVD (void) { } - - DiagMatrix singular_values_A (void) const { return sigmaA; } - DiagMatrix singular_values_B (void) const { return sigmaB; } - - ComplexMatrix left_singular_matrix_A (void) const; - ComplexMatrix left_singular_matrix_B (void) const; - - ComplexMatrix right_singular_matrix (void) const; - ComplexMatrix R_matrix (void) const; - - friend std::ostream& operator << (std::ostream& os, const ComplexGSVD& a); - -private: - - GSVD::type type_computed; - - DiagMatrix sigmaA, sigmaB; - ComplexMatrix left_smA, left_smB; - ComplexMatrix right_sm, R; - - octave_idx_type init (const ComplexMatrix& a, const ComplexMatrix& b, - GSVD::type gsvd_type = GSVD::economy); -}; - -#endif
--- a/liboctave/numeric/dbleGSVD.cc Thu Aug 04 07:50:31 2016 +0200 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,293 +0,0 @@ -// Copyright (C) 1996, 1997 John W. Eaton -// Copyright (C) 2006 Pascal Dupuis <Pascal.Dupuis@uclouvain.be> -// -// This program is free software; you can redistribute it and/or modify it under -// the terms of the GNU General Public License as published by the Free Software -// Foundation; either version 3 of the License, or (at your option) any later -// version. -// -// This program is distributed in the hope that it will be useful, but WITHOUT -// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or -// FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more -// details. -// -// You should have received a copy of the GNU General Public License along with -// this program; if not, see <http://www.gnu.org/licenses/>. - -#ifdef HAVE_CONFIG_H -# include <config.h> -#endif - -#include "dbleGSVD.h" -#include "f77-fcn.h" -#include "lo-error.h" - -/* - uncomment those lines to monitor k and l - #include "oct-obj.h" - #include "pager.h" -*/ - -extern "C" -{ - F77_RET_T - F77_FUNC (dggsvd, DGGSVD) - ( - F77_CONST_CHAR_ARG_DECL, // JOBU (input) CHARACTER*1 - F77_CONST_CHAR_ARG_DECL, // JOBV (input) CHARACTER*1 - F77_CONST_CHAR_ARG_DECL, // JOBQ (input) CHARACTER*1 - const octave_idx_type&, // M (input) INTEGER - const octave_idx_type&, // N (input) INTEGER - const octave_idx_type&, // P (input) INTEGER - octave_idx_type &, // K (output) INTEGER - octave_idx_type &, // L (output) INTEGER - double*, // A (input/output) DOUBLE PRECISION array, dimension (LDA,N) - const octave_idx_type&, // LDA (input) INTEGER - double*, // B (input/output) DOUBLE PRECISION array, dimension (LDB,N) - const octave_idx_type&, // LDB (input) INTEGER - double*, // ALPHA (output) DOUBLE PRECISION array, dimension (N) - double*, // BETA (output) DOUBLE PRECISION array, dimension (N) - double*, // U (output) DOUBLE PRECISION array, dimension (LDU,M) - const octave_idx_type&, // LDU (input) INTEGER - double*, // V (output) DOUBLE PRECISION array, dimension (LDV,P) - const octave_idx_type&, // LDV (input) INTEGER - double*, // Q (output) DOUBLE PRECISION array, dimension (LDQ,N) - const octave_idx_type&, // LDQ (input) INTEGER - double*, // WORK (workspace) DOUBLE PRECISION array - int*, // IWORK (workspace/output) INTEGER array, dimension (N) - octave_idx_type& // INFO (output)INTEGER - F77_CHAR_ARG_LEN_DECL - F77_CHAR_ARG_LEN_DECL - F77_CHAR_ARG_LEN_DECL - ); -} - -Matrix -GSVD::left_singular_matrix_A (void) const -{ - if (type_computed == GSVD::sigma_only) - { - (*current_liboctave_error_handler) - ("dbleGSVD: U not computed because type == GSVD::sigma_only"); - return Matrix (); - } - else - return left_smA; -} - -Matrix -GSVD::left_singular_matrix_B (void) const -{ - if (type_computed == GSVD::sigma_only) - { - (*current_liboctave_error_handler) - ("dbleGSVD: V not computed because type == GSVD::sigma_only"); - return Matrix (); - } - else - return left_smB; -} - -Matrix -GSVD::right_singular_matrix (void) const -{ - if (type_computed == GSVD::sigma_only) - { - (*current_liboctave_error_handler) - ("dbleGSVD: X not computed because type == GSVD::sigma_only"); - return Matrix (); - } - else - return right_sm; -} -Matrix -GSVD::R_matrix (void) const -{ - if (type_computed != GSVD::std) - { - (*current_liboctave_error_handler) - ("dbleGSVD: R not computed because type != GSVD::std"); - return Matrix (); - } - else - return R; -} - -octave_idx_type -GSVD::init (const Matrix& a, const Matrix& b, GSVD::type gsvd_type) -{ - octave_idx_type info; - - octave_idx_type m = a.rows (); - octave_idx_type n = a.cols (); - octave_idx_type p = b.rows (); - - Matrix atmp = a; - double *tmp_dataA = atmp.fortran_vec (); - - Matrix btmp = b; - double *tmp_dataB = btmp.fortran_vec (); - - // octave_idx_type min_mn = m < n ? m : n; - - char jobu = 'U'; - char jobv = 'V'; - char jobq = 'Q'; - - octave_idx_type nrow_u = m; - octave_idx_type nrow_v = p; - octave_idx_type nrow_q = n; - - octave_idx_type k, l; - - switch (gsvd_type) - { - - case GSVD::sigma_only: - - // Note: for this case, both jobu and jobv should be 'N', but - // there seems to be a bug in dgesvd from Lapack V2.0. To - // demonstrate the bug, set both jobu and jobv to 'N' and find - // the singular values of [eye(3), eye(3)]. The result is - // [-sqrt(2), -sqrt(2), -sqrt(2)]. - // - // For Lapack 3.0, this problem seems to be fixed. - - jobu = 'N'; - jobv = 'N'; - jobq = 'N'; - nrow_u = nrow_v = nrow_q = 1; - break; - - default: - break; - } - - type_computed = gsvd_type; - - if (! (jobu == 'N' || jobu == 'O')) { - left_smA.resize (nrow_u, m); - } - - double *u = left_smA.fortran_vec (); - - if (! (jobv == 'N' || jobv == 'O')) { - left_smB.resize (nrow_v, p); - } - - double *v = left_smB.fortran_vec (); - - if (! (jobq == 'N' || jobq == 'O')) { - right_sm.resize (nrow_q, n); - } - double *q = right_sm.fortran_vec (); - - octave_idx_type lwork = 3*n; - lwork = lwork > m ? lwork : m; - lwork = (lwork > p ? lwork : p) + n; - - Array<double> work (dim_vector (lwork, 1)); - Array<double> alpha (dim_vector (n, 1)); - Array<double> beta (dim_vector (n, 1)); - Array<int> iwork (dim_vector (n, 1)); - - F77_XFCN (dggsvd, DGGSVD, (F77_CONST_CHAR_ARG2 (&jobu, 1), - F77_CONST_CHAR_ARG2 (&jobv, 1), - F77_CONST_CHAR_ARG2 (&jobq, 1), - m, n, p, k, l, tmp_dataA, m, - tmp_dataB, p, alpha.fortran_vec (), - beta.fortran_vec (), u, nrow_u, - v, nrow_v, q, nrow_q, work.fortran_vec (), - iwork.fortran_vec (), info - F77_CHAR_ARG_LEN (1) - F77_CHAR_ARG_LEN (1) - F77_CHAR_ARG_LEN (1))); - - if (f77_exception_encountered) - (*current_liboctave_error_handler) ("unrecoverable error in dggsvd"); - - if (info < 0) { - (*current_liboctave_error_handler) ("dggsvd.f: argument %d illegal", -info); - } else { - if (info > 0) { - (*current_liboctave_error_handler) ("dggsvd.f: Jacobi-type procedure failed to converge."); - } else { - octave_idx_type i, j; - - if (GSVD::std == gsvd_type) { - R.resize(k+l, k+l); - int astart = n-k-l; - if (m - k - l >= 0) { - int astart = n-k-l; - /* - * R is stored in A(1:K+L,N-K-L+1:N) - */ - for (i = 0; i < k+l; i++) - for (j = 0; j < k+l; j++) - R.xelem(i, j) = atmp.xelem(i, astart + j); - } else { - /* - * (R11 R12 R13 ) is stored in A(1:M, N-K-L+1:N), - * ( 0 R22 R23 ) - */ - - for (i = 0; i < m; i++) - for (j = 0; j < k+l; j++) - R.xelem(i, j) = atmp.xelem(i, astart + j); - /* - * and R33 is stored in B(M-K+1:L,N+M-K-L+1:N) - */ - for (i = k+l-1; i >=m; i--) { - for (j = 0; j < m; j++) - R.xelem(i, j) = 0.0; - for (j = m; j < k+l; j++) - R.xelem(i, j) = btmp.xelem(i - k, astart + j); - } - } - } - /* - uncomment this to monitor k and l - octave_value tmp; - octave_stdout << "dbleGSVD k: "; - tmp = k; - tmp.print(octave_stdout); - octave_stdout << "\n"; - octave_stdout << "dbleGSVD l: "; - tmp = l; - tmp.print(octave_stdout); - octave_stdout << "\n"; - */ - - if (m-k-l >= 0) { - // Fills in C and S - sigmaA.resize (l, l); - sigmaB.resize (l, l); - for (i = 0; i < l; i++) { - sigmaA.dgxelem(i) = alpha.elem(k+i); - sigmaB.dgxelem(i) = beta.elem(k+i); - } - } else { - // Fills in C and S - sigmaA.resize (m-k, m-k); - sigmaB.resize (m-k, m-k); - for (i = 0; i < m-k; i++) { - sigmaA.dgxelem(i) = alpha.elem(k+i); - sigmaB.dgxelem(i) = beta.elem(k+i); - } - } - } - } - return info; -} - -std::ostream& -operator << (std::ostream& os, const GSVD& a) -{ - os << a.left_singular_matrix_A () << "\n"; - os << a.left_singular_matrix_B () << "\n"; - os << a.singular_values_A () << "\n"; - os << a.singular_values_B () << "\n"; - os << a.right_singular_matrix () << "\n"; - - return os; -}
--- a/liboctave/numeric/dbleGSVD.h Thu Aug 04 07:50:31 2016 +0200 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,91 +0,0 @@ -// Copyright (C) 1996, 1997 John W. Eaton -// Copyright (C) 2006 Pascal Dupuis <Pascal.Dupuis@uclouvain.be> -// -// This program is free software; you can redistribute it and/or modify it under -// the terms of the GNU General Public License as published by the Free Software -// Foundation; either version 3 of the License, or (at your option) any later -// version. -// -// This program is distributed in the hope that it will be useful, but WITHOUT -// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or -// FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more -// details. -// -// You should have received a copy of the GNU General Public License along with -// this program; if not, see <http://www.gnu.org/licenses/>. - -#if !defined (octave_GSVD_h) -#define octave_GSVD_h 1 - -#include "octave-config.h" - -#include "dDiagMatrix.h" -#include "dMatrix.h" - -class -GSVD -{ -public: - - enum type - { - std, - economy, - sigma_only - }; - - GSVD (void) : sigmaA (), sigmaB (), left_smA (), left_smB (), right_sm () { } - GSVD (const Matrix& a, const Matrix& b, type gsvd_type = GSVD::economy) { init (a, b, gsvd_type); } - - GSVD (const Matrix& a, const Matrix& b, octave_idx_type& info, type gsvd_type = GSVD::economy) - { - info = init (a, b, gsvd_type); - } - - GSVD (const GSVD& a) - : type_computed (a.type_computed), - sigmaA (a.sigmaA), sigmaB (a.sigmaB), - left_smA (a.left_smA), left_smB (a.left_smB), right_sm (a.right_sm), - R(a.R) { } - - GSVD& operator = (const GSVD& a) - { - if (this != &a) - { - type_computed = a.type_computed; - sigmaA = a.sigmaA; - sigmaB = a.sigmaB; - left_smA = a.left_smA; - left_smB = a.left_smB; - right_sm = a.right_sm; - R = a.R; - } - - return *this; - } - - ~GSVD (void) { } - - DiagMatrix singular_values_A (void) const { return sigmaA; } - DiagMatrix singular_values_B (void) const { return sigmaB; } - - Matrix left_singular_matrix_A (void) const; - Matrix left_singular_matrix_B (void) const; - - Matrix right_singular_matrix (void) const; - Matrix R_matrix (void) const; - - friend std::ostream& operator << (std::ostream& os, const GSVD& a); - -private: - - GSVD::type type_computed; - - DiagMatrix sigmaA, sigmaB; - Matrix left_smA, left_smB; - Matrix right_sm, R; - - octave_idx_type init (const Matrix& a, const Matrix& b, type gsvd_type = economy); -}; - -#endif
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/numeric/gsvd.cc Tue Aug 09 18:02:11 2016 +0200 @@ -0,0 +1,357 @@ +// Copyright (C) 1996, 1997 John W. Eaton +// Copyright (C) 2006 Pascal Dupuis <Pascal.Dupuis@uclouvain.be> +// Copyright (C) 2016 Barbara Lócsi +// +// This program is free software; you can redistribute it and/or modify it under +// the terms of the GNU General Public License as published by the Free Software +// Foundation; either version 3 of the License, or (at your option) any later +// version. +// +// This program is distributed in the hope that it will be useful, but WITHOUT +// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +// FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more +// details. +// +// You should have received a copy of the GNU General Public License along with +// this program; if not, see <http://www.gnu.org/licenses/>. + +#ifdef HAVE_CONFIG_H +# include <config.h> +#endif + +#include "gsvd.h" +#include "f77-fcn.h" +#include "lo-error.h" +#include "CMatrix.h" +#include "dDiagMatrix.h" +#include "dMatrix.h" + +extern "C" +{ + F77_RET_T + F77_FUNC (dggsvd, DGGSVD) + ( + F77_CONST_CHAR_ARG_DECL, // JOBU (input) CHARACTER*1 + F77_CONST_CHAR_ARG_DECL, // JOBV (input) CHARACTER*1 + F77_CONST_CHAR_ARG_DECL, // JOBQ (input) CHARACTER*1 + const F77_INT&, // M (input) INTEGER + const F77_INT&, // N (input) INTEGER + const F77_INT&, // P (input) INTEGER + F77_INT &, // K (output) INTEGER + F77_INT &, // L (output) INTEGER + F77_DBLE*, // A (input/output) DOUBLE PRECISION array, dimension (LDA,N) + const F77_INT&, // LDA (input) INTEGER + F77_DBLE*, // B (input/output) DOUBLE PRECISION array, dimension (LDB,N) + const F77_INT&, // LDB (input) INTEGER + F77_DBLE*, // ALPHA (output) DOUBLE PRECISION array, dimension (N) + F77_DBLE*, // BETA (output) DOUBLE PRECISION array, dimension (N) + F77_DBLE*, // U (output) DOUBLE PRECISION array, dimension (LDU,M) + const F77_INT&, // LDU (input) INTEGER + F77_DBLE*, // V (output) DOUBLE PRECISION array, dimension (LDV,P) + const F77_INT&, // LDV (input) INTEGER + F77_DBLE*, // Q (output) DOUBLE PRECISION array, dimension (LDQ,N) + const F77_INT&, // LDQ (input) INTEGER + F77_DBLE*, // WORK (workspace) DOUBLE PRECISION array + int*, // IWORK (workspace/output) INTEGER array, dimension (N) + F77_INT& // INFO (output)INTEGER + F77_CHAR_ARG_LEN_DECL + F77_CHAR_ARG_LEN_DECL + F77_CHAR_ARG_LEN_DECL + ); + + F77_RET_T + F77_FUNC (zggsvd, ZGGSVD) + ( + F77_CONST_CHAR_ARG_DECL, // JOBU (input) CHARACTER*1 + F77_CONST_CHAR_ARG_DECL, // JOBV (input) CHARACTER*1 + F77_CONST_CHAR_ARG_DECL, // JOBQ (input) CHARACTER*1 + const F77_INT&, // M (input) INTEGER + const F77_INT&, // N (input) INTEGER + const F77_INT&, // P (input) INTEGER + F77_INT &, // K (output) INTEGER + F77_INT &, // L (output) INTEGER + F77_DBLE_CMPLX*, // A (input/output) COMPLEX*16 array, dimension (LDA,N) + const F77_INT&, // LDA (input) INTEGER + F77_DBLE_CMPLX*, // B (input/output) COMPLEX*16 array, dimension (LDB,N) + const F77_INT&, // LDB (input) INTEGER + F77_DBLE*, // ALPHA (output) DOUBLE PRECISION array, dimension (N) + F77_DBLE*, // BETA (output) DOUBLE PRECISION array, dimension (N) + F77_DBLE_CMPLX*, // U (output) COMPLEX*16 array, dimension (LDU,M) + const F77_INT&, // LDU (input) INTEGER + F77_DBLE_CMPLX*, // V (output) COMPLEX*16 array, dimension (LDV,P) + const F77_INT&, // LDV (input) INTEGER + F77_DBLE_CMPLX*, // Q (output) COMPLEX*16 array, dimension (LDQ,N) + const F77_INT&, // LDQ (input) INTEGER + F77_DBLE_CMPLX*, // WORK (workspace) COMPLEX*16 array + F77_DBLE*, // RWORK (workspace) DOUBLE PRECISION array + int*, // IWORK (workspace/output) INTEGER array, dimension (N) + F77_INT& // INFO (output)INTEGER + F77_CHAR_ARG_LEN_DECL + F77_CHAR_ARG_LEN_DECL + F77_CHAR_ARG_LEN_DECL + ); +} + +template <> +void +gsvd<Matrix>::ggsvd (char& jobu, char& jobv, char& jobq, octave_idx_type m, + octave_idx_type n, octave_idx_type p, octave_idx_type& k, + octave_idx_type& l, double *tmp_dataA, octave_idx_type m1, + double *tmp_dataB, octave_idx_type p1, Matrix& alpha, + Matrix& beta, double *u, octave_idx_type nrow_u, double *v, + octave_idx_type nrow_v, double *q, octave_idx_type nrow_q, + Matrix& work, octave_idx_type* iwork, octave_idx_type& info) +{ + F77_XFCN (dggsvd, DGGSVD, (F77_CONST_CHAR_ARG2 (&jobu, 1), + F77_CONST_CHAR_ARG2 (&jobv, 1), + F77_CONST_CHAR_ARG2 (&jobq, 1), + m, n, p, k, l, tmp_dataA, m1, + tmp_dataB, p1, alpha.fortran_vec (), + beta.fortran_vec (), u, nrow_u, + v, nrow_v, q, nrow_q, work.fortran_vec (), + iwork, info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); +} + + +template <> +void +gsvd<ComplexMatrix>::ggsvd (char& jobu, char& jobv, char& jobq, + octave_idx_type m, octave_idx_type n, + octave_idx_type p, octave_idx_type& k, + octave_idx_type& l, Complex *tmp_dataA, + octave_idx_type m1, Complex *tmp_dataB, + octave_idx_type p1, Matrix& alpha, + Matrix& beta, Complex *u, octave_idx_type nrow_u, + Complex *v, octave_idx_type nrow_v, Complex *q, + octave_idx_type nrow_q, ComplexMatrix& work, + octave_idx_type* iwork, octave_idx_type& info) +{ + Matrix rwork(2*n, 1); + F77_XFCN (zggsvd, ZGGSVD, (F77_CONST_CHAR_ARG2 (&jobu, 1), + F77_CONST_CHAR_ARG2 (&jobv, 1), + F77_CONST_CHAR_ARG2 (&jobq, 1), + m, n, p, k, l, F77_DBLE_CMPLX_ARG (tmp_dataA), m1, + F77_DBLE_CMPLX_ARG (tmp_dataB), p1, + alpha.fortran_vec (), beta.fortran_vec (), + F77_DBLE_CMPLX_ARG (u), nrow_u, + F77_DBLE_CMPLX_ARG (v), nrow_v, + F77_DBLE_CMPLX_ARG (q), nrow_q, + F77_DBLE_CMPLX_ARG (work.fortran_vec ()), + rwork.fortran_vec (), iwork, info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); +} + +template <typename T> +T +gsvd<T>::left_singular_matrix_A (void) const +{ + if (type == gsvd::Type::sigma_only) + { + (*current_liboctave_error_handler) + ("gsvd: U not computed because type == gsvd::sigma_only"); + return T (); + } + else + return left_smA; +} + +template <typename T> +T +gsvd<T>::left_singular_matrix_B (void) const +{ + if (type == gsvd::Type::sigma_only) + { + (*current_liboctave_error_handler) + ("gsvd: V not computed because type == gsvd::sigma_only"); + return T (); + } + else + return left_smB; +} + +template <typename T> +T +gsvd<T>::right_singular_matrix (void) const +{ + if (type == gsvd::Type::sigma_only) + { + (*current_liboctave_error_handler) + ("gsvd: X not computed because type == gsvd::sigma_only"); + return T (); + } + else + return right_sm; +} + +template <typename T> +T +gsvd<T>::R_matrix (void) const +{ + if (type != gsvd::Type::std) + { + (*current_liboctave_error_handler) + ("gsvd: R not computed because type != gsvd::std"); + return T (); + } + else + return R; +} + +template <typename T> +gsvd<T>::gsvd (const T& a, const T& b, gsvd::Type gsvd_type) +{ + octave_idx_type info; + + octave_idx_type m = a.rows (); + octave_idx_type n = a.cols (); + octave_idx_type p = b.rows (); + + T atmp = a; + P *tmp_dataA = atmp.fortran_vec (); + + T btmp = b; + P *tmp_dataB = btmp.fortran_vec (); + + char jobu = 'U'; + char jobv = 'V'; + char jobq = 'Q'; + + octave_idx_type nrow_u = m; + octave_idx_type nrow_v = p; + octave_idx_type nrow_q = n; + + octave_idx_type k, l; + + switch (gsvd_type) + { + case gsvd<T>::Type::sigma_only: + + // Note: for this case, both jobu and jobv should be 'N', but + // there seems to be a bug in dgesvd from Lapack V2.0. To + // demonstrate the bug, set both jobu and jobv to 'N' and find + // the singular values of [eye(3), eye(3)]. The result is + // [-sqrt(2), -sqrt(2), -sqrt(2)]. + // + // For Lapack 3.0, this problem seems to be fixed. + + jobu = 'N'; + jobv = 'N'; + jobq = 'N'; + nrow_u = nrow_v = nrow_q = 1; + break; + + default: + break; + } + + type = gsvd_type; + + if (! (jobu == 'N' || jobu == 'O')) + left_smA.resize (nrow_u, m); + + P *u = left_smA.fortran_vec (); + + if (! (jobv == 'N' || jobv == 'O')) + left_smB.resize (nrow_v, p); + + P *v = left_smB.fortran_vec (); + + if (! (jobq == 'N' || jobq == 'O')) + right_sm.resize (nrow_q, n); + + P *q = right_sm.fortran_vec (); + + octave_idx_type lwork = 3*n; + lwork = lwork > m ? lwork : m; + lwork = (lwork > p ? lwork : p) + n; + + T work (lwork, 1); + Matrix alpha (n, 1); + Matrix beta (n, 1); + + std::vector<octave_idx_type> iwork (n); + + gsvd<T>::ggsvd (jobu, jobv, jobq, m, n, p, k, l, + tmp_dataA, m, tmp_dataB, p, alpha, beta, u, + nrow_u, v, nrow_v, q, nrow_q, work, iwork.data (), info); + + if (f77_exception_encountered) + (*current_liboctave_error_handler) ("unrecoverable error in *ggsvd"); + + if (info < 0) + (*current_liboctave_error_handler) ("*ggsvd.f: argument %d illegal", -info); + else + { + if (info > 0) + (*current_liboctave_error_handler) + ("*ggsvd.f: Jacobi-type procedure failed to converge."); + else + { + octave_idx_type i, j; + + if (gsvd<T>::Type::std == gsvd_type) + { + R.resize(k+l, k+l); + int astart = n-k-l; + if (m - k - l >= 0) + { + astart = n-k-l; + // R is stored in A(1:K+L,N-K-L+1:N) + for (i = 0; i < k+l; i++) + for (j = 0; j < k+l; j++) + R.xelem (i, j) = atmp.xelem (i, astart + j); + } + else + { + // (R11 R12 R13 ) is stored in A(1:M, N-K-L+1:N), + // ( 0 R22 R23 ) + + for (i = 0; i < m; i++) + for (j = 0; j < k+l; j++) + R.xelem (i, j) = atmp.xelem (i, astart + j); + // and R33 is stored in B(M-K+1:L,N+M-K-L+1:N) + for (i = k+l-1; i >=m; i--) { + for (j = 0; j < m; j++) + R.xelem(i, j) = 0.0; + for (j = m; j < k+l; j++) + R.xelem (i, j) = btmp.xelem (i - k, astart + j); + } + } + } + + if (m-k-l >= 0) + { + // Fills in C and S + sigmaA.resize (l, l); + sigmaB.resize (l, l); + for (i = 0; i < l; i++) + { + sigmaA.dgxelem(i) = alpha.elem(k+i); + sigmaB.dgxelem(i) = beta.elem(k+i); + } + } + else + { + // Fills in C and S + sigmaA.resize (m-k, m-k); + sigmaB.resize (m-k, m-k); + for (i = 0; i < m-k; i++) + { + sigmaA.dgxelem(i) = alpha.elem(k+i); + sigmaB.dgxelem(i) = beta.elem(k+i); + } + } + } + } +} + +// Instantiations we need. + +template class gsvd<Matrix>; + +template class gsvd<ComplexMatrix>;
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/numeric/gsvd.h Tue Aug 09 18:02:11 2016 +0200 @@ -0,0 +1,95 @@ +// Copyright (C) 1996, 1997 John W. Eaton +// Copyright (C) 2006 Pascal Dupuis <Pascal.Dupuis@uclouvain.be> +// Copyright (C) 2016 Barbara Lócsi +// +// This program is free software; you can redistribute it and/or modify it under +// the terms of the GNU General Public License as published by the Free Software +// Foundation; either version 3 of the License, or (at your option) any later +// version. +// +// This program is distributed in the hope that it will be useful, but WITHOUT +// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +// FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more +// details. +// +// You should have received a copy of the GNU General Public License along with +// this program; if not, see <http://www.gnu.org/licenses/>. + +#if !defined (octave_gsvd_h) +#define octave_gsvd_h 1 + +#include "octave-config.h" + +#include "dDiagMatrix.h" +#include "dMatrix.h" + +template <typename T> +class +gsvd +{ +public: + + enum class Type + { + std, + economy, + sigma_only + }; + + gsvd (void) : sigmaA (), sigmaB (), left_smA (), left_smB (), right_sm () { } + + gsvd (const T& a, const T& b, gsvd::Type gsvd_type = gsvd<T>::Type::economy); + + gsvd (const gsvd& a) + : type (a.type), + sigmaA (a.sigmaA), sigmaB (a.sigmaB), + left_smA (a.left_smA), left_smB (a.left_smB), right_sm (a.right_sm), + R(a.R) { } + + gsvd& operator = (const gsvd& a) + { + if (this != &a) + { + type = a.type; + sigmaA = a.sigmaA; + sigmaB = a.sigmaB; + left_smA = a.left_smA; + left_smB = a.left_smB; + right_sm = a.right_sm; + R = a.R; + } + + return *this; + } + + ~gsvd (void) { } + + DiagMatrix singular_values_A (void) const { return sigmaA; } + DiagMatrix singular_values_B (void) const { return sigmaB; } + + T left_singular_matrix_A (void) const; + T left_singular_matrix_B (void) const; + + T right_singular_matrix (void) const; + T R_matrix (void) const; + +private: + + gsvd::Type type; + + typedef typename T::element_type P; + + DiagMatrix sigmaA, sigmaB; + T left_smA, left_smB; + T right_sm, R; + + void ggsvd (char& jobu, char& jobv, char& jobq, octave_idx_type m, + octave_idx_type n, octave_idx_type p, octave_idx_type& k, + octave_idx_type& l, P *tmp_dataA, octave_idx_type m1, + P *tmp_dataB, octave_idx_type p1, Matrix& alpha, Matrix& beta, + P *u, octave_idx_type nrow_u, P *v, octave_idx_type nrow_v, P *q, + octave_idx_type nrow_q, T& work, octave_idx_type* iwork, + octave_idx_type& info); +}; + +#endif
--- a/liboctave/numeric/module.mk Thu Aug 04 07:50:31 2016 +0200 +++ b/liboctave/numeric/module.mk Tue Aug 09 18:02:11 2016 +0200 @@ -8,7 +8,6 @@ LIBOCTAVE_OPT_IN = $(LIBOCTAVE_OPT_INC:.h=.in) NUMERIC_INC = \ - liboctave/numeric/CmplxGSVD.h \ liboctave/numeric/CollocWt.h \ liboctave/numeric/DAE.h \ liboctave/numeric/DAEFunc.h \ @@ -17,9 +16,9 @@ liboctave/numeric/DASPK.h \ liboctave/numeric/DASRT.h \ liboctave/numeric/DASSL.h \ - liboctave/numeric/dbleGSVD.h \ liboctave/numeric/DET.h \ liboctave/numeric/EIG.h \ + liboctave/numeric/gsvd.h \ liboctave/numeric/LSODE.h \ liboctave/numeric/ODE.h \ liboctave/numeric/ODEFunc.h \ @@ -58,13 +57,12 @@ liboctave/numeric/svd.h NUMERIC_SRC = \ - liboctave/numeric/CmplxGSVD.cc \ liboctave/numeric/CollocWt.cc \ liboctave/numeric/DASPK.cc \ liboctave/numeric/DASRT.cc \ liboctave/numeric/DASSL.cc \ - liboctave/numeric/dbleGSVD.cc \ liboctave/numeric/EIG.cc \ + liboctave/numeric/gsvd.cc \ liboctave/numeric/LSODE.cc \ liboctave/numeric/ODES.cc \ liboctave/numeric/Quad.cc \
--- a/liboctave/operators/mx-defs.h Thu Aug 04 07:50:31 2016 +0200 +++ b/liboctave/operators/mx-defs.h Tue Aug 09 18:02:11 2016 +0200 @@ -66,7 +66,7 @@ class EIG; -class GSVD; +template <typename T> class gsvd; template <typename T> class hess;
--- a/liboctave/operators/mx-ext.h Thu Aug 04 07:50:31 2016 +0200 +++ b/liboctave/operators/mx-ext.h Tue Aug 09 18:02:11 2016 +0200 @@ -59,8 +59,7 @@ // Result of a Generalized Singular Value Decomposition. -#include "dbleGSVD.h" -#include "CmplxGSVD.h" +#include "gsvd.h" // Result of an LU decomposition.