Mercurial > octave-nkf
changeset 9662:0d3b248f4ab6
further improve mixed real-complex division
| author | Jaroslav Hajek <highegg@gmail.com> |
|---|---|
| date | Wed, 23 Sep 2009 11:10:52 +0200 |
| parents | afcf852256d2 |
| children | 7e5b4de5fbfe |
| files | liboctave/ChangeLog liboctave/dMatrix.cc liboctave/fMatrix.cc src/ChangeLog src/OPERATORS/op-fm-fcm.cc src/OPERATORS/op-m-cm.cc |
| diffstat | 6 files changed, 129 insertions(+), 16 deletions(-) [+] |
line wrap: on
line diff
--- a/liboctave/ChangeLog Wed Sep 23 10:00:16 2009 +0200 +++ b/liboctave/ChangeLog Wed Sep 23 11:10:52 2009 +0200 @@ -1,3 +1,14 @@ +2009-09-23 Jaroslav Hajek <highegg@gmail.com> + + * dMatrix.cc (stack_complex_matrix, unstack_complex_matrix): New + static funcs. + (Matrix::solve (..., const ComplexMatrix&, ...)): Use the above funcs. + Improve forwarding. + * fMatrix.cc (stack_complex_matrix, unstack_complex_matrix): New + static funcs. + (FloatMatrix::solve (..., const FloatComplexMatrix&, ...)): Use the + above funcs. Improve forwarding. + 2009-09-23 Jaroslav Hajek <highegg@gmail.com> * mx-defs.h (blas_trans_type): New enum.
--- a/liboctave/dMatrix.cc Wed Sep 23 10:00:16 2009 +0200 +++ b/liboctave/dMatrix.cc Wed Sep 23 11:10:52 2009 +0200 @@ -1911,6 +1911,13 @@ } Matrix +Matrix::solve (MatrixType &typ, const Matrix& b, octave_idx_type& info) const +{ + double rcon; + return solve (typ, b, info, rcon, 0); +} + +Matrix Matrix::solve (MatrixType &typ, const Matrix& b, octave_idx_type& info, double& rcon) const { @@ -1956,24 +1963,51 @@ ComplexMatrix Matrix::solve (MatrixType &typ, const ComplexMatrix& b) const { - ComplexMatrix tmp (*this); - return tmp.solve (typ, b); + octave_idx_type info; + double rcon; + return solve (typ, b, info, rcon, 0); } ComplexMatrix Matrix::solve (MatrixType &typ, const ComplexMatrix& b, octave_idx_type& info) const { - ComplexMatrix tmp (*this); - return tmp.solve (typ, b, info); + double rcon; + return solve (typ, b, info, rcon, 0); } ComplexMatrix Matrix::solve (MatrixType &typ, const ComplexMatrix& b, octave_idx_type& info, double& rcon) const { - ComplexMatrix tmp (*this); - return tmp.solve (typ, b, info, rcon); + return solve (typ, b, info, rcon, 0); +} + +static Matrix +stack_complex_matrix (const ComplexMatrix& cm) +{ + octave_idx_type m = cm.rows (), n = cm.cols (), nel = m*n; + Matrix retval (m, 2*n); + const Complex *cmd = cm.data (); + double *rd = retval.fortran_vec (); + for (octave_idx_type i = 0; i < nel; i++) + { + rd[i] = std::real (cmd[i]); + rd[nel+i] = std::imag (cmd[i]); + } + return retval; +} + +static ComplexMatrix +unstack_complex_matrix (const Matrix& sm) +{ + octave_idx_type m = sm.rows (), n = sm.cols () / 2, nel = m*n; + ComplexMatrix retval (m, n); + const double *smd = sm.data (); + Complex *rd = retval.fortran_vec (); + for (octave_idx_type i = 0; i < nel; i++) + rd[i] = Complex (smd[i], smd[nel+i]); + return retval; } ComplexMatrix @@ -1981,8 +2015,9 @@ double& rcon, solve_singularity_handler sing_handler, bool singular_fallback, blas_trans_type transt) const { - ComplexMatrix tmp (*this); - return tmp.solve (typ, b, info, rcon, sing_handler, singular_fallback, transt); + Matrix tmp = stack_complex_matrix (b); + tmp = solve (typ, tmp, info, rcon, sing_handler, singular_fallback, transt); + return unstack_complex_matrix (tmp); } ColumnVector
--- a/liboctave/fMatrix.cc Wed Sep 23 10:00:16 2009 +0200 +++ b/liboctave/fMatrix.cc Wed Sep 23 11:10:52 2009 +0200 @@ -1910,6 +1910,13 @@ } FloatMatrix +FloatMatrix::solve (MatrixType &typ, const FloatMatrix& b, octave_idx_type& info) const +{ + float rcon; + return solve (typ, b, info, rcon, 0); +} + +FloatMatrix FloatMatrix::solve (MatrixType &typ, const FloatMatrix& b, octave_idx_type& info, float& rcon) const { @@ -1955,24 +1962,51 @@ FloatComplexMatrix FloatMatrix::solve (MatrixType &typ, const FloatComplexMatrix& b) const { - FloatComplexMatrix tmp (*this); - return tmp.solve (typ, b); + octave_idx_type info; + float rcon; + return solve (typ, b, info, rcon, 0); } FloatComplexMatrix FloatMatrix::solve (MatrixType &typ, const FloatComplexMatrix& b, octave_idx_type& info) const { - FloatComplexMatrix tmp (*this); - return tmp.solve (typ, b, info); + float rcon; + return solve (typ, b, info, rcon, 0); } FloatComplexMatrix FloatMatrix::solve (MatrixType &typ, const FloatComplexMatrix& b, octave_idx_type& info, float& rcon) const { - FloatComplexMatrix tmp (*this); - return tmp.solve (typ, b, info, rcon); + return solve (typ, b, info, rcon, 0); +} + +static FloatMatrix +stack_complex_matrix (const FloatComplexMatrix& cm) +{ + octave_idx_type m = cm.rows (), n = cm.cols (), nel = m*n; + FloatMatrix retval (m, 2*n); + const FloatComplex *cmd = cm.data (); + float *rd = retval.fortran_vec (); + for (octave_idx_type i = 0; i < nel; i++) + { + rd[i] = std::real (cmd[i]); + rd[nel+i] = std::imag (cmd[i]); + } + return retval; +} + +static FloatComplexMatrix +unstack_complex_matrix (const FloatMatrix& sm) +{ + octave_idx_type m = sm.rows (), n = sm.cols () / 2, nel = m*n; + FloatComplexMatrix retval (m, n); + const float *smd = sm.data (); + FloatComplex *rd = retval.fortran_vec (); + for (octave_idx_type i = 0; i < nel; i++) + rd[i] = FloatComplex (smd[i], smd[nel+i]); + return retval; } FloatComplexMatrix @@ -1980,8 +2014,9 @@ float& rcon, solve_singularity_handler sing_handler, bool singular_fallback, blas_trans_type transt) const { - FloatComplexMatrix tmp (*this); - return tmp.solve (typ, b, info, rcon, sing_handler, singular_fallback, transt); + FloatMatrix tmp = stack_complex_matrix (b); + tmp = solve (typ, tmp, info, rcon, sing_handler, singular_fallback, transt); + return unstack_complex_matrix (tmp); } FloatColumnVector
--- a/src/ChangeLog Wed Sep 23 10:00:16 2009 +0200 +++ b/src/ChangeLog Wed Sep 23 11:10:52 2009 +0200 @@ -1,3 +1,8 @@ +2009-09-23 Jaroslav Hajek <highegg@gmail.com> + + * OPERATORS/op-m-cm.cc: Declare and install trans_ldiv operator. + * OPERATORS/op-fm-fcm.cc: Ditto. + 2009-09-23 Jaroslav Hajek <highegg@gmail.com> * ov.h (octave_value::op_trans_ldiv, op_herm_ldiv): New enum constants.
--- a/src/OPERATORS/op-fm-fcm.cc Wed Sep 23 10:00:16 2009 +0200 +++ b/src/OPERATORS/op-fm-fcm.cc Wed Sep 23 11:10:52 2009 +0200 @@ -83,6 +83,19 @@ return ret; } +DEFBINOP (trans_ldiv, float_matrix, float_complex_matrix) +{ + CAST_BINOP_ARGS (const octave_float_matrix&, + const octave_float_complex_matrix&); + MatrixType typ = v1.matrix_type (); + + FloatComplexMatrix ret = xleftdiv (v1.float_matrix_value (), + v2.float_complex_matrix_value (), typ, blas_trans); + + v1.matrix_type (typ); + return ret; +} + DEFNDCMPLXCMPOP_FN (lt, float_matrix, float_complex_matrix, float_array, float_complex_array, mx_el_lt) DEFNDCMPLXCMPOP_FN (le, float_matrix, float_complex_matrix, float_array,
--- a/src/OPERATORS/op-m-cm.cc Wed Sep 23 10:00:16 2009 +0200 +++ b/src/OPERATORS/op-m-cm.cc Wed Sep 23 11:10:52 2009 +0200 @@ -79,6 +79,18 @@ return ret; } +DEFBINOP (trans_ldiv, matrix, complex_matrix) +{ + CAST_BINOP_ARGS (const octave_matrix&, const octave_complex_matrix&); + MatrixType typ = v1.matrix_type (); + + ComplexMatrix ret = xleftdiv (v1.matrix_value (), + v2.complex_matrix_value (), typ, blas_trans); + + v1.matrix_type (typ); + return ret; +} + DEFNDCMPLXCMPOP_FN (lt, matrix, complex_matrix, array, complex_array, mx_el_lt) DEFNDCMPLXCMPOP_FN (le, matrix, complex_matrix, array, complex_array, mx_el_le) DEFNDCMPLXCMPOP_FN (eq, matrix, complex_matrix, array, complex_array, mx_el_eq) @@ -130,6 +142,8 @@ INSTALL_BINOP (op_el_ldiv, octave_matrix, octave_complex_matrix, el_ldiv); INSTALL_BINOP (op_el_and, octave_matrix, octave_complex_matrix, el_and); INSTALL_BINOP (op_el_or, octave_matrix, octave_complex_matrix, el_or); + INSTALL_BINOP (op_trans_ldiv, octave_matrix, octave_complex_matrix, trans_ldiv); + INSTALL_BINOP (op_herm_ldiv, octave_matrix, octave_complex_matrix, trans_ldiv); INSTALL_CATOP (octave_matrix, octave_complex_matrix, m_cm);