Mercurial > octave
changeset 9661:afcf852256d2
optimize / and '\ for triangular matrices
| author | Jaroslav Hajek <highegg@gmail.com> |
|---|---|
| date | Wed, 23 Sep 2009 10:00:16 +0200 |
| parents | 0256e187d13b |
| children | 0d3b248f4ab6 |
| files | liboctave/CMatrix.cc liboctave/CMatrix.h liboctave/ChangeLog liboctave/dMatrix.cc liboctave/dMatrix.h liboctave/fCMatrix.cc liboctave/fCMatrix.h liboctave/fMatrix.cc liboctave/fMatrix.h liboctave/mx-defs.h src/ChangeLog src/OPERATORS/op-cm-cm.cc src/OPERATORS/op-fcm-fcm.cc src/OPERATORS/op-fm-fm.cc src/OPERATORS/op-m-m.cc src/ov.cc src/ov.h src/pt-cbinop.cc src/xdiv.cc src/xdiv.h |
| diffstat | 20 files changed, 389 insertions(+), 227 deletions(-) [+] |
line wrap: on
line diff
--- a/liboctave/CMatrix.cc Mon Sep 21 14:24:27 2009 +0200 +++ b/liboctave/CMatrix.cc Wed Sep 23 10:00:16 2009 +0200 @@ -1852,7 +1852,7 @@ ComplexMatrix::utsolve (MatrixType &mattype, const ComplexMatrix& b, octave_idx_type& info, double& rcon, solve_singularity_handler sing_handler, - bool calc_cond) const + bool calc_cond, blas_trans_type transt) const { ComplexMatrix retval; @@ -1928,7 +1928,7 @@ Complex *result = retval.fortran_vec (); char uplo = 'U'; - char trans = 'N'; + char trans = get_blas_char (transt); char dia = 'N'; F77_XFCN (ztrtrs, ZTRTRS, (F77_CONST_CHAR_ARG2 (&uplo, 1), @@ -1953,7 +1953,7 @@ ComplexMatrix::ltsolve (MatrixType &mattype, const ComplexMatrix& b, octave_idx_type& info, double& rcon, solve_singularity_handler sing_handler, - bool calc_cond) const + bool calc_cond, blas_trans_type transt) const { ComplexMatrix retval; @@ -2029,7 +2029,7 @@ Complex *result = retval.fortran_vec (); char uplo = 'L'; - char trans = 'N'; + char trans = get_blas_char (transt); char dia = 'N'; F77_XFCN (ztrtrs, ZTRTRS, (F77_CONST_CHAR_ARG2 (&uplo, 1), @@ -2260,10 +2260,10 @@ ComplexMatrix ComplexMatrix::solve (MatrixType &typ, const Matrix& b, octave_idx_type& info, double& rcon, solve_singularity_handler sing_handler, - bool singular_fallback) const + bool singular_fallback, blas_trans_type transt) const { ComplexMatrix tmp (b); - return solve (typ, tmp, info, rcon, sing_handler, singular_fallback); + return solve (typ, tmp, info, rcon, sing_handler, singular_fallback, transt); } ComplexMatrix @@ -2293,7 +2293,7 @@ ComplexMatrix::solve (MatrixType &mattype, const ComplexMatrix& b, octave_idx_type& info, double& rcon, solve_singularity_handler sing_handler, - bool singular_fallback) const + bool singular_fallback, blas_trans_type transt) const { ComplexMatrix retval; int typ = mattype.type (); @@ -2303,9 +2303,13 @@ // Only calculate the condition number for LU/Cholesky if (typ == MatrixType::Upper || typ == MatrixType::Permuted_Upper) - retval = utsolve (mattype, b, info, rcon, sing_handler, false); + retval = utsolve (mattype, b, info, rcon, sing_handler, false, transt); else if (typ == MatrixType::Lower || typ == MatrixType::Permuted_Lower) - retval = ltsolve (mattype, b, info, rcon, sing_handler, false); + retval = ltsolve (mattype, b, info, rcon, sing_handler, false, transt); + else if (transt == blas_trans) + return transpose ().solve (mattype, b, info, rcon, sing_handler, singular_fallback); + else if (transt == blas_conj_trans) + retval = hermitian ().solve (mattype, b, info, rcon, sing_handler, singular_fallback); else if (typ == MatrixType::Full || typ == MatrixType::Hermitian) retval = fsolve (mattype, b, info, rcon, sing_handler, true); else if (typ != MatrixType::Rectangular) @@ -2350,9 +2354,9 @@ ComplexColumnVector ComplexMatrix::solve (MatrixType &typ, const ColumnVector& b, octave_idx_type& info, double& rcon, - solve_singularity_handler sing_handler) const + solve_singularity_handler sing_handler, blas_trans_type transt) const { - return solve (typ, ComplexColumnVector (b), info, rcon, sing_handler); + return solve (typ, ComplexColumnVector (b), info, rcon, sing_handler, transt); } ComplexColumnVector @@ -2381,11 +2385,11 @@ ComplexColumnVector ComplexMatrix::solve (MatrixType &typ, const ComplexColumnVector& b, octave_idx_type& info, double& rcon, - solve_singularity_handler sing_handler) const + solve_singularity_handler sing_handler, blas_trans_type transt) const { ComplexMatrix tmp (b); - return solve (typ, tmp, info, rcon, sing_handler).column(static_cast<octave_idx_type> (0)); + return solve (typ, tmp, info, rcon, sing_handler, transt).column(static_cast<octave_idx_type> (0)); } ComplexMatrix @@ -2411,10 +2415,10 @@ ComplexMatrix ComplexMatrix::solve (const Matrix& b, octave_idx_type& info, double& rcon, - solve_singularity_handler sing_handler) const + solve_singularity_handler sing_handler, blas_trans_type transt) const { ComplexMatrix tmp (b); - return solve (tmp, info, rcon, sing_handler); + return solve (tmp, info, rcon, sing_handler, transt); } ComplexMatrix @@ -2440,10 +2444,10 @@ ComplexMatrix ComplexMatrix::solve (const ComplexMatrix& b, octave_idx_type& info, double& rcon, - solve_singularity_handler sing_handler) const + solve_singularity_handler sing_handler, blas_trans_type transt) const { MatrixType mattype (*this); - return solve (mattype, b, info, rcon, sing_handler); + return solve (mattype, b, info, rcon, sing_handler, true, transt); } ComplexColumnVector @@ -2471,9 +2475,9 @@ ComplexColumnVector ComplexMatrix::solve (const ColumnVector& b, octave_idx_type& info, double& rcon, - solve_singularity_handler sing_handler) const + solve_singularity_handler sing_handler, blas_trans_type transt) const { - return solve (ComplexColumnVector (b), info, rcon, sing_handler); + return solve (ComplexColumnVector (b), info, rcon, sing_handler, transt); } ComplexColumnVector @@ -2501,10 +2505,10 @@ ComplexColumnVector ComplexMatrix::solve (const ComplexColumnVector& b, octave_idx_type& info, double& rcon, - solve_singularity_handler sing_handler) const + solve_singularity_handler sing_handler, blas_trans_type transt) const { MatrixType mattype (*this); - return solve (mattype, b, info, rcon, sing_handler); + return solve (mattype, b, info, rcon, sing_handler, transt); } ComplexMatrix
--- a/liboctave/CMatrix.h Mon Sep 21 14:24:27 2009 +0200 +++ b/liboctave/CMatrix.h Wed Sep 23 10:00:16 2009 +0200 @@ -190,13 +190,14 @@ ComplexMatrix utsolve (MatrixType &typ, const ComplexMatrix& b, octave_idx_type& info, double& rcon, solve_singularity_handler sing_handler, - bool calc_cond = false) const; + bool calc_cond = false, + blas_trans_type transt = blas_no_trans) const; // Lower triangular matrix solvers ComplexMatrix ltsolve (MatrixType &typ, const ComplexMatrix& b, octave_idx_type& info, double& rcon, solve_singularity_handler sing_handler, - bool calc_cond = false) const; + bool calc_cond = false, blas_trans_type transt = blas_no_trans) const; // Full matrix solvers (umfpack/cholesky) ComplexMatrix fsolve (MatrixType &typ, const ComplexMatrix& b, @@ -213,7 +214,8 @@ octave_idx_type& info, double& rcon) const; ComplexMatrix solve (MatrixType &typ, const Matrix& b, octave_idx_type& info, double& rcon, solve_singularity_handler sing_handler, - bool singular_fallback = true) const; + bool singular_fallback = true, + blas_trans_type transt = blas_no_trans) const; ComplexMatrix solve (MatrixType &typ, const ComplexMatrix& b) const; ComplexMatrix solve (MatrixType &typ, const ComplexMatrix& b, @@ -223,7 +225,8 @@ ComplexMatrix solve (MatrixType &typ, const ComplexMatrix& b, octave_idx_type& info, double& rcon, solve_singularity_handler sing_handler, - bool singular_fallback = true) const; + bool singular_fallback = true, + blas_trans_type transt = blas_no_trans) const; ComplexColumnVector solve (MatrixType &typ, const ColumnVector& b) const; ComplexColumnVector solve (MatrixType &typ, const ColumnVector& b, @@ -232,7 +235,8 @@ octave_idx_type& info, double& rcon) const; ComplexColumnVector solve (MatrixType &typ, const ColumnVector& b, octave_idx_type& info, double& rcon, - solve_singularity_handler sing_handler) const; + solve_singularity_handler sing_handler, + blas_trans_type transt = blas_no_trans) const; ComplexColumnVector solve (MatrixType &typ, const ComplexColumnVector& b) const; @@ -242,35 +246,39 @@ octave_idx_type& info, double& rcon) const; ComplexColumnVector solve (MatrixType &typ, const ComplexColumnVector& b, octave_idx_type& info, double& rcon, - solve_singularity_handler sing_handler) const; + solve_singularity_handler sing_handler, + blas_trans_type transt = blas_no_trans) const; // Generic interface to solver with probing of type ComplexMatrix solve (const Matrix& b) const; ComplexMatrix solve (const Matrix& b, octave_idx_type& info) const; ComplexMatrix solve (const Matrix& b, octave_idx_type& info, double& rcon) const; ComplexMatrix solve (const Matrix& b, octave_idx_type& info, double& rcon, - solve_singularity_handler sing_handler) const; + solve_singularity_handler sing_handler, + blas_trans_type transt = blas_no_trans) const; ComplexMatrix solve (const ComplexMatrix& b) const; ComplexMatrix solve (const ComplexMatrix& b, octave_idx_type& info) const; ComplexMatrix solve (const ComplexMatrix& b, octave_idx_type& info, double& rcon) const; ComplexMatrix solve (const ComplexMatrix& b, octave_idx_type& info, double& rcon, - solve_singularity_handler sing_handler) const; + solve_singularity_handler sing_handler, + blas_trans_type transt = blas_no_trans) const; ComplexColumnVector solve (const ColumnVector& b) const; ComplexColumnVector solve (const ColumnVector& b, octave_idx_type& info) const; ComplexColumnVector solve (const ColumnVector& b, octave_idx_type& info, double& rcon) const; ComplexColumnVector solve (const ColumnVector& b, octave_idx_type& info, double& rcon, - solve_singularity_handler sing_handler) const; + solve_singularity_handler sing_handler, + blas_trans_type transt = blas_no_trans) const; ComplexColumnVector solve (const ComplexColumnVector& b) const; ComplexColumnVector solve (const ComplexColumnVector& b, octave_idx_type& info) const; ComplexColumnVector solve (const ComplexColumnVector& b, octave_idx_type& info, double& rcon) const; ComplexColumnVector solve (const ComplexColumnVector& b, octave_idx_type& info, - double& rcon, - solve_singularity_handler sing_handler) const; + double& rcon, solve_singularity_handler sing_handler, + blas_trans_type transt = blas_no_trans) const; ComplexMatrix lssolve (const Matrix& b) const; ComplexMatrix lssolve (const Matrix& b, octave_idx_type& info) const;
--- a/liboctave/ChangeLog Mon Sep 21 14:24:27 2009 +0200 +++ b/liboctave/ChangeLog Wed Sep 23 10:00:16 2009 +0200 @@ -1,3 +1,17 @@ +2009-09-23 Jaroslav Hajek <highegg@gmail.com> + + * mx-defs.h (blas_trans_type): New enum. + (get_blas_char): New inline func. + * dMatrix.cc (Matrix::utsolve, Matrix::ltsolve, Matrix::solve): + Support transt parameter. + * fMatrix.cc (FloatMatrix::utsolve, FloatMatrix::ltsolve, + FloatMatrix::solve): Ditto. + * CMatrix.cc (ComplexMatrix::utsolve, ComplexMatrix::ltsolve, + ComplexMatrix::solve): Ditto. + * fCMatrix.cc (FloatComplexMatrix::utsolve, + FloatComplexMatrix::ltsolve, FloatComplexMatrix::solve): Ditto. + * dMatrix.h, fMatrix.h, CMatrix.h, fCMatrix.h: Update. + 2009-09-21 Jaroslav Hajek <highegg@gmail.com> * mx-op-defs.h (VS_BIN_OP, SV_BIN_OP, VV_BIN_OP): Simplify.
--- a/liboctave/dMatrix.cc Mon Sep 21 14:24:27 2009 +0200 +++ b/liboctave/dMatrix.cc Wed Sep 23 10:00:16 2009 +0200 @@ -1523,7 +1523,7 @@ Matrix Matrix::utsolve (MatrixType &mattype, const Matrix& b, octave_idx_type& info, double& rcon, solve_singularity_handler sing_handler, - bool calc_cond) const + bool calc_cond, blas_trans_type transt) const { Matrix retval; @@ -1599,7 +1599,7 @@ double *result = retval.fortran_vec (); char uplo = 'U'; - char trans = 'N'; + char trans = get_blas_char (transt); char dia = 'N'; F77_XFCN (dtrtrs, DTRTRS, (F77_CONST_CHAR_ARG2 (&uplo, 1), @@ -1623,7 +1623,7 @@ Matrix Matrix::ltsolve (MatrixType &mattype, const Matrix& b, octave_idx_type& info, double& rcon, solve_singularity_handler sing_handler, - bool calc_cond) const + bool calc_cond, blas_trans_type transt) const { Matrix retval; @@ -1699,7 +1699,7 @@ double *result = retval.fortran_vec (); char uplo = 'L'; - char trans = 'N'; + char trans = get_blas_char (transt); char dia = 'N'; F77_XFCN (dtrtrs, DTRTRS, (F77_CONST_CHAR_ARG2 (&uplo, 1), @@ -1920,7 +1920,7 @@ Matrix Matrix::solve (MatrixType &mattype, const Matrix& b, octave_idx_type& info, double& rcon, solve_singularity_handler sing_handler, - bool singular_fallback) const + bool singular_fallback, blas_trans_type transt) const { Matrix retval; int typ = mattype.type (); @@ -1930,9 +1930,11 @@ // Only calculate the condition number for LU/Cholesky if (typ == MatrixType::Upper || typ == MatrixType::Permuted_Upper) - retval = utsolve (mattype, b, info, rcon, sing_handler, false); + retval = utsolve (mattype, b, info, rcon, sing_handler, false, transt); else if (typ == MatrixType::Lower || typ == MatrixType::Permuted_Lower) - retval = ltsolve (mattype, b, info, rcon, sing_handler, false); + retval = ltsolve (mattype, b, info, rcon, sing_handler, false, transt); + else if (transt == blas_trans || transt == blas_conj_trans) + return transpose ().solve (mattype, b, info, rcon, sing_handler, singular_fallback); else if (typ == MatrixType::Full || typ == MatrixType::Hermitian) retval = fsolve (mattype, b, info, rcon, sing_handler, true); else if (typ != MatrixType::Rectangular) @@ -1977,10 +1979,10 @@ ComplexMatrix Matrix::solve (MatrixType &typ, const ComplexMatrix& b, octave_idx_type& info, double& rcon, solve_singularity_handler sing_handler, - bool singular_fallback) const + bool singular_fallback, blas_trans_type transt) const { ComplexMatrix tmp (*this); - return tmp.solve (typ, b, info, rcon, sing_handler, singular_fallback); + return tmp.solve (typ, b, info, rcon, sing_handler, singular_fallback, transt); } ColumnVector @@ -2007,10 +2009,10 @@ ColumnVector Matrix::solve (MatrixType &typ, const ColumnVector& b, octave_idx_type& info, - double& rcon, solve_singularity_handler sing_handler) const + double& rcon, solve_singularity_handler sing_handler, blas_trans_type transt) const { Matrix tmp (b); - return solve (typ, tmp, info, rcon, sing_handler).column(static_cast<octave_idx_type> (0)); + return solve (typ, tmp, info, rcon, sing_handler, transt).column(static_cast<octave_idx_type> (0)); } ComplexColumnVector @@ -2039,10 +2041,10 @@ ComplexColumnVector Matrix::solve (MatrixType &typ, const ComplexColumnVector& b, octave_idx_type& info, double& rcon, - solve_singularity_handler sing_handler) const + solve_singularity_handler sing_handler, blas_trans_type transt) const { ComplexMatrix tmp (*this); - return tmp.solve(typ, b, info, rcon, sing_handler); + return tmp.solve(typ, b, info, rcon, sing_handler, transt); } Matrix @@ -2068,10 +2070,10 @@ Matrix Matrix::solve (const Matrix& b, octave_idx_type& info, - double& rcon, solve_singularity_handler sing_handler) const + double& rcon, solve_singularity_handler sing_handler, blas_trans_type transt) const { MatrixType mattype (*this); - return solve (mattype, b, info, rcon, sing_handler); + return solve (mattype, b, info, rcon, sing_handler, true, transt); } ComplexMatrix @@ -2097,10 +2099,10 @@ ComplexMatrix Matrix::solve (const ComplexMatrix& b, octave_idx_type& info, double& rcon, - solve_singularity_handler sing_handler) const + solve_singularity_handler sing_handler, blas_trans_type transt) const { ComplexMatrix tmp (*this); - return tmp.solve (b, info, rcon, sing_handler); + return tmp.solve (b, info, rcon, sing_handler, transt); } ColumnVector @@ -2125,10 +2127,10 @@ ColumnVector Matrix::solve (const ColumnVector& b, octave_idx_type& info, double& rcon, - solve_singularity_handler sing_handler) const + solve_singularity_handler sing_handler, blas_trans_type transt) const { MatrixType mattype (*this); - return solve (mattype, b, info, rcon, sing_handler); + return solve (mattype, b, info, rcon, sing_handler, transt); } ComplexColumnVector @@ -2154,10 +2156,10 @@ ComplexColumnVector Matrix::solve (const ComplexColumnVector& b, octave_idx_type& info, double& rcon, - solve_singularity_handler sing_handler) const + solve_singularity_handler sing_handler, blas_trans_type transt) const { ComplexMatrix tmp (*this); - return tmp.solve (b, info, rcon, sing_handler); + return tmp.solve (b, info, rcon, sing_handler, transt); } Matrix
--- a/liboctave/dMatrix.h Mon Sep 21 14:24:27 2009 +0200 +++ b/liboctave/dMatrix.h Wed Sep 23 10:00:16 2009 +0200 @@ -162,12 +162,12 @@ // Upper triangular matrix solvers Matrix utsolve (MatrixType &typ, const Matrix& b, octave_idx_type& info, double& rcon, solve_singularity_handler sing_handler, - bool calc_cond = false) const; + bool calc_cond = false, blas_trans_type transt = blas_no_trans) const; // Lower triangular matrix solvers Matrix ltsolve (MatrixType &typ, const Matrix& b, octave_idx_type& info, double& rcon, solve_singularity_handler sing_handler, - bool calc_cond = false) const; + bool calc_cond = false, blas_trans_type transt = blas_no_trans) const; // Full matrix solvers (lu/cholesky) Matrix fsolve (MatrixType &typ, const Matrix& b, octave_idx_type& info, @@ -182,7 +182,7 @@ double& rcon) const; Matrix solve (MatrixType &typ, const Matrix& b, octave_idx_type& info, double& rcon, solve_singularity_handler sing_handler, - bool singular_fallback = true) const; + bool singular_fallback = true, blas_trans_type transt = blas_no_trans) const; ComplexMatrix solve (MatrixType &typ, const ComplexMatrix& b) const; ComplexMatrix solve (MatrixType &typ, const ComplexMatrix& b, @@ -192,7 +192,8 @@ ComplexMatrix solve (MatrixType &typ, const ComplexMatrix& b, octave_idx_type& info, double& rcon, solve_singularity_handler sing_handler, - bool singular_fallback = true) const; + bool singular_fallback = true, + blas_trans_type transt = blas_no_trans) const; ColumnVector solve (MatrixType &typ, const ColumnVector& b) const; ColumnVector solve (MatrixType &typ, const ColumnVector& b, @@ -201,7 +202,8 @@ octave_idx_type& info, double& rcon) const; ColumnVector solve (MatrixType &typ, const ColumnVector& b, octave_idx_type& info, double& rcon, - solve_singularity_handler sing_handler) const; + solve_singularity_handler sing_handler, + blas_trans_type transt = blas_no_trans) const; ComplexColumnVector solve (MatrixType &typ, const ComplexColumnVector& b) const; @@ -211,34 +213,38 @@ octave_idx_type& info, double& rcon) const; ComplexColumnVector solve (MatrixType &typ, const ComplexColumnVector& b, octave_idx_type& info, double& rcon, - solve_singularity_handler sing_handler) const; + solve_singularity_handler sing_handler, + blas_trans_type transt = blas_no_trans) const; // Generic interface to solver with probing of type Matrix solve (const Matrix& b) const; Matrix solve (const Matrix& b, octave_idx_type& info) const; Matrix solve (const Matrix& b, octave_idx_type& info, double& rcon) const; Matrix solve (const Matrix& b, octave_idx_type& info, double& rcon, - solve_singularity_handler sing_handler) const; + solve_singularity_handler sing_handler, + blas_trans_type transt = blas_no_trans) const; ComplexMatrix solve (const ComplexMatrix& b) const; ComplexMatrix solve (const ComplexMatrix& b, octave_idx_type& info) const; ComplexMatrix solve (const ComplexMatrix& b, octave_idx_type& info, double& rcon) const; ComplexMatrix solve (const ComplexMatrix& b, octave_idx_type& info, double& rcon, - solve_singularity_handler sing_handler) const; + solve_singularity_handler sing_handler, + blas_trans_type transt = blas_no_trans) const; ColumnVector solve (const ColumnVector& b) const; ColumnVector solve (const ColumnVector& b, octave_idx_type& info) const; ColumnVector solve (const ColumnVector& b, octave_idx_type& info, double& rcon) const; ColumnVector solve (const ColumnVector& b, octave_idx_type& info, double& rcon, - solve_singularity_handler sing_handler) const; + solve_singularity_handler sing_handler, + blas_trans_type transt = blas_no_trans) const; ComplexColumnVector solve (const ComplexColumnVector& b) const; ComplexColumnVector solve (const ComplexColumnVector& b, octave_idx_type& info) const; ComplexColumnVector solve (const ComplexColumnVector& b, octave_idx_type& info, double& rcon) const; ComplexColumnVector solve (const ComplexColumnVector& b, octave_idx_type& info, - double& rcon, - solve_singularity_handler sing_handler) const; + double& rcon, solve_singularity_handler sing_handler, + blas_trans_type transt = blas_no_trans) const; // Singular solvers Matrix lssolve (const Matrix& b) const;
--- a/liboctave/fCMatrix.cc Mon Sep 21 14:24:27 2009 +0200 +++ b/liboctave/fCMatrix.cc Wed Sep 23 10:00:16 2009 +0200 @@ -1845,7 +1845,7 @@ FloatComplexMatrix::utsolve (MatrixType &mattype, const FloatComplexMatrix& b, octave_idx_type& info, float& rcon, solve_singularity_handler sing_handler, - bool calc_cond) const + bool calc_cond, blas_trans_type transt) const { FloatComplexMatrix retval; @@ -1921,7 +1921,7 @@ FloatComplex *result = retval.fortran_vec (); char uplo = 'U'; - char trans = 'N'; + char trans = get_blas_char (transt); char dia = 'N'; F77_XFCN (ctrtrs, CTRTRS, (F77_CONST_CHAR_ARG2 (&uplo, 1), @@ -1946,7 +1946,7 @@ FloatComplexMatrix::ltsolve (MatrixType &mattype, const FloatComplexMatrix& b, octave_idx_type& info, float& rcon, solve_singularity_handler sing_handler, - bool calc_cond) const + bool calc_cond, blas_trans_type transt) const { FloatComplexMatrix retval; @@ -2022,7 +2022,7 @@ FloatComplex *result = retval.fortran_vec (); char uplo = 'L'; - char trans = 'N'; + char trans = get_blas_char (transt); char dia = 'N'; F77_XFCN (ctrtrs, CTRTRS, (F77_CONST_CHAR_ARG2 (&uplo, 1), @@ -2253,10 +2253,10 @@ FloatComplexMatrix FloatComplexMatrix::solve (MatrixType &typ, const FloatMatrix& b, octave_idx_type& info, float& rcon, solve_singularity_handler sing_handler, - bool singular_fallback) const + bool singular_fallback, blas_trans_type transt) const { FloatComplexMatrix tmp (b); - return solve (typ, tmp, info, rcon, sing_handler, singular_fallback); + return solve (typ, tmp, info, rcon, sing_handler, singular_fallback, transt); } FloatComplexMatrix @@ -2286,7 +2286,7 @@ FloatComplexMatrix::solve (MatrixType &mattype, const FloatComplexMatrix& b, octave_idx_type& info, float& rcon, solve_singularity_handler sing_handler, - bool singular_fallback) const + bool singular_fallback, blas_trans_type transt) const { FloatComplexMatrix retval; int typ = mattype.type (); @@ -2296,9 +2296,13 @@ // Only calculate the condition number for LU/Cholesky if (typ == MatrixType::Upper || typ == MatrixType::Permuted_Upper) - retval = utsolve (mattype, b, info, rcon, sing_handler, false); + retval = utsolve (mattype, b, info, rcon, sing_handler, false, transt); else if (typ == MatrixType::Lower || typ == MatrixType::Permuted_Lower) - retval = ltsolve (mattype, b, info, rcon, sing_handler, false); + retval = ltsolve (mattype, b, info, rcon, sing_handler, false, transt); + else if (transt == blas_trans) + return transpose ().solve (mattype, b, info, rcon, sing_handler, singular_fallback); + else if (transt == blas_conj_trans) + retval = hermitian ().solve (mattype, b, info, rcon, sing_handler, singular_fallback); else if (typ == MatrixType::Full || typ == MatrixType::Hermitian) retval = fsolve (mattype, b, info, rcon, sing_handler, true); else if (typ != MatrixType::Rectangular) @@ -2343,9 +2347,9 @@ FloatComplexColumnVector FloatComplexMatrix::solve (MatrixType &typ, const FloatColumnVector& b, octave_idx_type& info, float& rcon, - solve_singularity_handler sing_handler) const + solve_singularity_handler sing_handler, blas_trans_type transt) const { - return solve (typ, FloatComplexColumnVector (b), info, rcon, sing_handler); + return solve (typ, FloatComplexColumnVector (b), info, rcon, sing_handler, transt); } FloatComplexColumnVector @@ -2374,11 +2378,11 @@ FloatComplexColumnVector FloatComplexMatrix::solve (MatrixType &typ, const FloatComplexColumnVector& b, octave_idx_type& info, float& rcon, - solve_singularity_handler sing_handler) const + solve_singularity_handler sing_handler, blas_trans_type transt) const { FloatComplexMatrix tmp (b); - return solve (typ, tmp, info, rcon, sing_handler).column(static_cast<octave_idx_type> (0)); + return solve (typ, tmp, info, rcon, sing_handler, transt).column(static_cast<octave_idx_type> (0)); } FloatComplexMatrix @@ -2404,10 +2408,10 @@ FloatComplexMatrix FloatComplexMatrix::solve (const FloatMatrix& b, octave_idx_type& info, float& rcon, - solve_singularity_handler sing_handler) const + solve_singularity_handler sing_handler, blas_trans_type transt) const { FloatComplexMatrix tmp (b); - return solve (tmp, info, rcon, sing_handler); + return solve (tmp, info, rcon, sing_handler, transt); } FloatComplexMatrix @@ -2433,10 +2437,10 @@ FloatComplexMatrix FloatComplexMatrix::solve (const FloatComplexMatrix& b, octave_idx_type& info, float& rcon, - solve_singularity_handler sing_handler) const + solve_singularity_handler sing_handler, blas_trans_type transt) const { MatrixType mattype (*this); - return solve (mattype, b, info, rcon, sing_handler); + return solve (mattype, b, info, rcon, sing_handler, true, transt); } FloatComplexColumnVector @@ -2464,9 +2468,9 @@ FloatComplexColumnVector FloatComplexMatrix::solve (const FloatColumnVector& b, octave_idx_type& info, float& rcon, - solve_singularity_handler sing_handler) const + solve_singularity_handler sing_handler, blas_trans_type transt) const { - return solve (FloatComplexColumnVector (b), info, rcon, sing_handler); + return solve (FloatComplexColumnVector (b), info, rcon, sing_handler, transt); } FloatComplexColumnVector @@ -2494,10 +2498,10 @@ FloatComplexColumnVector FloatComplexMatrix::solve (const FloatComplexColumnVector& b, octave_idx_type& info, float& rcon, - solve_singularity_handler sing_handler) const + solve_singularity_handler sing_handler, blas_trans_type transt) const { MatrixType mattype (*this); - return solve (mattype, b, info, rcon, sing_handler); + return solve (mattype, b, info, rcon, sing_handler, transt); } FloatComplexMatrix
--- a/liboctave/fCMatrix.h Mon Sep 21 14:24:27 2009 +0200 +++ b/liboctave/fCMatrix.h Wed Sep 23 10:00:16 2009 +0200 @@ -190,13 +190,13 @@ FloatComplexMatrix utsolve (MatrixType &typ, const FloatComplexMatrix& b, octave_idx_type& info, float& rcon, solve_singularity_handler sing_handler, - bool calc_cond = false) const; + bool calc_cond = false, blas_trans_type transt = blas_no_trans) const; // Lower triangular matrix solvers FloatComplexMatrix ltsolve (MatrixType &typ, const FloatComplexMatrix& b, octave_idx_type& info, float& rcon, solve_singularity_handler sing_handler, - bool calc_cond = false) const; + bool calc_cond = false, blas_trans_type transt = blas_no_trans) const; // Full matrix solvers (umfpack/cholesky) FloatComplexMatrix fsolve (MatrixType &typ, const FloatComplexMatrix& b, @@ -213,7 +213,8 @@ octave_idx_type& info, float& rcon) const; FloatComplexMatrix solve (MatrixType &typ, const FloatMatrix& b, octave_idx_type& info, float& rcon, solve_singularity_handler sing_handler, - bool singular_fallback = true) const; + bool singular_fallback = true, + blas_trans_type transt = blas_no_trans) const; FloatComplexMatrix solve (MatrixType &typ, const FloatComplexMatrix& b) const; FloatComplexMatrix solve (MatrixType &typ, const FloatComplexMatrix& b, @@ -223,7 +224,8 @@ FloatComplexMatrix solve (MatrixType &typ, const FloatComplexMatrix& b, octave_idx_type& info, float& rcon, solve_singularity_handler sing_handler, - bool singular_fallback = true) const; + bool singular_fallback = true, + blas_trans_type transt = blas_no_trans) const; FloatComplexColumnVector solve (MatrixType &typ, const FloatColumnVector& b) const; FloatComplexColumnVector solve (MatrixType &typ, const FloatColumnVector& b, @@ -232,7 +234,8 @@ octave_idx_type& info, float& rcon) const; FloatComplexColumnVector solve (MatrixType &typ, const FloatColumnVector& b, octave_idx_type& info, float& rcon, - solve_singularity_handler sing_handler) const; + solve_singularity_handler sing_handler, + blas_trans_type transt = blas_no_trans) const; FloatComplexColumnVector solve (MatrixType &typ, const FloatComplexColumnVector& b) const; @@ -242,27 +245,31 @@ octave_idx_type& info, float& rcon) const; FloatComplexColumnVector solve (MatrixType &typ, const FloatComplexColumnVector& b, octave_idx_type& info, float& rcon, - solve_singularity_handler sing_handler) const; + solve_singularity_handler sing_handler, + blas_trans_type transt = blas_no_trans) const; // Generic interface to solver with probing of type FloatComplexMatrix solve (const FloatMatrix& b) const; FloatComplexMatrix solve (const FloatMatrix& b, octave_idx_type& info) const; FloatComplexMatrix solve (const FloatMatrix& b, octave_idx_type& info, float& rcon) const; FloatComplexMatrix solve (const FloatMatrix& b, octave_idx_type& info, float& rcon, - solve_singularity_handler sing_handler) const; + solve_singularity_handler sing_handler, + blas_trans_type transt = blas_no_trans) const; FloatComplexMatrix solve (const FloatComplexMatrix& b) const; FloatComplexMatrix solve (const FloatComplexMatrix& b, octave_idx_type& info) const; FloatComplexMatrix solve (const FloatComplexMatrix& b, octave_idx_type& info, float& rcon) const; FloatComplexMatrix solve (const FloatComplexMatrix& b, octave_idx_type& info, float& rcon, - solve_singularity_handler sing_handler) const; + solve_singularity_handler sing_handler, + blas_trans_type transt = blas_no_trans) const; FloatComplexColumnVector solve (const FloatColumnVector& b) const; FloatComplexColumnVector solve (const FloatColumnVector& b, octave_idx_type& info) const; FloatComplexColumnVector solve (const FloatColumnVector& b, octave_idx_type& info, float& rcon) const; FloatComplexColumnVector solve (const FloatColumnVector& b, octave_idx_type& info, float& rcon, - solve_singularity_handler sing_handler) const; + solve_singularity_handler sing_handler, + blas_trans_type transt = blas_no_trans) const; FloatComplexColumnVector solve (const FloatComplexColumnVector& b) const; FloatComplexColumnVector solve (const FloatComplexColumnVector& b, octave_idx_type& info) const; @@ -270,7 +277,8 @@ float& rcon) const; FloatComplexColumnVector solve (const FloatComplexColumnVector& b, octave_idx_type& info, float& rcon, - solve_singularity_handler sing_handler) const; + solve_singularity_handler sing_handler, + blas_trans_type transt = blas_no_trans) const; FloatComplexMatrix lssolve (const FloatMatrix& b) const; FloatComplexMatrix lssolve (const FloatMatrix& b, octave_idx_type& info) const;
--- a/liboctave/fMatrix.cc Mon Sep 21 14:24:27 2009 +0200 +++ b/liboctave/fMatrix.cc Wed Sep 23 10:00:16 2009 +0200 @@ -1522,7 +1522,7 @@ FloatMatrix FloatMatrix::utsolve (MatrixType &mattype, const FloatMatrix& b, octave_idx_type& info, float& rcon, solve_singularity_handler sing_handler, - bool calc_cond) const + bool calc_cond, blas_trans_type transt) const { FloatMatrix retval; @@ -1598,7 +1598,7 @@ float *result = retval.fortran_vec (); char uplo = 'U'; - char trans = 'N'; + char trans = get_blas_char (transt); char dia = 'N'; F77_XFCN (strtrs, STRTRS, (F77_CONST_CHAR_ARG2 (&uplo, 1), @@ -1622,7 +1622,7 @@ FloatMatrix FloatMatrix::ltsolve (MatrixType &mattype, const FloatMatrix& b, octave_idx_type& info, float& rcon, solve_singularity_handler sing_handler, - bool calc_cond) const + bool calc_cond, blas_trans_type transt) const { FloatMatrix retval; @@ -1698,7 +1698,7 @@ float *result = retval.fortran_vec (); char uplo = 'L'; - char trans = 'N'; + char trans = get_blas_char (transt); char dia = 'N'; F77_XFCN (strtrs, STRTRS, (F77_CONST_CHAR_ARG2 (&uplo, 1), @@ -1919,7 +1919,7 @@ FloatMatrix FloatMatrix::solve (MatrixType &mattype, const FloatMatrix& b, octave_idx_type& info, float& rcon, solve_singularity_handler sing_handler, - bool singular_fallback) const + bool singular_fallback, blas_trans_type transt) const { FloatMatrix retval; int typ = mattype.type (); @@ -1929,9 +1929,11 @@ // Only calculate the condition number for LU/Cholesky if (typ == MatrixType::Upper || typ == MatrixType::Permuted_Upper) - retval = utsolve (mattype, b, info, rcon, sing_handler, false); + retval = utsolve (mattype, b, info, rcon, sing_handler, false, transt); else if (typ == MatrixType::Lower || typ == MatrixType::Permuted_Lower) - retval = ltsolve (mattype, b, info, rcon, sing_handler, false); + retval = ltsolve (mattype, b, info, rcon, sing_handler, false, transt); + else if (transt == blas_trans || transt == blas_conj_trans) + return transpose ().solve (mattype, b, info, rcon, sing_handler, singular_fallback); else if (typ == MatrixType::Full || typ == MatrixType::Hermitian) retval = fsolve (mattype, b, info, rcon, sing_handler, true); else if (typ != MatrixType::Rectangular) @@ -1976,10 +1978,10 @@ FloatComplexMatrix FloatMatrix::solve (MatrixType &typ, const FloatComplexMatrix& b, octave_idx_type& info, float& rcon, solve_singularity_handler sing_handler, - bool singular_fallback) const + bool singular_fallback, blas_trans_type transt) const { FloatComplexMatrix tmp (*this); - return tmp.solve (typ, b, info, rcon, sing_handler, singular_fallback); + return tmp.solve (typ, b, info, rcon, sing_handler, singular_fallback, transt); } FloatColumnVector @@ -2006,10 +2008,10 @@ FloatColumnVector FloatMatrix::solve (MatrixType &typ, const FloatColumnVector& b, octave_idx_type& info, - float& rcon, solve_singularity_handler sing_handler) const + float& rcon, solve_singularity_handler sing_handler, blas_trans_type transt) const { FloatMatrix tmp (b); - return solve (typ, tmp, info, rcon, sing_handler).column(static_cast<octave_idx_type> (0)); + return solve (typ, tmp, info, rcon, sing_handler, transt).column(static_cast<octave_idx_type> (0)); } FloatComplexColumnVector @@ -2038,10 +2040,10 @@ FloatComplexColumnVector FloatMatrix::solve (MatrixType &typ, const FloatComplexColumnVector& b, octave_idx_type& info, float& rcon, - solve_singularity_handler sing_handler) const + solve_singularity_handler sing_handler, blas_trans_type transt) const { FloatComplexMatrix tmp (*this); - return tmp.solve(typ, b, info, rcon, sing_handler); + return tmp.solve(typ, b, info, rcon, sing_handler, transt); } FloatMatrix @@ -2067,10 +2069,10 @@ FloatMatrix FloatMatrix::solve (const FloatMatrix& b, octave_idx_type& info, - float& rcon, solve_singularity_handler sing_handler) const + float& rcon, solve_singularity_handler sing_handler, blas_trans_type transt) const { MatrixType mattype (*this); - return solve (mattype, b, info, rcon, sing_handler); + return solve (mattype, b, info, rcon, sing_handler, true, transt); } FloatComplexMatrix @@ -2096,10 +2098,10 @@ FloatComplexMatrix FloatMatrix::solve (const FloatComplexMatrix& b, octave_idx_type& info, float& rcon, - solve_singularity_handler sing_handler) const + solve_singularity_handler sing_handler, blas_trans_type transt) const { FloatComplexMatrix tmp (*this); - return tmp.solve (b, info, rcon, sing_handler); + return tmp.solve (b, info, rcon, sing_handler, transt); } FloatColumnVector @@ -2124,10 +2126,10 @@ FloatColumnVector FloatMatrix::solve (const FloatColumnVector& b, octave_idx_type& info, float& rcon, - solve_singularity_handler sing_handler) const + solve_singularity_handler sing_handler, blas_trans_type transt) const { MatrixType mattype (*this); - return solve (mattype, b, info, rcon, sing_handler); + return solve (mattype, b, info, rcon, sing_handler, transt); } FloatComplexColumnVector @@ -2153,10 +2155,10 @@ FloatComplexColumnVector FloatMatrix::solve (const FloatComplexColumnVector& b, octave_idx_type& info, float& rcon, - solve_singularity_handler sing_handler) const + solve_singularity_handler sing_handler, blas_trans_type transt) const { FloatComplexMatrix tmp (*this); - return tmp.solve (b, info, rcon, sing_handler); + return tmp.solve (b, info, rcon, sing_handler, transt); } FloatMatrix
--- a/liboctave/fMatrix.h Mon Sep 21 14:24:27 2009 +0200 +++ b/liboctave/fMatrix.h Wed Sep 23 10:00:16 2009 +0200 @@ -163,12 +163,12 @@ // Upper triangular matrix solvers FloatMatrix utsolve (MatrixType &typ, const FloatMatrix& b, octave_idx_type& info, float& rcon, solve_singularity_handler sing_handler, - bool calc_cond = false) const; + bool calc_cond = false, blas_trans_type transt = blas_no_trans) const; // Lower triangular matrix solvers FloatMatrix ltsolve (MatrixType &typ, const FloatMatrix& b, octave_idx_type& info, float& rcon, solve_singularity_handler sing_handler, - bool calc_cond = false) const; + bool calc_cond = false, blas_trans_type transt = blas_no_trans) const; // Full matrix solvers (lu/cholesky) FloatMatrix fsolve (MatrixType &typ, const FloatMatrix& b, octave_idx_type& info, @@ -183,7 +183,7 @@ float& rcon) const; FloatMatrix solve (MatrixType &typ, const FloatMatrix& b, octave_idx_type& info, float& rcon, solve_singularity_handler sing_handler, - bool singular_fallback = true) const; + bool singular_fallback = true, blas_trans_type transt = blas_no_trans) const; FloatComplexMatrix solve (MatrixType &typ, const FloatComplexMatrix& b) const; FloatComplexMatrix solve (MatrixType &typ, const FloatComplexMatrix& b, @@ -193,7 +193,8 @@ FloatComplexMatrix solve (MatrixType &typ, const FloatComplexMatrix& b, octave_idx_type& info, float& rcon, solve_singularity_handler sing_handler, - bool singular_fallback = true) const; + bool singular_fallback = true, + blas_trans_type transt = blas_no_trans) const; FloatColumnVector solve (MatrixType &typ, const FloatColumnVector& b) const; FloatColumnVector solve (MatrixType &typ, const FloatColumnVector& b, @@ -202,7 +203,8 @@ octave_idx_type& info, float& rcon) const; FloatColumnVector solve (MatrixType &typ, const FloatColumnVector& b, octave_idx_type& info, float& rcon, - solve_singularity_handler sing_handler) const; + solve_singularity_handler sing_handler, + blas_trans_type transt = blas_no_trans) const; FloatComplexColumnVector solve (MatrixType &typ, const FloatComplexColumnVector& b) const; @@ -212,34 +214,37 @@ octave_idx_type& info, float& rcon) const; FloatComplexColumnVector solve (MatrixType &typ, const FloatComplexColumnVector& b, octave_idx_type& info, float& rcon, - solve_singularity_handler sing_handler) const; + solve_singularity_handler sing_handler, blas_trans_type transt = blas_no_trans) const; // Generic interface to solver with probing of type FloatMatrix solve (const FloatMatrix& b) const; FloatMatrix solve (const FloatMatrix& b, octave_idx_type& info) const; FloatMatrix solve (const FloatMatrix& b, octave_idx_type& info, float& rcon) const; FloatMatrix solve (const FloatMatrix& b, octave_idx_type& info, float& rcon, - solve_singularity_handler sing_handler) const; + solve_singularity_handler sing_handler, + blas_trans_type transt = blas_no_trans) const; FloatComplexMatrix solve (const FloatComplexMatrix& b) const; FloatComplexMatrix solve (const FloatComplexMatrix& b, octave_idx_type& info) const; FloatComplexMatrix solve (const FloatComplexMatrix& b, octave_idx_type& info, float& rcon) const; FloatComplexMatrix solve (const FloatComplexMatrix& b, octave_idx_type& info, float& rcon, - solve_singularity_handler sing_handler) const; + solve_singularity_handler sing_handler, + blas_trans_type transt = blas_no_trans) const; FloatColumnVector solve (const FloatColumnVector& b) const; FloatColumnVector solve (const FloatColumnVector& b, octave_idx_type& info) const; FloatColumnVector solve (const FloatColumnVector& b, octave_idx_type& info, float& rcon) const; FloatColumnVector solve (const FloatColumnVector& b, octave_idx_type& info, float& rcon, - solve_singularity_handler sing_handler) const; + solve_singularity_handler sing_handler, + blas_trans_type transt = blas_no_trans) const; FloatComplexColumnVector solve (const FloatComplexColumnVector& b) const; FloatComplexColumnVector solve (const FloatComplexColumnVector& b, octave_idx_type& info) const; FloatComplexColumnVector solve (const FloatComplexColumnVector& b, octave_idx_type& info, - float& rcon) const; + float& rcon) const; FloatComplexColumnVector solve (const FloatComplexColumnVector& b, octave_idx_type& info, - float& rcon, - solve_singularity_handler sing_handler) const; + float& rcon, solve_singularity_handler sing_handler, + blas_trans_type transt = blas_no_trans) const; // Singular solvers FloatMatrix lssolve (const FloatMatrix& b) const;
--- a/liboctave/mx-defs.h Mon Sep 21 14:24:27 2009 +0200 +++ b/liboctave/mx-defs.h Wed Sep 23 10:00:16 2009 +0200 @@ -126,6 +126,20 @@ typedef float (*f_fc_Mapper)(const FloatComplex&); typedef FloatComplex (*fc_fc_Mapper)(const FloatComplex&); +enum blas_trans_type +{ + blas_no_trans = 'N', + blas_trans = 'T', + blas_conj_trans = 'C' +}; + +inline char +get_blas_char (blas_trans_type transt) +{ + return static_cast<char> (transt); +} + + #endif #endif
--- a/src/ChangeLog Mon Sep 21 14:24:27 2009 +0200 +++ b/src/ChangeLog Wed Sep 23 10:00:16 2009 +0200 @@ -1,3 +1,18 @@ +2009-09-23 Jaroslav Hajek <highegg@gmail.com> + + * ov.h (octave_value::op_trans_ldiv, op_herm_ldiv): New enum constants. + * ov.cc (decompose_binary_op, binary_op_fcn_name): Support them. + * xdiv.h: Include mx-defs.h, delete forward decls. + * xdiv.cc (xleftdiv): Support transt parameter. + (xdiv): Optimize. + * pt-cbinop.cc (simplify_ldiv_op): New static func. + (maybe_compound_binary_expression): Try it. + * OPERATORS/op-m-m.cc: Define and install trans_ldiv handler. + * OPERATORS/op-fm-fm.cc: Ditto. + * OPERATORS/op-cm-cm.cc: Define and install trans_ldiv and herm_ldiv + handlers. + * OPERATORS/op-fcm-fcm.cc: Ditto. + 2009-09-19 Jaroslav Hajek <highegg@gmail.com> * ov.h (octave_value_extract): New template function.
--- a/src/OPERATORS/op-cm-cm.cc Mon Sep 21 14:24:27 2009 +0200 +++ b/src/OPERATORS/op-cm-cm.cc Wed Sep 23 10:00:16 2009 +0200 @@ -137,6 +137,30 @@ true, true, v2.complex_matrix_value ())); } +DEFBINOP (trans_ldiv, complex_matrix, complex_matrix) +{ + CAST_BINOP_ARGS (const octave_complex_matrix&, const octave_complex_matrix&); + MatrixType typ = v1.matrix_type (); + + ComplexMatrix ret = xleftdiv (v1.complex_matrix_value (), + v2.complex_matrix_value (), typ, blas_trans); + + v1.matrix_type (typ); + return ret; +} + +DEFBINOP (herm_ldiv, complex_matrix, complex_matrix) +{ + CAST_BINOP_ARGS (const octave_complex_matrix&, const octave_complex_matrix&); + MatrixType typ = v1.matrix_type (); + + ComplexMatrix ret = xleftdiv (v1.complex_matrix_value (), + v2.complex_matrix_value (), typ, blas_conj_trans); + + v1.matrix_type (typ); + return ret; +} + DEFNDCMPLXCMPOP_FN (lt, complex_matrix, complex_matrix, complex_array, complex_array, mx_el_lt) DEFNDCMPLXCMPOP_FN (le, complex_matrix, complex_matrix, complex_array, complex_array, mx_el_le) DEFNDCMPLXCMPOP_FN (eq, complex_matrix, complex_matrix, complex_array, complex_array, mx_el_eq) @@ -199,6 +223,9 @@ INSTALL_BINOP (op_mul_trans, octave_complex_matrix, octave_complex_matrix, mul_trans); INSTALL_BINOP (op_herm_mul, octave_complex_matrix, octave_complex_matrix, herm_mul); INSTALL_BINOP (op_mul_herm, octave_complex_matrix, octave_complex_matrix, mul_herm); + INSTALL_BINOP (op_trans_ldiv, octave_complex_matrix, octave_complex_matrix, trans_ldiv); + INSTALL_BINOP (op_herm_ldiv, octave_complex_matrix, octave_complex_matrix, herm_ldiv); + INSTALL_BINOP (op_lt, octave_complex_matrix, octave_complex_matrix, lt); INSTALL_BINOP (op_le, octave_complex_matrix, octave_complex_matrix, le); INSTALL_BINOP (op_eq, octave_complex_matrix, octave_complex_matrix, eq);
--- a/src/OPERATORS/op-fcm-fcm.cc Mon Sep 21 14:24:27 2009 +0200 +++ b/src/OPERATORS/op-fcm-fcm.cc Wed Sep 23 10:00:16 2009 +0200 @@ -141,6 +141,32 @@ true, true, v2.float_complex_matrix_value ())); } +DEFBINOP (trans_ldiv, float_complex_matrix, float_complex_matrix) +{ + CAST_BINOP_ARGS (const octave_float_complex_matrix&, + const octave_float_complex_matrix&); + MatrixType typ = v1.matrix_type (); + + FloatComplexMatrix ret = xleftdiv (v1.float_complex_matrix_value (), + v2.float_complex_matrix_value (), typ, blas_trans); + + v1.matrix_type (typ); + return ret; +} + +DEFBINOP (herm_ldiv, float_complex_matrix, float_complex_matrix) +{ + CAST_BINOP_ARGS (const octave_float_complex_matrix&, + const octave_float_complex_matrix&); + MatrixType typ = v1.matrix_type (); + + FloatComplexMatrix ret = xleftdiv (v1.float_complex_matrix_value (), + v2.float_complex_matrix_value (), typ, blas_conj_trans); + + v1.matrix_type (typ); + return ret; +} + DEFNDCMPLXCMPOP_FN (lt, float_complex_matrix, float_complex_matrix, float_complex_array, float_complex_array, mx_el_lt) DEFNDCMPLXCMPOP_FN (le, float_complex_matrix, float_complex_matrix, @@ -239,6 +265,11 @@ octave_float_complex_matrix, herm_mul); INSTALL_BINOP (op_mul_herm, octave_float_complex_matrix, octave_float_complex_matrix, mul_herm); + INSTALL_BINOP (op_trans_ldiv, octave_float_complex_matrix, + octave_float_complex_matrix, trans_ldiv); + INSTALL_BINOP (op_herm_ldiv, octave_float_complex_matrix, + octave_float_complex_matrix, herm_ldiv); + INSTALL_BINOP (op_lt, octave_float_complex_matrix, octave_float_complex_matrix, lt); INSTALL_BINOP (op_le, octave_float_complex_matrix,
--- a/src/OPERATORS/op-fm-fm.cc Mon Sep 21 14:24:27 2009 +0200 +++ b/src/OPERATORS/op-fm-fm.cc Wed Sep 23 10:00:16 2009 +0200 @@ -110,6 +110,18 @@ true, v2.float_matrix_value ())); } +DEFBINOP (trans_ldiv, float_matrix, float_matrix) +{ + CAST_BINOP_ARGS (const octave_float_matrix&, const octave_float_matrix&); + MatrixType typ = v1.matrix_type (); + + FloatMatrix ret = xleftdiv (v1.float_matrix_value (), + v2.float_matrix_value (), typ, blas_trans); + + v1.matrix_type (typ); + return ret; +} + DEFNDBINOP_FN (lt, float_matrix, float_matrix, float_array, float_array, mx_el_lt) DEFNDBINOP_FN (le, float_matrix, float_matrix, float_array, @@ -217,6 +229,8 @@ INSTALL_BINOP (op_mul_trans, octave_float_matrix, octave_float_matrix, mul_trans); INSTALL_BINOP (op_herm_mul, octave_float_matrix, octave_float_matrix, trans_mul); INSTALL_BINOP (op_mul_herm, octave_float_matrix, octave_float_matrix, mul_trans); + INSTALL_BINOP (op_trans_ldiv, octave_float_matrix, octave_float_matrix, trans_ldiv); + INSTALL_BINOP (op_herm_ldiv, octave_float_matrix, octave_float_matrix, trans_ldiv); INSTALL_CATOP (octave_float_matrix, octave_float_matrix, fm_fm); INSTALL_CATOP (octave_matrix, octave_float_matrix, m_fm);
--- a/src/OPERATORS/op-m-m.cc Mon Sep 21 14:24:27 2009 +0200 +++ b/src/OPERATORS/op-m-m.cc Wed Sep 23 10:00:16 2009 +0200 @@ -106,6 +106,17 @@ return octave_value(xgemm (false, v1.matrix_value (), true, v2.matrix_value ())); } +DEFBINOP (trans_ldiv, matrix, matrix) +{ + CAST_BINOP_ARGS (const octave_matrix&, const octave_matrix&); + MatrixType typ = v1.matrix_type (); + + Matrix ret = xleftdiv (v1.matrix_value (), v2.matrix_value (), typ, blas_trans); + + v1.matrix_type (typ); + return ret; +} + DEFNDBINOP_FN (lt, matrix, matrix, array, array, mx_el_lt) DEFNDBINOP_FN (le, matrix, matrix, array, array, mx_el_le) DEFNDBINOP_FN (eq, matrix, matrix, array, array, mx_el_eq) @@ -190,6 +201,8 @@ INSTALL_BINOP (op_mul_trans, octave_matrix, octave_matrix, mul_trans); INSTALL_BINOP (op_herm_mul, octave_matrix, octave_matrix, trans_mul); INSTALL_BINOP (op_mul_herm, octave_matrix, octave_matrix, mul_trans); + INSTALL_BINOP (op_trans_ldiv, octave_matrix, octave_matrix, trans_ldiv); + INSTALL_BINOP (op_herm_ldiv, octave_matrix, octave_matrix, trans_ldiv); INSTALL_CATOP (octave_matrix, octave_matrix, m_m);
--- a/src/ov.cc Mon Sep 21 14:24:27 2009 +0200 +++ b/src/ov.cc Wed Sep 23 10:00:16 2009 +0200 @@ -377,6 +377,14 @@ retval = "timesherm"; break; + case op_trans_ldiv: + retval = "transldiv"; + break; + + case op_herm_ldiv: + retval = "hermldiv"; + break; + case op_el_and_not: retval = "andnot"; break; @@ -1996,6 +2004,16 @@ v1, do_unary_op (octave_value::op_hermitian, v2)); break; + case octave_value::op_trans_ldiv: + retval = do_binary_op (octave_value::op_ldiv, + do_unary_op (octave_value::op_transpose, v1), + v2); + break; + case octave_value::op_herm_ldiv: + retval = do_binary_op (octave_value::op_ldiv, + do_unary_op (octave_value::op_hermitian, v1), + v2); + break; case octave_value::op_el_not_and: retval = do_binary_op (octave_value::op_el_and, do_unary_op (octave_value::op_not, v1),
--- a/src/ov.h Mon Sep 21 14:24:27 2009 +0200 +++ b/src/ov.h Wed Sep 23 10:00:16 2009 +0200 @@ -112,6 +112,8 @@ op_mul_trans, op_herm_mul, op_mul_herm, + op_trans_ldiv, + op_herm_ldiv, op_el_not_and, op_el_not_or, op_el_and_not,
--- a/src/pt-cbinop.cc Mon Sep 21 14:24:27 2009 +0200 +++ b/src/pt-cbinop.cc Wed Sep 23 10:00:16 2009 +0200 @@ -106,6 +106,24 @@ return retop; } +// Possibly convert left division to trans_ldiv or herm_ldiv. + +static octave_value::compound_binary_op +simplify_ldiv_op (tree_expression *&a, tree_expression *&b) +{ + octave_value::compound_binary_op retop; + octave_value::unary_op opa = strip_trans_herm (a); + + if (opa == octave_value::op_hermitian) + retop = octave_value::op_herm_ldiv; + else if (opa == octave_value::op_transpose) + retop = octave_value::op_trans_ldiv; + else + retop = octave_value::unknown_compound_binary_op; + + return retop; +} + // Possibly contract and/or with negation. static octave_value::compound_binary_op @@ -152,6 +170,10 @@ ct = simplify_mul_op (ca, cb); break; + case octave_value::op_ldiv: + ct = simplify_ldiv_op (ca, cb); + break; + case octave_value::op_el_and: case octave_value::op_el_or: ct = simplify_and_or_op (ca, cb, t);
--- a/src/xdiv.cc Mon Sep 21 14:24:27 2009 +0200 +++ b/src/xdiv.cc Wed Sep 23 10:00:16 2009 +0200 @@ -132,17 +132,13 @@ if (! mx_div_conform (a, b)) return Matrix (); - Matrix atmp = a.transpose (); - Matrix btmp = b.transpose (); - MatrixType btyp = typ.transpose (); - octave_idx_type info; double rcond = 0.0; Matrix result - = btmp.solve (btyp, atmp, info, rcond, solve_singularity_warning); + = b.solve (typ, a.transpose (), info, rcond, + solve_singularity_warning, true, blas_trans); - typ = btyp.transpose (); return result.transpose (); } @@ -153,18 +149,14 @@ if (! mx_div_conform (a, b)) return ComplexMatrix (); - Matrix atmp = a.transpose (); - ComplexMatrix btmp = b.hermitian (); - MatrixType btyp = typ.transpose (); - octave_idx_type info; double rcond = 0.0; ComplexMatrix result - = btmp.solve (btyp, atmp, info, rcond, solve_singularity_warning); + = b.solve (typ, a.transpose (), info, rcond, + solve_singularity_warning, true, blas_trans); - typ = btyp.transpose (); - return result.hermitian (); + return result.transpose (); } // -*- 3 -*- @@ -174,18 +166,14 @@ if (! mx_div_conform (a, b)) return ComplexMatrix (); - ComplexMatrix atmp = a.hermitian (); - Matrix btmp = b.transpose (); - MatrixType btyp = typ.transpose (); - octave_idx_type info; double rcond = 0.0; ComplexMatrix result - = btmp.solve (btyp, atmp, info, rcond, solve_singularity_warning); + = b.solve (typ, a.transpose (), info, rcond, + solve_singularity_warning, true, blas_trans); - typ = btyp.transpose (); - return result.hermitian (); + return result.transpose (); } // -*- 4 -*- @@ -195,18 +183,14 @@ if (! mx_div_conform (a, b)) return ComplexMatrix (); - ComplexMatrix atmp = a.hermitian (); - ComplexMatrix btmp = b.hermitian (); - MatrixType btyp = typ.transpose (); - octave_idx_type info; double rcond = 0.0; ComplexMatrix result - = btmp.solve (btyp, atmp, info, rcond, solve_singularity_warning); + = b.solve (typ, a.transpose (), info, rcond, + solve_singularity_warning, true, blas_trans); - typ = btyp.transpose (); - return result.hermitian (); + return result.transpose (); } // Funny element by element division operations. @@ -366,19 +350,19 @@ // -*- 1 -*- Matrix -xleftdiv (const Matrix& a, const Matrix& b, MatrixType &typ) +xleftdiv (const Matrix& a, const Matrix& b, MatrixType &typ, blas_trans_type transt) { if (! mx_leftdiv_conform (a, b)) return Matrix (); octave_idx_type info; double rcond = 0.0; - return a.solve (typ, b, info, rcond, solve_singularity_warning); + return a.solve (typ, b, info, rcond, solve_singularity_warning, true, transt); } // -*- 2 -*- ComplexMatrix -xleftdiv (const Matrix& a, const ComplexMatrix& b, MatrixType &typ) +xleftdiv (const Matrix& a, const ComplexMatrix& b, MatrixType &typ, blas_trans_type transt) { if (! mx_leftdiv_conform (a, b)) return ComplexMatrix (); @@ -386,31 +370,31 @@ octave_idx_type info; double rcond = 0.0; - return a.solve (typ, b, info, rcond, solve_singularity_warning); + return a.solve (typ, b, info, rcond, solve_singularity_warning, true, transt); } // -*- 3 -*- ComplexMatrix -xleftdiv (const ComplexMatrix& a, const Matrix& b, MatrixType &typ) +xleftdiv (const ComplexMatrix& a, const Matrix& b, MatrixType &typ, blas_trans_type transt) { if (! mx_leftdiv_conform (a, b)) return ComplexMatrix (); octave_idx_type info; double rcond = 0.0; - return a.solve (typ, b, info, rcond, solve_singularity_warning); + return a.solve (typ, b, info, rcond, solve_singularity_warning, true, transt); } // -*- 4 -*- ComplexMatrix -xleftdiv (const ComplexMatrix& a, const ComplexMatrix& b, MatrixType &typ) +xleftdiv (const ComplexMatrix& a, const ComplexMatrix& b, MatrixType &typ, blas_trans_type transt) { if (! mx_leftdiv_conform (a, b)) return ComplexMatrix (); octave_idx_type info; double rcond = 0.0; - return a.solve (typ, b, info, rcond, solve_singularity_warning); + return a.solve (typ, b, info, rcond, solve_singularity_warning, true, transt); } static void @@ -446,17 +430,13 @@ if (! mx_div_conform (a, b)) return FloatMatrix (); - FloatMatrix atmp = a.transpose (); - FloatMatrix btmp = b.transpose (); - MatrixType btyp = typ.transpose (); - octave_idx_type info; float rcond = 0.0; FloatMatrix result - = btmp.solve (btyp, atmp, info, rcond, solve_singularity_warning); + = b.solve (typ, a.transpose (), info, rcond, + solve_singularity_warning, true, blas_trans); - typ = btyp.transpose (); return result.transpose (); } @@ -467,18 +447,14 @@ if (! mx_div_conform (a, b)) return FloatComplexMatrix (); - FloatMatrix atmp = a.transpose (); - FloatComplexMatrix btmp = b.hermitian (); - MatrixType btyp = typ.transpose (); - octave_idx_type info; float rcond = 0.0; - FloatComplexMatrix result - = btmp.solve (btyp, atmp, info, rcond, solve_singularity_warning); + FloatComplexMatrix result + = b.solve (typ, a.transpose (), info, rcond, + solve_singularity_warning, true, blas_trans); - typ = btyp.transpose (); - return result.hermitian (); + return result.transpose (); } // -*- 3 -*- @@ -488,18 +464,14 @@ if (! mx_div_conform (a, b)) return FloatComplexMatrix (); - FloatComplexMatrix atmp = a.hermitian (); - FloatMatrix btmp = b.transpose (); - MatrixType btyp = typ.transpose (); - octave_idx_type info; float rcond = 0.0; - FloatComplexMatrix result - = btmp.solve (btyp, atmp, info, rcond, solve_singularity_warning); + FloatComplexMatrix result + = b.solve (typ, a.transpose (), info, rcond, + solve_singularity_warning, true, blas_trans); - typ = btyp.transpose (); - return result.hermitian (); + return result.transpose (); } // -*- 4 -*- @@ -509,18 +481,14 @@ if (! mx_div_conform (a, b)) return FloatComplexMatrix (); - FloatComplexMatrix atmp = a.hermitian (); - FloatComplexMatrix btmp = b.hermitian (); - MatrixType btyp = typ.transpose (); - octave_idx_type info; float rcond = 0.0; - FloatComplexMatrix result - = btmp.solve (btyp, atmp, info, rcond, solve_singularity_warning); + FloatComplexMatrix result + = b.solve (typ, a.transpose (), info, rcond, + solve_singularity_warning, true, blas_trans); - typ = btyp.transpose (); - return result.hermitian (); + return result.transpose (); } // Funny element by element division operations. @@ -680,19 +648,19 @@ // -*- 1 -*- FloatMatrix -xleftdiv (const FloatMatrix& a, const FloatMatrix& b, MatrixType &typ) +xleftdiv (const FloatMatrix& a, const FloatMatrix& b, MatrixType &typ, blas_trans_type transt) { if (! mx_leftdiv_conform (a, b)) return FloatMatrix (); octave_idx_type info; float rcond = 0.0; - return a.solve (typ, b, info, rcond, solve_singularity_warning); + return a.solve (typ, b, info, rcond, solve_singularity_warning, true, transt); } // -*- 2 -*- FloatComplexMatrix -xleftdiv (const FloatMatrix& a, const FloatComplexMatrix& b, MatrixType &typ) +xleftdiv (const FloatMatrix& a, const FloatComplexMatrix& b, MatrixType &typ, blas_trans_type transt) { if (! mx_leftdiv_conform (a, b)) return FloatComplexMatrix (); @@ -700,31 +668,31 @@ octave_idx_type info; float rcond = 0.0; - return a.solve (typ, b, info, rcond, solve_singularity_warning); + return a.solve (typ, b, info, rcond, solve_singularity_warning, true, transt); } // -*- 3 -*- FloatComplexMatrix -xleftdiv (const FloatComplexMatrix& a, const FloatMatrix& b, MatrixType &typ) +xleftdiv (const FloatComplexMatrix& a, const FloatMatrix& b, MatrixType &typ, blas_trans_type transt) { if (! mx_leftdiv_conform (a, b)) return FloatComplexMatrix (); octave_idx_type info; float rcond = 0.0; - return a.solve (typ, b, info, rcond, solve_singularity_warning); + return a.solve (typ, b, info, rcond, solve_singularity_warning, true, transt); } // -*- 4 -*- FloatComplexMatrix -xleftdiv (const FloatComplexMatrix& a, const FloatComplexMatrix& b, MatrixType &typ) +xleftdiv (const FloatComplexMatrix& a, const FloatComplexMatrix& b, MatrixType &typ, blas_trans_type transt) { if (! mx_leftdiv_conform (a, b)) return FloatComplexMatrix (); octave_idx_type info; float rcond = 0.0; - return a.solve (typ, b, info, rcond, solve_singularity_warning); + return a.solve (typ, b, info, rcond, solve_singularity_warning, true, transt); } // Diagonal matrix division.
--- a/src/xdiv.h Mon Sep 21 14:24:27 2009 +0200 +++ b/src/xdiv.h Wed Sep 23 10:00:16 2009 +0200 @@ -25,15 +25,9 @@ #if !defined (octave_xdiv_h) #define octave_xdiv_h 1 -#include "oct-cmplx.h" +#include "mx-defs.h" #include "MatrixType.h" -class Matrix; -class ComplexMatrix; - -class NDArray; -class ComplexNDArray; - extern Matrix xdiv (const Matrix& a, const Matrix& b, MatrixType &typ); extern ComplexMatrix xdiv (const Matrix& a, const ComplexMatrix& b, MatrixType &typ); @@ -52,19 +46,14 @@ extern ComplexNDArray x_el_div (const Complex a, const NDArray& b); extern ComplexNDArray x_el_div (const Complex a, const ComplexNDArray& b); -extern Matrix xleftdiv (const Matrix& a, const Matrix& b, MatrixType &typ); +extern Matrix xleftdiv (const Matrix& a, const Matrix& b, MatrixType &typ, + blas_trans_type transt = blas_no_trans); extern ComplexMatrix xleftdiv (const Matrix& a, const ComplexMatrix& b, - MatrixType &typ); + MatrixType &typ, blas_trans_type transt = blas_no_trans); extern ComplexMatrix xleftdiv (const ComplexMatrix& a, const Matrix& b, - MatrixType &typ); + MatrixType &typ, blas_trans_type transt = blas_no_trans); extern ComplexMatrix xleftdiv (const ComplexMatrix& a, const ComplexMatrix& b, - MatrixType &typ); - -class FloatMatrix; -class FloatComplexMatrix; - -class FloatNDArray; -class FloatComplexNDArray; + MatrixType &typ, blas_trans_type transt = blas_no_trans); extern FloatMatrix xdiv (const FloatMatrix& a, const FloatMatrix& b, MatrixType &typ); extern FloatComplexMatrix xdiv (const FloatMatrix& a, const FloatComplexMatrix& b, @@ -84,19 +73,15 @@ extern FloatComplexNDArray x_el_div (const FloatComplex a, const FloatNDArray& b); extern FloatComplexNDArray x_el_div (const FloatComplex a, const FloatComplexNDArray& b); -extern FloatMatrix xleftdiv (const FloatMatrix& a, const FloatMatrix& b, MatrixType &typ); +extern FloatMatrix xleftdiv (const FloatMatrix& a, const FloatMatrix& b, MatrixType &typ, + blas_trans_type transt = blas_no_trans); extern FloatComplexMatrix xleftdiv (const FloatMatrix& a, const FloatComplexMatrix& b, - MatrixType &typ); + MatrixType &typ, blas_trans_type transt = blas_no_trans); extern FloatComplexMatrix xleftdiv (const FloatComplexMatrix& a, const FloatMatrix& b, - MatrixType &typ); + MatrixType &typ, blas_trans_type transt = blas_no_trans); extern FloatComplexMatrix xleftdiv (const FloatComplexMatrix& a, const FloatComplexMatrix& b, - MatrixType &typ); - + MatrixType &typ, blas_trans_type transt = blas_no_trans); -class DiagMatrix; -class FloatDiagMatrix; -class ComplexDiagMatrix; -class FloatComplexDiagMatrix; extern Matrix xdiv (const Matrix& a, const DiagMatrix& b); extern ComplexMatrix xdiv (const ComplexMatrix& a, const DiagMatrix& b);