# HG changeset patch # User Carlo de Falco # Date 1488739737 -3600 # Sun Mar 05 19:48:57 2017 +0100 # Node ID ca9f21b242e68bad07d1dbb552c29c3de4d3388d # Parent b506b43f999e5f625d8b1f0ea304b8c06090a978 [mq]: cggsvd diff --git a/liboctave/numeric/gsvd.cc b/liboctave/numeric/gsvd.cc --- a/liboctave/numeric/gsvd.cc +++ b/liboctave/numeric/gsvd.cc @@ -50,20 +50,20 @@ double *tmp_dataA, F77_INT m1, double *tmp_dataB, F77_INT p1, Matrix& alpha, Matrix& beta, double *u, F77_INT nrow_u, double *v, F77_INT nrow_v, double *q, - F77_INT nrow_q, Matrix& work, F77_INT* iwork, - F77_INT& info) + F77_INT nrow_q, Matrix& work, F77_INT lwork, + F77_INT* iwork, F77_INT& 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))); + F77_XFCN (dggsvd3, DGGSVD3, (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 (), + lwork, iwork, info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); } template <> @@ -75,19 +75,19 @@ FloatMatrix& beta, float *u, F77_INT nrow_u, float *v, F77_INT nrow_v, float *q, F77_INT nrow_q, FloatMatrix& work, - F77_INT* iwork, F77_INT& info) + F77_INT lwork, F77_INT* iwork, F77_INT& 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))); + F77_XFCN (sggsvd3, SGGSVD3, (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 (), + lwork, iwork, info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); } template <> @@ -99,23 +99,23 @@ Matrix& beta, Complex *u, F77_INT nrow_u, Complex *v, F77_INT nrow_v, Complex *q, F77_INT nrow_q, ComplexMatrix& work, - F77_INT* iwork, F77_INT& info) + F77_INT lwork, F77_INT* iwork, F77_INT& 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))); + F77_XFCN (zggsvd3, ZGGSVD3, (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 ()), + lwork, rwork.fortran_vec (), iwork, info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); } template <> @@ -129,24 +129,24 @@ FloatComplex *u, F77_INT nrow_u, FloatComplex *v, F77_INT nrow_v, FloatComplex *q, F77_INT nrow_q, - FloatComplexMatrix& work, + FloatComplexMatrix& work, F77_INT lwork, F77_INT* iwork, F77_INT& info) { 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))); + F77_XFCN (cggsvd3, CGGSVD3, (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 ()), lwork, + rwork.fortran_vec (), iwork, info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); } template @@ -269,20 +269,28 @@ P *q = right_sm.fortran_vec (); - F77_INT lwork = 3*n; - lwork = lwork > m ? lwork : m; - lwork = (lwork > p ? lwork : p) + n; - - T work (lwork, 1); real_matrix alpha (n, 1); real_matrix beta (n, 1); std::vector iwork (n); + F77_INT lwork = -1; + T work_tmp (1, 1); + gsvd::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); + nrow_u, v, nrow_v, q, nrow_q, work_tmp, lwork, + iwork.data (), info); + lwork = static_cast (std::abs (work_tmp(1, 1))); + T work (lwork, 1); + info = 0; + + gsvd::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, lwork, + iwork.data (), info); + if (f77_exception_encountered) (*current_liboctave_error_handler) ("unrecoverable error in *ggsvd"); diff --git a/liboctave/numeric/gsvd.h b/liboctave/numeric/gsvd.h --- a/liboctave/numeric/gsvd.h +++ b/liboctave/numeric/gsvd.h @@ -104,7 +104,8 @@ P *u, octave_f77_int_type nrow_u, P *v, octave_f77_int_type nrow_v, P *q, octave_f77_int_type nrow_q, - T& work, octave_f77_int_type* iwork, + T& work, octave_f77_int_type lwork, + octave_f77_int_type* iwork, octave_f77_int_type& info); }; } diff --git a/liboctave/numeric/lo-lapack-proto.h b/liboctave/numeric/lo-lapack-proto.h --- a/liboctave/numeric/lo-lapack-proto.h +++ b/liboctave/numeric/lo-lapack-proto.h @@ -810,7 +810,7 @@ // GGSVD F77_RET_T - F77_FUNC (dggsvd, DGGSVD) + F77_FUNC (dggsvd3, DGGSVD3) (F77_CONST_CHAR_ARG_DECL, // JOBU F77_CONST_CHAR_ARG_DECL, // JOBV F77_CONST_CHAR_ARG_DECL, // JOBQ @@ -832,6 +832,7 @@ F77_DBLE*, // Q(LDQ,N) const F77_INT&, // LDQ F77_DBLE*, // WORK + const F77_INT&, // LWORK F77_INT*, // IWORK(N) F77_INT& // INFO F77_CHAR_ARG_LEN_DECL @@ -839,8 +840,8 @@ F77_CHAR_ARG_LEN_DECL); F77_RET_T - F77_FUNC (sggsvd, SGGSVD) - (F77_CONST_CHAR_ARG_DECL, // JOBU + F77_FUNC (sggsvd3, SGGSVD3) + (F77_CONST_CHAR_ARG_DECL, // JOBU F77_CONST_CHAR_ARG_DECL, // JOBV F77_CONST_CHAR_ARG_DECL, // JOBQ const F77_INT&, // M @@ -861,6 +862,7 @@ F77_REAL*, // Q const F77_INT&, // LDQ F77_REAL*, // WORK + const F77_INT&, // LWORK F77_INT*, // IWORK F77_INT& // INFO F77_CHAR_ARG_LEN_DECL @@ -868,7 +870,7 @@ F77_CHAR_ARG_LEN_DECL); F77_RET_T - F77_FUNC (zggsvd, ZGGSVD) + F77_FUNC (zggsvd3, ZGGSVD3) (F77_CONST_CHAR_ARG_DECL, // JOBU F77_CONST_CHAR_ARG_DECL, // JOBV F77_CONST_CHAR_ARG_DECL, // JOBQ @@ -890,6 +892,7 @@ F77_DBLE_CMPLX*, // Q(LDQ,N) const F77_INT&, // LDQ F77_DBLE_CMPLX*, // WORK + const F77_INT&, // LWORK F77_DBLE*, // RWORK F77_INT*, // IWORK(N) F77_INT& // INFO @@ -898,7 +901,7 @@ F77_CHAR_ARG_LEN_DECL); F77_RET_T - F77_FUNC (cggsvd, CGGSVD) + F77_FUNC (cggsvd3, CGGSVD3) (F77_CONST_CHAR_ARG_DECL, // JOBU F77_CONST_CHAR_ARG_DECL, // JOBV F77_CONST_CHAR_ARG_DECL, // JOBQ @@ -920,6 +923,7 @@ F77_CMPLX*, // Q const F77_INT&, // LDQ F77_CMPLX*, // WORK + const F77_INT&, // LWORK F77_REAL*, // RWORK F77_INT*, // IWORK F77_INT& // INFO