Mercurial > octave
changeset 22437:0aee8b620864
gsvd: move new class into the octave::math namespace.
| author | Carnë Draug <carandraug@octave.org> |
|---|---|
| date | Sun, 04 Sep 2016 17:37:35 +0100 |
| parents | 09005ac7d56c |
| children | 5a838a892adb |
| files | libinterp/corefcn/gsvd.cc liboctave/numeric/gsvd.cc liboctave/numeric/gsvd.h |
| diffstat | 3 files changed, 366 insertions(+), 345 deletions(-) [+] |
line wrap: on
line diff
--- a/libinterp/corefcn/gsvd.cc Sun Sep 04 17:12:03 2016 +0100 +++ b/libinterp/corefcn/gsvd.cc Sun Sep 04 17:37:35 2016 +0100 @@ -34,21 +34,21 @@ template <typename T> -static typename gsvd<T>::Type +static typename octave::math::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); + ? octave::math::gsvd<T>::Type::sigma_only + : (nargout > 5) ? octave::math::gsvd<T>::Type::std + : octave::math::gsvd<T>::Type::economy); } - // Named like this to avoid conflicts with the gsvd class. template <typename T> static octave_value_list function_gsvd (const T& A, const T& B, const octave_idx_type nargout) { - gsvd<T> result (A, B, gsvd_type<T> (nargout)); + octave::math::gsvd<T> result (A, B, gsvd_type<T> (nargout)); octave_value_list retval (nargout); if (nargout < 2)
--- a/liboctave/numeric/gsvd.cc Sun Sep 04 17:12:03 2016 +0100 +++ b/liboctave/numeric/gsvd.cc Sun Sep 04 17:37:35 2016 +0100 @@ -32,326 +32,340 @@ #include "dDiagMatrix.h" #include "fDiagMatrix.h" - -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<FloatMatrix>::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, - float *tmp_dataA, octave_idx_type m1, - float *tmp_dataB, octave_idx_type p1, - FloatMatrix& alpha, FloatMatrix& beta, float *u, - octave_idx_type nrow_u, float *v, - octave_idx_type nrow_v, float *q, - octave_idx_type nrow_q, FloatMatrix& work, - octave_idx_type* iwork, octave_idx_type& info) -{ - F77_XFCN (sggsvd, SGGSVD, (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 <> -void -gsvd<FloatComplexMatrix>::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, FloatComplex *tmp_dataA, - octave_idx_type m1, FloatComplex *tmp_dataB, - octave_idx_type p1, FloatMatrix& alpha, - FloatMatrix& beta, FloatComplex *u, - octave_idx_type nrow_u, FloatComplex *v, - octave_idx_type nrow_v, FloatComplex *q, - octave_idx_type nrow_q, - FloatComplexMatrix& work, - octave_idx_type* iwork, octave_idx_type& info) +namespace octave { - FloatMatrix rwork(2*n, 1); - F77_XFCN (cggsvd, CGGSVD, (F77_CONST_CHAR_ARG2 (&jobu, 1), - F77_CONST_CHAR_ARG2 (&jobv, 1), - F77_CONST_CHAR_ARG2 (&jobq, 1), - m, n, p, k, l, F77_CMPLX_ARG (tmp_dataA), m1, - F77_CMPLX_ARG (tmp_dataB), p1, - alpha.fortran_vec (), beta.fortran_vec (), - F77_CMPLX_ARG (u), nrow_u, - F77_CMPLX_ARG (v), nrow_v, - F77_CMPLX_ARG (q), nrow_q, - F77_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) + namespace math + { + 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) { - (*current_liboctave_error_handler) - ("gsvd: U not computed because type == gsvd::sigma_only"); - return T (); + 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))); } - else - return left_smA; -} -template <typename T> -T -gsvd<T>::left_singular_matrix_B (void) const -{ - if (type == gsvd::Type::sigma_only) + template <> + void + gsvd<FloatMatrix>::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, + float *tmp_dataA, octave_idx_type m1, + float *tmp_dataB, octave_idx_type p1, + FloatMatrix& alpha, FloatMatrix& beta, float *u, + octave_idx_type nrow_u, float *v, + octave_idx_type nrow_v, float *q, + octave_idx_type nrow_q, FloatMatrix& work, + octave_idx_type* iwork, octave_idx_type& info) { - (*current_liboctave_error_handler) - ("gsvd: V not computed because type == gsvd::sigma_only"); - return T (); + F77_XFCN (sggsvd, SGGSVD, (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))); } - 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) + 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) { - (*current_liboctave_error_handler) - ("gsvd: R not computed because type != gsvd::std"); - return T (); + 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))); } - 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) + template <> + void + gsvd<FloatComplexMatrix>::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, + FloatComplex *tmp_dataA, + octave_idx_type m1, + FloatComplex *tmp_dataB, + octave_idx_type p1, FloatMatrix& alpha, + FloatMatrix& beta, FloatComplex *u, + octave_idx_type nrow_u, FloatComplex *v, + octave_idx_type nrow_v, FloatComplex *q, + octave_idx_type nrow_q, + FloatComplexMatrix& work, + octave_idx_type* iwork, + octave_idx_type& info) { - case gsvd<T>::Type::sigma_only: + FloatMatrix rwork(2*n, 1); + F77_XFCN (cggsvd, CGGSVD, (F77_CONST_CHAR_ARG2 (&jobu, 1), + F77_CONST_CHAR_ARG2 (&jobv, 1), + F77_CONST_CHAR_ARG2 (&jobq, 1), + m, n, p, k, l, F77_CMPLX_ARG (tmp_dataA), m1, + F77_CMPLX_ARG (tmp_dataB), p1, + alpha.fortran_vec (), beta.fortran_vec (), + F77_CMPLX_ARG (u), nrow_u, + F77_CMPLX_ARG (v), nrow_v, + F77_CMPLX_ARG (q), nrow_q, + F77_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))); + } - // 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. + 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; + } - jobu = 'N'; - jobv = 'N'; - jobq = 'N'; - nrow_u = nrow_v = nrow_q = 1; - break; - - default: - break; + 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; } - 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); + 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; + } - P *v = left_smB.fortran_vec (); + 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; + } - if (! (jobq == 'N' || jobq == 'O')) - right_sm.resize (nrow_q, n); + template <typename T> + gsvd<T>::gsvd (const T& a, const T& b, gsvd::Type gsvd_type) + { + octave_idx_type info; - P *q = right_sm.fortran_vec (); + 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 (); - octave_idx_type lwork = 3*n; - lwork = lwork > m ? lwork : m; - lwork = (lwork > p ? lwork : p) + 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; - T work (lwork, 1); - real_matrix alpha (n, 1); - real_matrix beta (n, 1); + 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. - std::vector<octave_idx_type> iwork (n); + jobu = 'N'; + jobv = 'N'; + jobq = 'N'; + nrow_u = nrow_v = nrow_q = 1; + break; - 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); + default: + break; + } + + type = gsvd_type; + + if (! (jobu == 'N' || jobu == 'O')) + left_smA.resize (nrow_u, m); + + P *u = left_smA.fortran_vec (); - if (f77_exception_encountered) - (*current_liboctave_error_handler) ("unrecoverable error in *ggsvd"); + 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; - 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."); + T work (lwork, 1); + real_matrix alpha (n, 1); + real_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 { - octave_idx_type i, j; - - if (gsvd<T>::Type::std == gsvd_type) + if (info > 0) + (*current_liboctave_error_handler) + ("*ggsvd.f: Jacobi-type procedure failed to converge."); + else { - R.resize(k+l, k+l); - int astart = n-k-l; - if (m - k - l >= 0) + octave_idx_type i, j; + + if (gsvd<T>::Type::std == gsvd_type) { - 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); + 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 { - // (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--) + // Fills in C and S + sigmaA.resize (m-k, m-k); + sigmaB.resize (m-k, m-k); + for (i = 0; i < m-k; 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); + sigmaA.dgxelem(i) = alpha.elem(k+i); + sigmaB.dgxelem(i) = beta.elem(k+i); } } } - - 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<FloatMatrix>; + template class gsvd<ComplexMatrix>; + template class gsvd<FloatComplexMatrix>; + } } - -// Instantiations we need. -template class gsvd<Matrix>; -template class gsvd<FloatMatrix>; -template class gsvd<ComplexMatrix>; -template class gsvd<FloatComplexMatrix>;
--- a/liboctave/numeric/gsvd.h Sun Sep 04 17:12:03 2016 +0100 +++ b/liboctave/numeric/gsvd.h Sun Sep 04 17:37:35 2016 +0100 @@ -20,76 +20,83 @@ #include "octave-config.h" -template <typename T> -class -gsvd +namespace octave { -public: - - enum class Type + namespace math { - std, - economy, - sigma_only - }; + template <typename T> + class + gsvd + { + public: - gsvd (void) : sigmaA (), sigmaB (), left_smA (), left_smB (), right_sm () { } + enum class Type + { + std, + economy, + sigma_only + }; - gsvd (const T& a, const T& b, gsvd::Type gsvd_type = gsvd<T>::Type::economy); + gsvd (void) : sigmaA (), sigmaB (), left_smA (), left_smB (), right_sm () + { } - 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 (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) + gsvd& operator = (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; + 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; } - return *this; - } + ~gsvd (void) { } - ~gsvd (void) { } - - typename T::real_diag_matrix_type - singular_values_A (void) const { return sigmaA; } + typename T::real_diag_matrix_type + singular_values_A (void) const { return sigmaA; } - typename T::real_diag_matrix_type - singular_values_B (void) const { return sigmaB; } + typename T::real_diag_matrix_type + singular_values_B (void) const { return sigmaB; } - T left_singular_matrix_A (void) const; - T left_singular_matrix_B (void) const; + 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; + T right_singular_matrix (void) const; + T R_matrix (void) const; -private: - typedef typename T::value_type P; - typedef typename T::real_matrix_type real_matrix; + private: + typedef typename T::value_type P; + typedef typename T::real_matrix_type real_matrix; - gsvd::Type type; - typename T::real_diag_matrix_type sigmaA, sigmaB; - T left_smA, left_smB; - T right_sm, R; + gsvd::Type type; + typename T::real_diag_matrix_type 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, real_matrix& alpha, - real_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); -}; + 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, real_matrix& alpha, + real_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 -