# HG changeset patch # User Jason Riedy # Date 1236635354 14400 # Node ID 42aff15e059b8da1164d2e4bce7a58fb6f2df749 # Parent f4f4d65faaa0d95b35109cf63a24f6e49e528aed Implement diag \ sparse and sparse / diag. Not pretty, but somewhat efficient and preserves sparsity. diff -r f4f4d65faaa0 -r 42aff15e059b src/ChangeLog --- a/src/ChangeLog Mon Mar 09 17:49:13 2009 -0400 +++ b/src/ChangeLog Mon Mar 09 17:49:14 2009 -0400 @@ -101,6 +101,37 @@ 2009-03-08 Jason Riedy + * sparse-xdiv.h (xleftdiv): Declare overrides for + xleftdiv (diagonal, sparse), both real and complex. + (xdiv): Declare overrides for xdiv (sparse, diagonal), both real + and complex. + + * sparse-xdiv.cc (do_rightdiv_sm_dm): New template function. + Implementation for xdiv (sparse, diagonal). + (xdiv): Call do_rightdiv_sm_dm to implement overrides, real and + complex. + (do_leftdiv_dm_sm): New template function. Implementation for + xleftdiv (diagonal, sparse). + (xleftdiv): Call do_leftdiv_dm_sm to implement overrides, real and + complex. + + * OPERATORS/op-dm-sm.cc (ldiv_dm_sm): Octave binding for real left + division, diagonal into sparse. + (div_sm_dm): Octave binding for real right division, sparse by + diagonal. + (install_dm_sm_ops): Install above bindings. + + * OPERATORS/op-dm-scm.cc (ldiv_dm_scm): Octave binding for real + diagonal \ complex sparse. + (ldiv_cdm_sm): Octave binding for complex diagonal \ real sparse. + (ldiv_cdm_scm): Octave binding for complex diagonal \ complex sparse. + (div_scm_dm): Octave binding for complex sparse / real diagonal. + (div_sm_cdm): Octave binding for real sparse / complex diagonal. + (div_scm_cdm): Octave binding for complex sparse / complex diagonal. + (install_dm_scm_ops): Install above bindings. + +2009-03-08 Jason Riedy + * OPERATORS/op-dm-scm.cc: New file. Implement multiplication between diagonal matrices (both real and complex) and complex sparse matrices. diff -r f4f4d65faaa0 -r 42aff15e059b src/OPERATORS/op-dm-scm.cc --- a/src/OPERATORS/op-dm-scm.cc Mon Mar 09 17:49:13 2009 -0400 +++ b/src/OPERATORS/op-dm-scm.cc Mon Mar 09 17:49:14 2009 -0400 @@ -35,6 +35,8 @@ #include "ov-re-sparse.h" #include "ov-cx-sparse.h" +#include "sparse-xdiv.h" + // diagonal matrix by sparse matrix ops DEFBINOP (mul_dm_scm, diag_matrix, sparse_complex_matrix) @@ -106,6 +108,36 @@ } } +DEFBINOP (ldiv_dm_scm, diag_matrix, sparse_complex_matrix) +{ + CAST_BINOP_ARGS (const octave_diag_matrix&, + const octave_sparse_complex_matrix&); + + MatrixType typ = v2.matrix_type (); + return xleftdiv (v1.diag_matrix_value (), v2.sparse_complex_matrix_value (), + typ); +} + +DEFBINOP (ldiv_cdm_sm, complex_diag_matrix, sparse_matrix) +{ + CAST_BINOP_ARGS (const octave_complex_diag_matrix&, + const octave_sparse_matrix&); + + MatrixType typ = v2.matrix_type (); + return xleftdiv (v1.complex_diag_matrix_value (), v2.sparse_matrix_value (), + typ); +} + +DEFBINOP (ldiv_cdm_scm, complex_diag_matrix, sparse_complex_matrix) +{ + CAST_BINOP_ARGS (const octave_complex_diag_matrix&, + const octave_sparse_complex_matrix&); + + MatrixType typ = v2.matrix_type (); + return xleftdiv (v1.complex_diag_matrix_value (), v2.sparse_complex_matrix_value (), + typ); +} + // sparse matrix by diagonal matrix ops DEFBINOP (mul_scm_dm, sparse_complex_matrix, diag_matrix) @@ -184,6 +216,66 @@ } } +DEFBINOP (div_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) + { + double d = v2.scalar_value (); + + if (d == 0.0) + gripe_divide_by_zero (); + + return octave_value (v1.sparse_complex_matrix_value () / d); + } + else + { + MatrixType typ = v2.matrix_type (); + return xdiv (v1.sparse_complex_matrix_value (), v2.diag_matrix_value (), typ); + } +} + +DEFBINOP (div_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) + { + std::complex d = v2.complex_value (); + + if (d == 0.0) + gripe_divide_by_zero (); + + return octave_value (v1.sparse_matrix_value () / d); + } + else + { + MatrixType typ = v2.matrix_type (); + return xdiv (v1.sparse_matrix_value (), v2.complex_diag_matrix_value (), typ); + } +} + +DEFBINOP (div_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) + { + std::complex d = v2.complex_value (); + + if (d == 0.0) + gripe_divide_by_zero (); + + return octave_value (v1.sparse_complex_matrix_value () / d); + } + else + { + MatrixType typ = v2.matrix_type (); + return xdiv (v1.sparse_complex_matrix_value (), v2.complex_diag_matrix_value (), typ); + } +} + void install_dm_scm_ops (void) { @@ -193,10 +285,19 @@ mul_cdm_sm); INSTALL_BINOP (op_mul, octave_complex_diag_matrix, octave_sparse_complex_matrix, mul_cdm_scm); + INSTALL_BINOP (op_ldiv, octave_diag_matrix, octave_sparse_complex_matrix, ldiv_dm_scm); + INSTALL_BINOP (op_ldiv, octave_complex_diag_matrix, octave_sparse_matrix, ldiv_cdm_sm); + INSTALL_BINOP (op_ldiv, octave_complex_diag_matrix, octave_sparse_complex_matrix, + ldiv_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, mul_sm_cdm); INSTALL_BINOP (op_mul, octave_sparse_complex_matrix, octave_complex_diag_matrix, mul_scm_cdm); + + 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); } diff -r f4f4d65faaa0 -r 42aff15e059b src/OPERATORS/op-dm-sm.cc --- a/src/OPERATORS/op-dm-sm.cc Mon Mar 09 17:49:13 2009 -0400 +++ b/src/OPERATORS/op-dm-sm.cc Mon Mar 09 17:49:14 2009 -0400 @@ -33,6 +33,8 @@ #include "ov-re-diag.h" #include "ov-re-sparse.h" +#include "sparse-xdiv.h" + // diagonal matrix by sparse matrix ops DEFBINOP (mul_dm_sm, diag_matrix, sparse_matrix) @@ -58,6 +60,14 @@ } } +DEFBINOP (ldiv_dm_sm, diag_matrix, sparse_matrix) +{ + CAST_BINOP_ARGS (const octave_diag_matrix&, const octave_sparse_matrix&); + + MatrixType typ = v2.matrix_type (); + return xleftdiv (v1.diag_matrix_value (), v2.sparse_matrix_value (), typ); +} + // sparse matrix by diagonal matrix ops DEFBINOP (mul_sm_dm, sparse_matrix, diag_matrix) @@ -83,12 +93,36 @@ } } +DEFBINOP (div_sm_dm, sparse_matrix, diag_matrix) +{ + CAST_BINOP_ARGS (const octave_sparse_matrix&, const octave_diag_matrix&); + + if (v2.rows() == 1 && v2.columns() == 1) + { + double d = v2.scalar_value (); + + if (d == 0.0) + gripe_divide_by_zero (); + + return octave_value (v1.sparse_matrix_value () / d); + } + else + { + MatrixType typ = v2.matrix_type (); + return xdiv (v1.sparse_matrix_value (), v2.diag_matrix_value (), typ); + } +} + void install_dm_sm_ops (void) { INSTALL_BINOP (op_mul, octave_diag_matrix, octave_sparse_matrix, mul_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_div, octave_sparse_matrix, octave_diag_matrix, div_sm_dm); } diff -r f4f4d65faaa0 -r 42aff15e059b src/sparse-xdiv.cc --- a/src/sparse-xdiv.cc Mon Mar 09 17:49:13 2009 -0400 +++ b/src/sparse-xdiv.cc Mon Mar 09 17:49:14 2009 -0400 @@ -33,7 +33,9 @@ #include "error.h" #include "dSparse.h" +#include "dDiagMatrix.h" #include "CSparse.h" +#include "CDiagMatrix.h" #include "oct-spparms.h" #include "sparse-xdiv.h" @@ -74,6 +76,10 @@ INSTANTIATE_MX_LEFTDIV_CONFORM (SparseMatrix, ComplexMatrix); INSTANTIATE_MX_LEFTDIV_CONFORM (SparseComplexMatrix, Matrix); INSTANTIATE_MX_LEFTDIV_CONFORM (SparseComplexMatrix, ComplexMatrix); +INSTANTIATE_MX_LEFTDIV_CONFORM (DiagMatrix, SparseMatrix); +INSTANTIATE_MX_LEFTDIV_CONFORM (DiagMatrix, SparseComplexMatrix); +INSTANTIATE_MX_LEFTDIV_CONFORM (ComplexDiagMatrix, SparseMatrix); +INSTANTIATE_MX_LEFTDIV_CONFORM (ComplexDiagMatrix, SparseComplexMatrix); template bool @@ -105,15 +111,23 @@ INSTANTIATE_MX_DIV_CONFORM (Matrix, SparseComplexMatrix); INSTANTIATE_MX_DIV_CONFORM (ComplexMatrix, SparseMatrix); INSTANTIATE_MX_DIV_CONFORM (ComplexMatrix, SparseComplexMatrix); +INSTANTIATE_MX_DIV_CONFORM (SparseMatrix, DiagMatrix); +INSTANTIATE_MX_DIV_CONFORM (SparseMatrix, ComplexDiagMatrix); +INSTANTIATE_MX_DIV_CONFORM (SparseComplexMatrix, DiagMatrix); +INSTANTIATE_MX_DIV_CONFORM (SparseComplexMatrix, ComplexDiagMatrix); -// Right division functions. +// Right division functions. X / Y = X * inv(Y) = (inv (Y') * X')' // -// op2 / op1: m cm sm scm +// Y / X: m cm sm scm // +-- +---+----+----+----+ // sparse matrix | 1 | 3 | 5 | 7 | // +---+----+----+----+ // sparse complex_matrix | 2 | 4 | 6 | 8 | // +---+----+----+----+ +// diagonal matrix | 9 | 11 | +// +----+----+ +// complex diag. matrix | 10 | 12 | +// +----+----+ // -*- 1 -*- Matrix @@ -275,6 +289,74 @@ return result.hermitian (); } +template +RT do_rightdiv_sm_dm (const SM& a, const DM& d) +{ + const octave_idx_type d_nr = d.rows (); + + const octave_idx_type a_nr = a.rows (); + const octave_idx_type a_nc = a.cols (); + + using std::min; + const octave_idx_type nc = min (d_nr, a_nc); + + if ( ! mx_div_conform (a, d)) + return RT (); + + const octave_idx_type nz = a.nnz (); + RT r (a_nr, nc, nz); + + const typename DM::element_type zero = typename DM::element_type (); + + octave_idx_type k_result = 0; + for (octave_idx_type j = 0; j < nc; ++j) + { + OCTAVE_QUIT; + const typename DM::element_type s = d.dgelem (j); + const octave_idx_type colend = a.cidx (j+1); + r.xcidx (j) = k_result; + if (s != zero) + for (octave_idx_type k = a.cidx (j); k < colend; ++k) + { + r.xdata (k_result) = a.data (k) / s; + r.xridx (k_result) = a.ridx (k); + ++k_result; + } + } + r.xcidx (nc) = k_result; + + r.maybe_compress (true); + return r; +} + +// -*- 9 -*- +SparseMatrix +xdiv (const SparseMatrix& a, const DiagMatrix& b, MatrixType &) +{ + return do_rightdiv_sm_dm (a, b); +} + +// -*- 10 -*- +SparseComplexMatrix +xdiv (const SparseMatrix& a, const ComplexDiagMatrix& b, MatrixType &) +{ + return do_rightdiv_sm_dm (a, b); +} + +// -*- 11 -*- +SparseComplexMatrix +xdiv (const SparseComplexMatrix& a, const DiagMatrix& b, MatrixType &) +{ + return do_rightdiv_sm_dm (a, b); +} + +// -*- 12 -*- +SparseComplexMatrix +xdiv (const SparseComplexMatrix& a, const ComplexDiagMatrix& b, MatrixType &) +{ + return do_rightdiv_sm_dm (a, b); +} + // Funny element by element division operations. // // op2 \ op1: s cs @@ -363,18 +445,18 @@ return result; } -// Left division functions. +// Left division functions. X \ Y = inv(X) * Y // -// op2 \ op1: m cm +// Y \ X : sm scm dm dcm // +-- +---+----+ // matrix | 1 | 5 | // +---+----+ // complex_matrix | 2 | 6 | -// +---+----+ -// sparse matrix | 3 | 7 | -// +---+----+ -// sparse complex_matrix | 4 | 8 | -// +---+----+ +// +---+----+----+----+ +// sparse matrix | 3 | 7 | 9 | 11 | +// +---+----+----+----+ +// sparse complex_matrix | 4 | 8 | 10 | 12 | +// +---+----+----+----+ // -*- 1 -*- Matrix @@ -473,6 +555,80 @@ return a.solve (typ, b, info, rcond, solve_singularity_warning); } +template +RT do_leftdiv_dm_sm (const DM& d, const SM& a) +{ + const octave_idx_type a_nr = a.rows (); + const octave_idx_type a_nc = a.cols (); + + const octave_idx_type d_nc = d.cols (); + + using std::min; + const octave_idx_type nr = min (d_nc, a_nr); + + if ( ! mx_leftdiv_conform (d, a)) + return RT (); + + const octave_idx_type nz = a.nnz (); + RT r (nr, a_nc, nz); + + const typename DM::element_type zero = typename DM::element_type (); + + octave_idx_type k_result = 0; + for (octave_idx_type j = 0; j < a_nc; ++j) + { + OCTAVE_QUIT; + const octave_idx_type colend = a.cidx (j+1); + r.xcidx (j) = k_result; + for (octave_idx_type k = a.cidx (j); k < colend; ++k) + { + const octave_idx_type i = a.ridx (k); + if (i < nr) + { + const typename DM::element_type s = d.dgelem (i); + if (s != zero) + { + r.xdata (k_result) = a.data (k) / s; + r.xridx (k_result) = i; + ++k_result; + } + } + } + } + r.xcidx (a_nc) = k_result; + + r.maybe_compress (true); + return r; +} + +// -*- 9 -*- +SparseMatrix +xleftdiv (const DiagMatrix& d, const SparseMatrix& a, MatrixType&) +{ + return do_leftdiv_dm_sm (d, a); +} + +// -*- 10 -*- +SparseComplexMatrix +xleftdiv (const DiagMatrix& d, const SparseComplexMatrix& a, MatrixType&) +{ + return do_leftdiv_dm_sm (d, a); +} + +// -*- 11 -*- +SparseComplexMatrix +xleftdiv (const ComplexDiagMatrix& d, const SparseMatrix& a, MatrixType&) +{ + return do_leftdiv_dm_sm (d, a); +} + +// -*- 12 -*- +SparseComplexMatrix +xleftdiv (const ComplexDiagMatrix& d, const SparseComplexMatrix& a, MatrixType&) +{ + return do_leftdiv_dm_sm (d, a); +} + /* ;;; Local Variables: *** ;;; mode: C++ *** diff -r f4f4d65faaa0 -r 42aff15e059b src/sparse-xdiv.h --- a/src/sparse-xdiv.h Mon Mar 09 17:49:13 2009 -0400 +++ b/src/sparse-xdiv.h Mon Mar 09 17:49:14 2009 -0400 @@ -27,6 +27,8 @@ #include "oct-cmplx.h" #include "MatrixType.h" +class DiagMatrix; +class ComplexDiagMatrix; class SparseMatrix; class SparseComplexMatrix; @@ -47,6 +49,15 @@ extern SparseComplexMatrix xdiv (const SparseComplexMatrix& a, const SparseComplexMatrix& b, MatrixType &typ); +extern SparseMatrix xdiv (const SparseMatrix& a, + const DiagMatrix& b, MatrixType &typ); +extern SparseComplexMatrix xdiv (const SparseMatrix& a, + const ComplexDiagMatrix& b, MatrixType &typ); +extern SparseComplexMatrix xdiv (const SparseComplexMatrix& a, + const DiagMatrix& b, MatrixType &typ); +extern SparseComplexMatrix xdiv (const SparseComplexMatrix& a, + const ComplexDiagMatrix& b, MatrixType &typ); + extern Matrix x_el_div (double a, const SparseMatrix& b); extern ComplexMatrix x_el_div (double a, const SparseComplexMatrix& b); extern ComplexMatrix x_el_div (const Complex a, const SparseMatrix& b); @@ -71,6 +82,14 @@ extern SparseComplexMatrix xleftdiv (const SparseComplexMatrix& a, const SparseComplexMatrix& b, MatrixType &typ); +extern SparseMatrix xleftdiv (const DiagMatrix&, const SparseMatrix&, MatrixType&); +extern SparseComplexMatrix xleftdiv (const ComplexDiagMatrix&, const SparseMatrix&, + MatrixType&); +extern SparseComplexMatrix xleftdiv (const DiagMatrix&, const SparseComplexMatrix&, + MatrixType&); +extern SparseComplexMatrix xleftdiv (const ComplexDiagMatrix&, const SparseComplexMatrix&, + MatrixType&); + #endif /* diff -r f4f4d65faaa0 -r 42aff15e059b test/ChangeLog --- a/test/ChangeLog Mon Mar 09 17:49:13 2009 -0400 +++ b/test/ChangeLog Mon Mar 09 17:49:14 2009 -0400 @@ -8,6 +8,10 @@ 2009-03-08 Jason Riedy + * test_diag_perm.m: Add tests for inverse scaling and sparse structure. + +2009-03-08 Jason Riedy + * test_diag_perm.m: Add tests for preserving sparse structure when scaling. 2009-02-25 John W. Eaton diff -r f4f4d65faaa0 -r 42aff15e059b test/test_diag_perm.m --- a/test/test_diag_perm.m Mon Mar 09 17:49:13 2009 -0400 +++ b/test/test_diag_perm.m Mon Mar 09 17:49:14 2009 -0400 @@ -157,3 +157,38 @@ %! A(4, 1) = Inf; %! assert (Dr * A * Dc, A .* kron (dr, dc), eps); +## sparse inverse row scaling with a zero factor +%!test +%! n = 8; +%! A = sprand (n, n, .5); +%! scalefact = rand (n, 1); +%! Dr = diag (scalefact); +%! scalefact(n-1) = Inf; +%! Dr(n-1, n-1) = 0; +%! assert (full (Dr \ A), full (A) ./ repmat (scalefact, 1, n)); + +## narrow sparse inverse row scaling +%!test +%! n = 8; +%! A = sprand (n, n, .5); +%! scalefact = rand (n-2, 1); +%! Dr = diag (scalefact, n, n-2); +%! assert (full (Dr \ A), Dr \ full(A)) + +## sparse inverse column scaling with a zero factor +%!test +%! n = 11; +%! A = sprand (n, n, .5); +%! scalefact = rand (1, n); +%! Dc = diag (scalefact); +%! scalefact(n-1) = Inf; +%! Dc(n-1, n-1) = 0; +%! assert (full (A / Dc), full(A) / Dc) + +## short sparse inverse column scaling +%!test +%! n = 7; +%! A = sprand (n, n, .5); +%! scalefact = rand (1, n-2) + I () * rand(1, n-2); +%! Dc = diag (scalefact, n-2, n); +%! assert (full (A / Dc), full(A) / Dc)