# HG changeset patch # User Jason Riedy # Date 1236635354 14400 # Node ID 1bba53c0a38dce5200737372289fd10ee9116fff # Parent 42aff15e059b8da1164d2e4bce7a58fb6f2df749 Implement diag + sparse, diag - sparse, sparse + diag, sparse - diag. Date: Mon, 9 Mar 2009 17:45:22 -0400 This does not use the typical sparse-mx-ops generator. I suspect the sematics of elementwise multiplication and division to be rather controversial, so they are not included. If comparison operations are added, the implementation should be shifted over to use the typical generator. The template in Sparse-diag-op-defs.h likely could use const bools rather than functional argument operations. I haven't measured which is optimized more effectively. Also, the Octave binding layer in op-dm-scm.cc likely could use all sorts of macro or template trickery, but it's far easier to let Emacs handle it for now. That would be worth revisiting if further elementwise sparse and diagonal operations are added. diff -r 42aff15e059b -r 1bba53c0a38d liboctave/CSparse.cc --- a/liboctave/CSparse.cc Mon Mar 09 17:49:14 2009 -0400 +++ b/liboctave/CSparse.cc Mon Mar 09 17:49:14 2009 -0400 @@ -7622,6 +7622,68 @@ return octave_impl::do_mul_sm_dm (a, d); } +SparseComplexMatrix +operator + (const ComplexDiagMatrix& d, const SparseMatrix& a) +{ + return octave_impl::do_add_dm_sm (d, a); +} +SparseComplexMatrix +operator + (const DiagMatrix& d, const SparseComplexMatrix& a) +{ + return octave_impl::do_add_dm_sm (d, a); +} +SparseComplexMatrix +operator + (const ComplexDiagMatrix& d, const SparseComplexMatrix& a) +{ + return octave_impl::do_add_dm_sm (d, a); +} +SparseComplexMatrix +operator + (const SparseMatrix& a, const ComplexDiagMatrix& d) +{ + return octave_impl::do_add_sm_dm (a, d); +} +SparseComplexMatrix +operator + (const SparseComplexMatrix& a, const DiagMatrix& d) +{ + return octave_impl::do_add_sm_dm (a, d); +} +SparseComplexMatrix +operator + (const SparseComplexMatrix&a, const ComplexDiagMatrix& d) +{ + return octave_impl::do_add_sm_dm (a, d); +} + +SparseComplexMatrix +operator - (const ComplexDiagMatrix& d, const SparseMatrix& a) +{ + return octave_impl::do_sub_dm_sm (d, a); +} +SparseComplexMatrix +operator - (const DiagMatrix& d, const SparseComplexMatrix& a) +{ + return octave_impl::do_sub_dm_sm (d, a); +} +SparseComplexMatrix +operator - (const ComplexDiagMatrix& d, const SparseComplexMatrix& a) +{ + return octave_impl::do_sub_dm_sm (d, a); +} +SparseComplexMatrix +operator - (const SparseMatrix& a, const ComplexDiagMatrix& d) +{ + return octave_impl::do_sub_sm_dm (a, d); +} +SparseComplexMatrix +operator - (const SparseComplexMatrix& a, const DiagMatrix& d) +{ + return octave_impl::do_sub_sm_dm (a, d); +} +SparseComplexMatrix +operator - (const SparseComplexMatrix&a, const ComplexDiagMatrix& d) +{ + return octave_impl::do_sub_sm_dm (a, d); +} + // FIXME -- it would be nice to share code among the min/max // functions below. diff -r 42aff15e059b -r 1bba53c0a38d liboctave/CSparse.h --- a/liboctave/CSparse.h Mon Mar 09 17:49:14 2009 -0400 +++ b/liboctave/CSparse.h Mon Mar 09 17:49:14 2009 -0400 @@ -484,6 +484,20 @@ extern OCTAVE_API SparseComplexMatrix operator * (const ComplexDiagMatrix&, const SparseComplexMatrix&); extern OCTAVE_API SparseComplexMatrix operator * (const SparseComplexMatrix&, const ComplexDiagMatrix&); +extern OCTAVE_API SparseComplexMatrix operator + (const ComplexDiagMatrix&, const SparseMatrix&); +extern OCTAVE_API SparseComplexMatrix operator + (const DiagMatrix&, const SparseComplexMatrix&); +extern OCTAVE_API SparseComplexMatrix operator + (const ComplexDiagMatrix&, const SparseComplexMatrix&); +extern OCTAVE_API SparseComplexMatrix operator + (const SparseMatrix&, const ComplexDiagMatrix&); +extern OCTAVE_API SparseComplexMatrix operator + (const SparseComplexMatrix&, const DiagMatrix&); +extern OCTAVE_API SparseComplexMatrix operator + (const SparseComplexMatrix&, const ComplexDiagMatrix&); + +extern OCTAVE_API SparseComplexMatrix operator - (const ComplexDiagMatrix&, const SparseMatrix&); +extern OCTAVE_API SparseComplexMatrix operator - (const DiagMatrix&, const SparseComplexMatrix&); +extern OCTAVE_API SparseComplexMatrix operator - (const ComplexDiagMatrix&, const SparseComplexMatrix&); +extern OCTAVE_API SparseComplexMatrix operator - (const SparseMatrix&, const ComplexDiagMatrix&); +extern OCTAVE_API SparseComplexMatrix operator - (const SparseComplexMatrix&, const DiagMatrix&); +extern OCTAVE_API SparseComplexMatrix operator - (const SparseComplexMatrix&, const ComplexDiagMatrix&); + extern OCTAVE_API SparseComplexMatrix min (const Complex& c, const SparseComplexMatrix& m); extern OCTAVE_API SparseComplexMatrix min (const SparseComplexMatrix& m, diff -r 42aff15e059b -r 1bba53c0a38d liboctave/ChangeLog --- a/liboctave/ChangeLog Mon Mar 09 17:49:14 2009 -0400 +++ b/liboctave/ChangeLog Mon Mar 09 17:49:14 2009 -0400 @@ -1,3 +1,40 @@ +2009-03-09 Jason Riedy + + * Sparse-diag-op-defs.h (octave_impl::inner_do_add_sm_dm): New + template function. Implementation for adding sparse and diagonal + matrices. Takes two functional arguments, opa and opd, to + generate both subtraction variants. + (octave_impl::do_commutative_add_dm_sm): New template function. + Ensure A+D and D+A use the same generated code. + (octave_impl::do_add_dm_sm): New template function. Check + arguments for diag + sparse and call inner routine. + (octave_impl::do_sub_dm_sm): New template function. Check + arguments for diag - sparse and call inner routine. + (octave_impl::do_add_sm_dm): New template function. Check + arguments for sparse + diag and call inner routine. + (octave_impl::do_sub_sm_dm): New template function. Check + arguments for sparse - diag and call inner routine. + + * dSparse.h (operator +): Declare overrides for real diag + + sparse. + (operator -): Declare overrides for real diag - sparse, sparse - + diag. + + * dSparse.cc (operator +): Define overrides for real diag + + sparse. + (operator -): Define overrides for real diag - sparse, sparse - + diag. + + * CSparse.h (operator +): Declare overrides for complex and real + combinations of diag + sparse. + (operator -): Declare overrides for complex and real combinations + of diag - sparse, sparse - diag. + + * CSparse.cc (operator +): Define overrides for complex and real + combinations of diag + sparse. + (operator -): Define overrides for complex and real combinations + of diag - sparse, sparse - diag. + 2009-03-08 Jason Riedy * Sparse-diag-op-defs.h (octave_impl::do_mul_dm_sm) diff -r 42aff15e059b -r 1bba53c0a38d liboctave/Sparse-diag-op-defs.h --- a/liboctave/Sparse-diag-op-defs.h Mon Mar 09 17:49:14 2009 -0400 +++ b/liboctave/Sparse-diag-op-defs.h Mon Mar 09 17:49:14 2009 -0400 @@ -23,6 +23,8 @@ #if !defined (octave_sparse_diag_op_defs_h) #define octave_sparse_diag_op_defs_h 1 +// Matrix multiplication + template RT do_mul_dm_sm (const DM& d, const SM& a) { @@ -103,4 +105,122 @@ } } +// FIXME: functors such as this should be gathered somewhere +template +struct identity_val + : public std::unary_function +{ + T operator () (const T x) { return x; } +}; + +// Matrix addition + +template +RT inner_do_add_sm_dm (const SM& a, const DM& d, OpA opa, OpD opd) +{ + using std::min; + const octave_idx_type nr = d.rows (); + const octave_idx_type nc = d.cols (); + const octave_idx_type n = min (nr, nc); + + const octave_idx_type a_nr = a.rows (); + const octave_idx_type a_nc = a.cols (); + + const octave_idx_type nz = a.nnz (); + const typename SM::element_type zero = typename SM::element_type (); + const typename DM::element_type dzero = typename DM::element_type (); + RT r (a_nr, a_nc, nz + n); + octave_idx_type k = 0; + + for (octave_idx_type j = 0; j < nc; ++j) + { + OCTAVE_QUIT; + const octave_idx_type colend = a.cidx (j+1); + bool found_diag = false; + r.xcidx (j) = k; + for (octave_idx_type k_src = a.cidx (j); k_src < colend; ++k_src, ++k) + { + const octave_idx_type i = a.ridx (k_src); + r.xridx (k) = i; + if (i != j) + r.xdata (k) = opa (a.data (k_src)); + else + { + r.xdata (k) = opa (a.data (k_src)) + opd (d.dgelem (j)); + found_diag = true; + } + } + if (!found_diag) + { + r.xridx (k) = j; + r.xdata (k) = opd (d.dgelem (j)); + ++k; + } + } + r.xcidx (nc) = k; + + r.maybe_compress (true); + return r; +} + +template +RT do_commutative_add_dm_sm (const DM& d, const SM& a) +{ + // Extra function to ensure this is only emitted once. + return inner_do_add_sm_dm (a, d, + identity_val (), + identity_val ()); +} + +template +RT do_add_dm_sm (const DM& d, const SM& a) +{ + if (a.rows () != d.rows () || a.cols () != d.cols ()) + { + gripe_nonconformant ("operator +", d.rows (), d.cols (), a.rows (), a.cols ()); + return RT (); + } + else + return do_commutative_add_dm_sm (d, a); +} + +template +RT do_sub_dm_sm (const DM& d, const SM& a) +{ + if (a.rows () != d.rows () || a.cols () != d.cols ()) + { + gripe_nonconformant ("operator -", d.rows (), d.cols (), a.rows (), a.cols ()); + return RT (); + } + else + return inner_do_add_sm_dm (a, d, std::negate (), + identity_val ()); +} + +template +RT do_add_sm_dm (const SM& a, const DM& d) +{ + if (a.rows () != d.rows () || a.cols () != d.cols ()) + { + gripe_nonconformant ("operator +", a.rows (), a.cols (), d.rows (), d.cols ()); + return RT (); + } + else + return do_commutative_add_dm_sm (d, a); +} + +template +RT do_sub_sm_dm (const SM& a, const DM& d) +{ + if (a.rows () != d.rows () || a.cols () != d.cols ()) + { + gripe_nonconformant ("operator -", a.rows (), a.cols (), d.rows (), d.cols ()); + return RT (); + } + else + return inner_do_add_sm_dm (a, d, + identity_val (), + std::negate ()); +} + #endif // octave_sparse_diag_op_defs_h diff -r 42aff15e059b -r 1bba53c0a38d liboctave/dSparse.cc --- a/liboctave/dSparse.cc Mon Mar 09 17:49:14 2009 -0400 +++ b/liboctave/dSparse.cc Mon Mar 09 17:49:14 2009 -0400 @@ -7716,6 +7716,30 @@ return octave_impl::do_mul_sm_dm (a, d); } +SparseMatrix +operator + (const DiagMatrix& d, const SparseMatrix& a) +{ + return octave_impl::do_add_dm_sm (d, a); +} + +SparseMatrix +operator - (const DiagMatrix& d, const SparseMatrix& a) +{ + return octave_impl::do_sub_dm_sm (d, a); +} + +SparseMatrix +operator + (const SparseMatrix& a, const DiagMatrix& d) +{ + return octave_impl::do_add_sm_dm (a, d); +} + +SparseMatrix +operator - (const SparseMatrix& a, const DiagMatrix& d) +{ + return octave_impl::do_sub_sm_dm (a, d); +} + // FIXME -- it would be nice to share code among the min/max // functions below. diff -r 42aff15e059b -r 1bba53c0a38d liboctave/dSparse.h --- a/liboctave/dSparse.h Mon Mar 09 17:49:14 2009 -0400 +++ b/liboctave/dSparse.h Mon Mar 09 17:49:14 2009 -0400 @@ -454,6 +454,11 @@ extern OCTAVE_API SparseMatrix operator * (const DiagMatrix&, const SparseMatrix&); extern OCTAVE_API SparseMatrix operator * (const SparseMatrix&, const DiagMatrix&); +extern OCTAVE_API SparseMatrix operator + (const DiagMatrix&, const SparseMatrix&); +extern OCTAVE_API SparseMatrix operator + (const SparseMatrix&, const DiagMatrix&); +extern OCTAVE_API SparseMatrix operator - (const DiagMatrix&, const SparseMatrix&); +extern OCTAVE_API SparseMatrix operator - (const SparseMatrix&, const DiagMatrix&); + extern OCTAVE_API SparseMatrix min (double d, const SparseMatrix& m); extern OCTAVE_API SparseMatrix min (const SparseMatrix& m, double d); extern OCTAVE_API SparseMatrix min (const SparseMatrix& a, const SparseMatrix& b); diff -r 42aff15e059b -r 1bba53c0a38d src/ChangeLog --- a/src/ChangeLog Mon Mar 09 17:49:14 2009 -0400 +++ b/src/ChangeLog Mon Mar 09 17:49:14 2009 -0400 @@ -99,6 +99,28 @@ * Makefile.in (maintainer-clean): Remove y.tab.h here. (distclean): Not here. +2009-03-09 Jason Riedy + + * OPERATORS/op-dm-sm.cc (add_dm_sm): Octave binding for diag + sparse. + (sub_dm_sm): Octave binding for diag - sparse. + (add_sm_dm): Octave binding for diag + sparse. + (sub_sm_dm): Octave binding for diag - sparse. + (install_dm_sm_ops): Install above bindings. + + * OPERATORS/op-dm-scm.cc (add_cdm_sm): Octave binding for diag + sparse. + (add_dm_scm): Octave binding for diag + sparse. + (add_cdm_scm): Octave binding for diag + sparse. + (sub_cdm_sm): Octave binding for diag - sparse. + (sub_dm_scm): Octave binding for diag - sparse. + (sub_cdm_csm): Octave binding for diag - sparse. + (add_sm_cdm): Octave binding for diag + sparse. + (add_scm_dm): Octave binding for diag + sparse. + (add_scm_cdm): Octave binding for diag + sparse. + (sub_sm_cdm): Octave binding for diag - sparse. + (sub_scm_dm): Octave binding for diag - sparse. + (sub_scm_cdm): Octave binding for diag - sparse. + (install_dm_scm_ops): Install above bindings. + 2009-03-08 Jason Riedy * sparse-xdiv.h (xleftdiv): Declare overrides for diff -r 42aff15e059b -r 1bba53c0a38d src/OPERATORS/op-dm-scm.cc --- a/src/OPERATORS/op-dm-scm.cc Mon Mar 09 17:49:14 2009 -0400 +++ b/src/OPERATORS/op-dm-scm.cc Mon Mar 09 17:49:14 2009 -0400 @@ -1,6 +1,6 @@ /* -Copyright (C) 2009 Jason Riedy +Copyright (C) 2009 Jason Riedy, Jaroslav Hajek This file is part of Octave. @@ -138,6 +138,102 @@ typ); } +DEFBINOP (add_dm_scm, diag_matrix, sparse_complex_matrix) +{ + CAST_BINOP_ARGS (const octave_diag_matrix&, const octave_sparse_complex_matrix&); + + if (v2.rows() == 1 && v2.columns() == 1) + // If v2 is a scalar in disguise, return a diagonal matrix rather than + // a sparse matrix. + { + std::complex d = v2.complex_value (); + + return octave_value (v1.diag_matrix_value () + d); + } + else + return v1.diag_matrix_value () + v2.sparse_complex_matrix_value (); +} + +DEFBINOP (add_cdm_sm, complex_diag_matrix, sparse_matrix) +{ + CAST_BINOP_ARGS (const octave_complex_diag_matrix&, const octave_sparse_matrix&); + + if (v2.rows() == 1 && v2.columns() == 1) + // If v2 is a scalar in disguise, return a diagonal matrix rather than + // a sparse matrix. + { + double d = v2.scalar_value (); + + return octave_value (v1.complex_diag_matrix_value () + d); + } + else + return v1.complex_diag_matrix_value () + v2.sparse_matrix_value (); +} + +DEFBINOP (add_cdm_scm, complex_diag_matrix, sparse_complex_matrix) +{ + CAST_BINOP_ARGS (const octave_complex_diag_matrix&, const octave_sparse_complex_matrix&); + + if (v2.rows() == 1 && v2.columns() == 1) + // If v2 is a scalar in disguise, return a diagonal matrix rather than + // a sparse matrix. + { + std::complex d = v2.complex_value (); + + return octave_value (v1.complex_diag_matrix_value () + d); + } + else + return v1.complex_diag_matrix_value () + v2.sparse_complex_matrix_value (); +} + +DEFBINOP (sub_dm_scm, diag_matrix, sparse_complex_matrix) +{ + CAST_BINOP_ARGS (const octave_diag_matrix&, const octave_sparse_complex_matrix&); + + if (v2.rows() == 1 && v2.columns() == 1) + // If v2 is a scalar in disguise, return a diagonal matrix rather than + // a sparse matrix. + { + std::complex d = v2.complex_value (); + + return octave_value (v1.diag_matrix_value () + (-d)); + } + else + return v1.diag_matrix_value () - v2.sparse_complex_matrix_value (); +} + +DEFBINOP (sub_cdm_sm, complex_diag_matrix, sparse_matrix) +{ + CAST_BINOP_ARGS (const octave_complex_diag_matrix&, const octave_sparse_matrix&); + + if (v2.rows() == 1 && v2.columns() == 1) + // If v2 is a scalar in disguise, return a diagonal matrix rather than + // a sparse matrix. + { + double d = v2.scalar_value (); + + return octave_value (v1.complex_diag_matrix_value () + (-d)); + } + else + return v1.complex_diag_matrix_value () - v2.sparse_matrix_value (); +} + +DEFBINOP (sub_cdm_scm, complex_diag_matrix, sparse_complex_matrix) +{ + CAST_BINOP_ARGS (const octave_complex_diag_matrix&, const octave_sparse_complex_matrix&); + + if (v2.rows() == 1 && v2.columns() == 1) + // If v2 is a scalar in disguise, return a diagonal matrix rather than + // a sparse matrix. + { + std::complex d = v2.complex_value (); + + return octave_value (v1.complex_diag_matrix_value () + (-d)); + } + else + return v1.complex_diag_matrix_value () - v2.sparse_complex_matrix_value (); +} + // sparse matrix by diagonal matrix ops DEFBINOP (mul_scm_dm, sparse_complex_matrix, diag_matrix) @@ -276,6 +372,102 @@ } } +DEFBINOP (add_sm_cdm, sparse_matrix, complex_diag_matrix) +{ + CAST_BINOP_ARGS (const octave_sparse_matrix&, const octave_complex_diag_matrix&); + + if (v2.rows() == 1 && v2.columns() == 1) + // If v2 is a scalar in disguise, return a diagonal matrix rather than + // a sparse matrix. + { + std::complex d = v2.complex_value (); + + return octave_value (v1.sparse_matrix_value () + d); + } + else + return v1.sparse_matrix_value () + v2.complex_diag_matrix_value (); +} + +DEFBINOP (add_scm_dm, sparse_complex_matrix, diag_matrix) +{ + CAST_BINOP_ARGS (const octave_sparse_complex_matrix&, const octave_diag_matrix&); + + if (v2.rows() == 1 && v2.columns() == 1) + // If v2 is a scalar in disguise, return a diagonal matrix rather than + // a sparse matrix. + { + double d = v2.scalar_value (); + + return octave_value (v1.sparse_complex_matrix_value () + d); + } + else + return v1.sparse_complex_matrix_value () + v2.diag_matrix_value (); +} + +DEFBINOP (add_scm_cdm, sparse_complex_matrix, complex_diag_matrix) +{ + CAST_BINOP_ARGS (const octave_sparse_complex_matrix&, const octave_complex_diag_matrix&); + + if (v2.rows() == 1 && v2.columns() == 1) + // If v2 is a scalar in disguise, return a diagonal matrix rather than + // a sparse matrix. + { + std::complex d = v2.complex_value (); + + return octave_value (v1.sparse_complex_matrix_value () + d); + } + else + return v1.sparse_complex_matrix_value () + v2.complex_diag_matrix_value (); +} + +DEFBINOP (sub_sm_cdm, sparse_matrix, complex_diag_matrix) +{ + CAST_BINOP_ARGS (const octave_sparse_matrix&, const octave_complex_diag_matrix&); + + if (v2.rows() == 1 && v2.columns() == 1) + // If v2 is a scalar in disguise, return a diagonal matrix rather than + // a sparse matrix. + { + std::complex d = v2.complex_value (); + + return octave_value (v1.sparse_matrix_value () + (-d)); + } + else + return v1.sparse_matrix_value () - v2.complex_diag_matrix_value (); +} + +DEFBINOP (sub_scm_dm, sparse_complex_matrix, diag_matrix) +{ + CAST_BINOP_ARGS (const octave_sparse_complex_matrix&, const octave_diag_matrix&); + + if (v2.rows() == 1 && v2.columns() == 1) + // If v2 is a scalar in disguise, return a diagonal matrix rather than + // a sparse matrix. + { + double d = v2.scalar_value (); + + return octave_value (v1.sparse_complex_matrix_value () + (-d)); + } + else + return v1.sparse_complex_matrix_value () - v2.diag_matrix_value (); +} + +DEFBINOP (sub_scm_cdm, sparse_complex_matrix, complex_diag_matrix) +{ + CAST_BINOP_ARGS (const octave_sparse_complex_matrix&, const octave_complex_diag_matrix&); + + if (v2.rows() == 1 && v2.columns() == 1) + // If v2 is a scalar in disguise, return a diagonal matrix rather than + // a sparse matrix. + { + std::complex d = v2.complex_value (); + + return octave_value (v1.sparse_complex_matrix_value () + (-d)); + } + else + return v1.sparse_complex_matrix_value () - v2.complex_diag_matrix_value (); +} + void install_dm_scm_ops (void) { @@ -290,6 +482,15 @@ INSTALL_BINOP (op_ldiv, octave_complex_diag_matrix, octave_sparse_complex_matrix, ldiv_cdm_scm); + INSTALL_BINOP (op_add, octave_diag_matrix, octave_sparse_complex_matrix, add_dm_scm); + INSTALL_BINOP (op_add, octave_complex_diag_matrix, octave_sparse_matrix, add_cdm_sm); + INSTALL_BINOP (op_add, octave_complex_diag_matrix, octave_sparse_complex_matrix, + add_cdm_scm); + INSTALL_BINOP (op_sub, octave_diag_matrix, octave_sparse_complex_matrix, sub_dm_scm); + INSTALL_BINOP (op_sub, octave_complex_diag_matrix, octave_sparse_matrix, sub_cdm_sm); + INSTALL_BINOP (op_sub, octave_complex_diag_matrix, octave_sparse_complex_matrix, + sub_cdm_scm); + INSTALL_BINOP (op_mul, octave_sparse_complex_matrix, octave_diag_matrix, mul_scm_dm); INSTALL_BINOP (op_mul, octave_sparse_matrix, octave_complex_diag_matrix, @@ -300,4 +501,11 @@ INSTALL_BINOP (op_div, octave_sparse_complex_matrix, octave_diag_matrix, div_scm_dm); INSTALL_BINOP (op_div, octave_sparse_matrix, octave_complex_diag_matrix, div_sm_cdm); INSTALL_BINOP (op_div, octave_sparse_complex_matrix, octave_complex_diag_matrix, div_scm_cdm); + + INSTALL_BINOP (op_add, octave_sparse_complex_matrix, octave_diag_matrix, add_scm_dm); + INSTALL_BINOP (op_add, octave_sparse_matrix, octave_complex_diag_matrix, add_sm_cdm); + INSTALL_BINOP (op_add, octave_sparse_complex_matrix, octave_complex_diag_matrix, add_scm_cdm); + INSTALL_BINOP (op_sub, octave_sparse_complex_matrix, octave_diag_matrix, sub_scm_dm); + INSTALL_BINOP (op_sub, octave_sparse_matrix, octave_complex_diag_matrix, sub_sm_cdm); + INSTALL_BINOP (op_sub, octave_sparse_complex_matrix, octave_complex_diag_matrix, sub_scm_cdm); } diff -r 42aff15e059b -r 1bba53c0a38d src/OPERATORS/op-dm-sm.cc --- a/src/OPERATORS/op-dm-sm.cc Mon Mar 09 17:49:14 2009 -0400 +++ b/src/OPERATORS/op-dm-sm.cc Mon Mar 09 17:49:14 2009 -0400 @@ -1,6 +1,6 @@ /* -Copyright (C) 2009 Jason Riedy +Copyright (C) 2009 Jason Riedy, Jaroslav Hajek This file is part of Octave. @@ -68,6 +68,38 @@ return xleftdiv (v1.diag_matrix_value (), v2.sparse_matrix_value (), typ); } +DEFBINOP_OP (add_dm_sm, diag_matrix, sparse_matrix) +{ + CAST_BINOP_ARGS (const octave_diag_matrix&, const octave_sparse_matrix&); + + if (v2.rows() == 1 && v2.columns() == 1) + // If v2 is a scalar in disguise, return a diagonal matrix rather than + // a sparse matrix. + { + double d = v2.scalar_value (); + + return octave_value (v1.diag_matrix_value () + d); + } + else + return v1.diag_matrix_value () + v2.sparse_matrix_value (); +} + +DEFBINOP (sub_dm_sm, diag_matrix, sparse_matrix) +{ + CAST_BINOP_ARGS (const octave_diag_matrix&, const octave_sparse_matrix&); + + if (v2.rows() == 1 && v2.columns() == 1) + // If v2 is a scalar in disguise, return a diagonal matrix rather than + // a sparse matrix. + { + double d = v2.scalar_value (); + + return octave_value (v1.diag_matrix_value () - d); + } + else + return v1.diag_matrix_value () - v2.sparse_matrix_value (); +} + // sparse matrix by diagonal matrix ops DEFBINOP (mul_sm_dm, sparse_matrix, diag_matrix) @@ -113,16 +145,52 @@ } } +DEFBINOP (add_sm_dm, sparse_matrix, diag_matrix) +{ + CAST_BINOP_ARGS (const octave_sparse_matrix&, const octave_diag_matrix&); + + if (v1.rows() == 1 && v1.columns() == 1) + // If v1 is a scalar in disguise, return a diagonal matrix rather than + // a sparse matrix. + { + double d = v1.scalar_value (); + + return octave_value (d + v2.diag_matrix_value ()); + } + else + return v1.sparse_matrix_value () + v2.diag_matrix_value (); +} + +DEFBINOP (sub_sm_dm, sparse_matrix, diag_matrix) +{ + CAST_BINOP_ARGS (const octave_sparse_matrix&, const octave_diag_matrix&); + + if (v1.rows() == 1 && v1.columns() == 1) + // If v1 is a scalar in disguise, return a diagonal matrix rather than + // a sparse matrix. + { + double d = v1.scalar_value (); + + return octave_value (d - v2.diag_matrix_value ()); + } + else + return v1.sparse_matrix_value () - v2.diag_matrix_value (); +} + void install_dm_sm_ops (void) { INSTALL_BINOP (op_mul, octave_diag_matrix, octave_sparse_matrix, mul_dm_sm); + INSTALL_BINOP (op_add, octave_diag_matrix, octave_sparse_matrix, add_dm_sm); + INSTALL_BINOP (op_sub, octave_diag_matrix, octave_sparse_matrix, sub_dm_sm); INSTALL_BINOP (op_ldiv, octave_diag_matrix, octave_sparse_matrix, ldiv_dm_sm); INSTALL_BINOP (op_mul, octave_sparse_matrix, octave_diag_matrix, mul_sm_dm); + INSTALL_BINOP (op_add, octave_sparse_matrix, octave_diag_matrix, add_sm_dm); + INSTALL_BINOP (op_sub, octave_sparse_matrix, octave_diag_matrix, sub_sm_dm); INSTALL_BINOP (op_div, octave_sparse_matrix, octave_diag_matrix, div_sm_dm); } diff -r 42aff15e059b -r 1bba53c0a38d test/ChangeLog --- a/test/ChangeLog Mon Mar 09 17:49:14 2009 -0400 +++ b/test/ChangeLog Mon Mar 09 17:49:14 2009 -0400 @@ -6,6 +6,10 @@ * test_diag_perm.m: Add a test for conversion to sparse form. +2009-03-09 Jason Riedy + + * test_diag_perm.m: Add tests for diag + sparse. + 2009-03-08 Jason Riedy * test_diag_perm.m: Add tests for inverse scaling and sparse structure. diff -r 42aff15e059b -r 1bba53c0a38d test/test_diag_perm.m --- a/test/test_diag_perm.m Mon Mar 09 17:49:14 2009 -0400 +++ b/test/test_diag_perm.m Mon Mar 09 17:49:14 2009 -0400 @@ -192,3 +192,28 @@ %! scalefact = rand (1, n-2) + I () * rand(1, n-2); %! Dc = diag (scalefact, n-2, n); %! assert (full (A / Dc), full(A) / Dc) + +## adding sparse and diagonal stays sparse +%!test +%! n = 9; +%! A = sprand (n, n, .5); +%! D = 2 * eye (n); +%! assert (typeinfo (A + D), "sparse matrix") +%! assert (typeinfo (A - D), "sparse matrix") +%! D = D * I () + D; +%! assert (typeinfo (A - D), "sparse complex matrix") +%! A = A * I () + A; +%! assert (typeinfo (D - A), "sparse complex matrix") + +## adding sparse and diagonal stays sparse +%!test +%! n = 9; +%! A = sprand (n, n, .5); +%! D = 2 * eye (n); +%! assert (full (A + D), full (A) + D) +%! assert (full (A - D), full (A) - D) +%! D = D * I () + D; +%! assert (full (D + A), D + full (A)) +%! A = A * I () + A; +%! A(6, 4) = nan (); +%! assert (full (D - A), D - full (A))