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,21 @@ 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 <> @@ -74,20 +75,20 @@ F77_INT p1, FloatMatrix& alpha, 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 nrow_q, FloatMatrix& work, + 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 +100,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 +130,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,19 +270,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(0, 0))); + 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,8 +901,8 @@ F77_CHAR_ARG_LEN_DECL); F77_RET_T - F77_FUNC (cggsvd, CGGSVD) - (F77_CONST_CHAR_ARG_DECL, // JOBU + F77_FUNC (cggsvd3, CGGSVD3) + (F77_CONST_CHAR_ARG_DECL, // JOBU F77_CONST_CHAR_ARG_DECL, // JOBV F77_CONST_CHAR_ARG_DECL, // JOBQ const F77_INT&, // M @@ -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