Mercurial > hg > octave-terminal
changeset 8968:91d53dc37f79
Add perm * sparse, perm \ sparse, sparse * perm, and sparse / perm operations.
Nothing terribly fancy in any of this. There probably is some
mechanism for using the permutation vectors and some assign or index
method in the sparse classes, but I've never understood all the
intricacies. I'm opting for a simple implementation at the cost of
possibly duplicating some functionality.
author | Jason Riedy <jason@acm.org> |
---|---|
date | Tue, 10 Mar 2009 21:54:44 -0400 |
parents | 5bbbf482909a |
children | 3ecbc236e2e0 |
files | liboctave/CSparse.cc liboctave/CSparse.h liboctave/ChangeLog liboctave/Sparse-perm-op-defs.h liboctave/dSparse.cc liboctave/dSparse.h src/ChangeLog src/Makefile.in src/OPERATORS/op-pm-scm.cc src/OPERATORS/op-pm-sm.cc test/ChangeLog test/test_diag_perm.m |
diffstat | 12 files changed, 499 insertions(+), 1 deletions(-) [+] |
line wrap: on
line diff
--- a/liboctave/CSparse.cc +++ b/liboctave/CSparse.cc @@ -52,6 +52,8 @@ #include "Sparse-diag-op-defs.h" +#include "Sparse-perm-op-defs.h" + // Define whether to use a basic QR solver or one that uses a Dulmange // Mendelsohn factorization to seperate the problem into under-determined, // well-determined and over-determined parts and solves them seperately @@ -7684,6 +7686,20 @@ return do_sub_sm_dm<SparseComplexMatrix> (a, d); } +// perm * sparse and sparse * perm + +SparseComplexMatrix +operator * (const PermMatrix& p, const SparseComplexMatrix& a) +{ + return octinternal_do_mul_pm_sm (p, a); +} + +SparseComplexMatrix +operator * (const SparseComplexMatrix& a, const PermMatrix& p) +{ + return octinternal_do_mul_sm_pm (a, p); +} + // FIXME -- it would be nice to share code among the min/max // functions below.
--- a/liboctave/CSparse.h +++ b/liboctave/CSparse.h @@ -37,6 +37,7 @@ #include "Sparse-op-defs.h" #include "MatrixType.h" +class PermMatrix; class DiagMatrix; class ComplexDiagMatrix; class SparseMatrix; @@ -498,6 +499,11 @@ extern OCTAVE_API SparseComplexMatrix operator - (const SparseComplexMatrix&, const DiagMatrix&); extern OCTAVE_API SparseComplexMatrix operator - (const SparseComplexMatrix&, const ComplexDiagMatrix&); +extern OCTAVE_API SparseComplexMatrix operator * (const PermMatrix&, + const SparseComplexMatrix&); +extern OCTAVE_API SparseComplexMatrix operator * (const SparseComplexMatrix&, + const PermMatrix&); + extern OCTAVE_API SparseComplexMatrix min (const Complex& c, const SparseComplexMatrix& m); extern OCTAVE_API SparseComplexMatrix min (const SparseComplexMatrix& m,
--- a/liboctave/ChangeLog +++ b/liboctave/ChangeLog @@ -89,6 +89,32 @@ 2009-03-10 Jason Riedy <jason@acm.org> + * Sparse-perm-op-defs.h (octinternal_do_mul_colpm_sm): New + template function. Logic for the column permutation * sparse + matrix operator. + (octinternal_do_mul_pm_sm): New template function. Logic for the + permutation matrix * sparse matrix operator. Note that there is + no special row perm * sparse routine; the permutation is inverted + and the col perm routine is called. + (octinternal_do_mul_sm_rowpm): New template function. Logic for + the sparse matrix * row permutation operator. + (octinternal_do_mul_sm_colpm): New template function. Logic for + the sparse matrix * column permutation operator. + (octinternal_do_mul_sm_pm): New template function. Logic for the + sparse matrix * permutation matrix operator. + + * dSparse.h (operator *): Declare sparse * permutation and + permutation * sparse. + * dSparse.cc (operator *): Define sparse * permutation and + permutation * sparse. + + * CSparse.h (operator *): Declare sparse * permutation and + permutation * sparse. + * CSparse.cc (operator *): Define sparse * permutation and + permutation * sparse. + +2009-03-10 Jason Riedy <jason@acm.org> + * sparse-base-lu.cc (Pc_vec): The column permutation should be Ufact.cols ()-long, not Lfact.rows ()-long.
new file mode 100644 --- /dev/null +++ b/liboctave/Sparse-perm-op-defs.h @@ -0,0 +1,164 @@ +/* -*- C++ -*- + +Copyright (C) 2009 Jason Riedy + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 3 of the License, or (at your +option) any later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, see +<http://www.gnu.org/licenses/>. + +*/ + +#if !defined (octave_sparse_perm_op_defs_h) +#define octave_sparse_perm_op_defs_h 1 + +// Matrix multiplication + +template <typename SM> +SM octinternal_do_mul_colpm_sm (const octave_idx_type *pcol, const SM& a) +// Relabel the rows according to pcol. +{ + const octave_idx_type nr = a.rows (); + const octave_idx_type nc = a.cols (); + const octave_idx_type nent = a.nnz (); + SM r (nr, nc, nent); + + for (octave_idx_type k = 0; k < nent; ++k) + { + OCTAVE_QUIT; + r.xridx (k) = pcol[a.ridx (k)]; + r.xdata (k) = a.data (k); + } + for (octave_idx_type j = 0; j <= nc; ++j) + r.xcidx (j) = a.cidx (j); + + r.maybe_compress (false); + return r; +} + +template <typename SM> +SM octinternal_do_mul_pm_sm (const PermMatrix& p, const SM& a) +{ + const octave_idx_type nr = a.rows (); + if (p.cols () != nr) + { + gripe_nonconformant ("operator *", p.rows (), p.cols (), a.rows (), a.cols ()); + return SM (); + } + + if (p.is_row_perm ()) + { + // Form the column permutation and then call the colpm_sm routine. + const octave_idx_type *prow = p.pvec ().data (); + OCTAVE_LOCAL_BUFFER(octave_idx_type, pcol, nr); + for (octave_idx_type i = 0; i < nr; ++i) + pcol[prow[i]] = i; + return octinternal_do_mul_colpm_sm (pcol, a); + } + else + return octinternal_do_mul_colpm_sm (p.pvec ().data (), a); +} + +template <typename SM> +SM octinternal_do_mul_sm_rowpm (const SM& a, const octave_idx_type *prow) +// For a row permutation, iterate across the source a and stuff the +// results into the correct destination column in r. +{ + const octave_idx_type nr = a.rows (); + const octave_idx_type nc = a.cols (); + const octave_idx_type nent = a.nnz (); + SM r (nr, nc, nent); + + for (octave_idx_type j_src = 0; j_src < nc; ++j_src) + r.xcidx (prow[j_src]) = a.cidx (j_src+1) - a.cidx (j_src); + octave_idx_type k = 0; + for (octave_idx_type j = 0; j < nc; ++j) + { + const octave_idx_type tmp = r.xcidx (j); + r.xcidx (j) = k; + k += tmp; + } + r.xcidx (nc) = nent; + + octave_idx_type k_src = 0; + for (octave_idx_type j_src = 0; j_src < nc; ++j_src) + { + OCTAVE_QUIT; + const octave_idx_type j = prow[j_src]; + const octave_idx_type kend_src = a.cidx (j_src + 1); + for (k = r.xcidx (j); k_src < kend_src; ++k, ++k_src) + { + r.xridx (k) = a.ridx (k_src); + r.xdata (k) = a.data (k_src); + } + } + assert (k_src == nent); + + r.maybe_compress (false); + return r; +} + +template <typename SM> +SM octinternal_do_mul_sm_colpm (const SM& a, const octave_idx_type *pcol) +// For a column permutation, iterate across the destination r and pull +// data from the correct column of a. +{ + const octave_idx_type nr = a.rows (); + const octave_idx_type nc = a.cols (); + const octave_idx_type nent = a.nnz (); + SM r (nr, nc, nent); + + for (octave_idx_type j = 0; j < nc; ++j) + { + const octave_idx_type j_src = pcol[j]; + r.xcidx (j+1) = r.xcidx (j) + (a.cidx (j_src+1) - a.cidx (j_src)); + } + assert (r.xcidx (nc) == nent); + + octave_idx_type k = 0; + for (octave_idx_type j = 0; j < nc; ++j) + { + OCTAVE_QUIT; + const octave_idx_type j_src = pcol[j]; + octave_idx_type k_src; + const octave_idx_type kend_src = a.cidx (j_src + 1); + for (k_src = a.cidx (j_src); k_src < kend_src; ++k_src, ++k) + { + r.xridx (k) = a.ridx (k_src); + r.xdata (k) = a.data (k_src); + } + } + assert (k == nent); + + r.maybe_compress (false); + return r; +} + +template <typename SM> +SM octinternal_do_mul_sm_pm (const SM& a, const PermMatrix& p) +{ + const octave_idx_type nc = a.cols (); + if (p.rows () != nc) + { + gripe_nonconformant ("operator *", a.rows (), a.cols (), p.rows (), p.cols ()); + return SM (); + } + + if (p.is_row_perm ()) + return octinternal_do_mul_sm_rowpm (a, p.pvec ().data ()); + else + return octinternal_do_mul_sm_colpm (a, p.pvec ().data ()); +} + +#endif // octave_sparse_perm_op_defs_h
--- a/liboctave/dSparse.cc +++ b/liboctave/dSparse.cc @@ -53,6 +53,8 @@ #include "Sparse-diag-op-defs.h" +#include "Sparse-perm-op-defs.h" + // Define whether to use a basic QR solver or one that uses a Dulmange // Mendelsohn factorization to seperate the problem into under-determined, // well-determined and over-determined parts and solves them seperately @@ -7740,6 +7742,20 @@ return do_sub_sm_dm<SparseMatrix> (a, d); } +// perm * sparse and sparse * perm + +SparseMatrix +operator * (const PermMatrix& p, const SparseMatrix& a) +{ + return octinternal_do_mul_pm_sm (p, a); +} + +SparseMatrix +operator * (const SparseMatrix& a, const PermMatrix& p) +{ + return octinternal_do_mul_sm_pm (a, p); +} + // FIXME -- it would be nice to share code among the min/max // functions below.
--- a/liboctave/dSparse.h +++ b/liboctave/dSparse.h @@ -36,6 +36,7 @@ #include "Sparse-op-defs.h" #include "MatrixType.h" +class PermMatrix; class DiagMatrix; class SparseComplexMatrix; class SparseBoolMatrix; @@ -459,6 +460,9 @@ extern OCTAVE_API SparseMatrix operator - (const DiagMatrix&, const SparseMatrix&); extern OCTAVE_API SparseMatrix operator - (const SparseMatrix&, const DiagMatrix&); +extern OCTAVE_API SparseMatrix operator * (const PermMatrix&, const SparseMatrix&); +extern OCTAVE_API SparseMatrix operator * (const SparseMatrix&, const PermMatrix&); + 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);
--- a/src/ChangeLog +++ b/src/ChangeLog @@ -47,6 +47,25 @@ 2009-03-10 Jason Riedy <jason@acm.org> + * OPERATORS/op-pm-sm.cc (mul_pm_sm): New Octave binding for + perm * sparse. + (ldiv_pm_sm): New Octave binding for perm \ sparse. + (mul_sm_pm): New Octave binding for sparse * perm. + (div_sm_pm): New Octave binding for sparse / perm. + (install_pm_sm_ops): Install the above bindings. + + * OPERATORS/op-pm-scm.cc (mul_pm_scm): New Octave binding for + perm * sparse complex. + (ldiv_pm_scm): New Octave binding for perm \ sparse complex. + (mul_scm_pm): New Octave binding for sparse complex * perm. + (div_scm_pm): New Octave binding for sparse complex / perm. + (install_pm_scm_ops): Install the above bindings. + + * Makefile.in (PERM_OP_XSRC): Add op-pm-sm.cc and op-pm-scm.cc for + operations between permutations and sparse matrices. + +2009-03-10 Jason Riedy <jason@acm.org> + * DLD-FUNCTIONS/find.cc (find_nonzero_elem_idx): New override for find on PermMatrix. (find): Add a branch testing arg.is_perm_matrix () and calling the
--- a/src/Makefile.in +++ b/src/Makefile.in @@ -174,7 +174,8 @@ op-fdm-fm.cc op-fdm-fs.cc op-fm-fcdm.cc op-fm-fdm.cc PERM_OP_XSRC := op-cm-pm.cc op-fcm-pm.cc op-fm-pm.cc op-pm-fcm.cc \ - op-pm-fm.cc op-m-pm.cc op-pm-cm.cc op-pm-m.cc op-pm-pm.cc + op-pm-fm.cc op-m-pm.cc op-pm-cm.cc op-pm-m.cc op-pm-pm.cc \ + op-pm-sm.cc op-pm-scm.cc OP_XSRC := op-b-b.cc op-b-bm.cc op-bm-b.cc op-bm-bm.cc op-cell.cc \ op-chm.cc op-class.cc op-list.cc op-range.cc op-str-m.cc \
new file mode 100644 --- /dev/null +++ b/src/OPERATORS/op-pm-scm.cc @@ -0,0 +1,97 @@ +/* + +Copyright (C) 2009 Jason Riedy + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 3 of the License, or (at your +option) any later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, see +<http://www.gnu.org/licenses/>. + +*/ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include "gripes.h" +#include "oct-obj.h" +#include "ov.h" +#include "ov-typeinfo.h" +#include "ops.h" + +#include "ov-perm.h" +#include "ov-cx-sparse.h" + +// permutation matrix by sparse matrix ops + +DEFBINOP (mul_pm_scm, perm_matrix, sparse_complex_matrix) +{ + CAST_BINOP_ARGS (const octave_perm_matrix&, const octave_sparse_complex_matrix&); + + if (v2.rows() == 1 && v2.columns() == 1) + { + std::complex<double> d = v2.complex_value (); + + return octave_value (v1.sparse_matrix_value () * d); + } + else if (v1.rows() == 1 && v1.columns() == 1) + return octave_value (v2.sparse_complex_matrix_value ()); + else + return v1.perm_matrix_value () * v2.sparse_complex_matrix_value (); +} + +DEFBINOP (ldiv_pm_scm, perm_matrix, sparse_complex_matrix) +{ + CAST_BINOP_ARGS (const octave_perm_matrix&, const octave_sparse_complex_matrix&); + + return v1.perm_matrix_value ().inverse () * v2.sparse_complex_matrix_value (); +} + +// sparse matrix by diagonal matrix ops + +DEFBINOP (mul_scm_pm, sparse_complex_matrix, perm_matrix) +{ + CAST_BINOP_ARGS (const octave_sparse_complex_matrix&, const octave_perm_matrix&); + + if (v1.rows() == 1 && v1.columns() == 1) + { + std::complex<double> d = v1.scalar_value (); + + return octave_value (d * v2.sparse_matrix_value ()); + } + else if (v2.rows() == 1 && v2.columns() == 1) + return octave_value (v1.sparse_complex_matrix_value ()); + else + return v1.sparse_complex_matrix_value () * v2.perm_matrix_value (); +} + +DEFBINOP (div_scm_pm, sparse_complex_matrix, perm_matrix) +{ + CAST_BINOP_ARGS (const octave_sparse_complex_matrix&, const octave_perm_matrix&); + + return v1.sparse_complex_matrix_value () * v2.perm_matrix_value ().inverse (); +} + +void +install_pm_scm_ops (void) +{ + INSTALL_BINOP (op_mul, octave_perm_matrix, octave_sparse_complex_matrix, + mul_pm_scm); + INSTALL_BINOP (op_ldiv, octave_perm_matrix, octave_sparse_complex_matrix, + ldiv_pm_scm); + INSTALL_BINOP (op_mul, octave_sparse_complex_matrix, octave_perm_matrix, + mul_scm_pm); + INSTALL_BINOP (op_div, octave_sparse_complex_matrix, octave_perm_matrix, + div_scm_pm); +}
new file mode 100644 --- /dev/null +++ b/src/OPERATORS/op-pm-sm.cc @@ -0,0 +1,97 @@ +/* + +Copyright (C) 2009 Jason Riedy + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 3 of the License, or (at your +option) any later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, see +<http://www.gnu.org/licenses/>. + +*/ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include "gripes.h" +#include "oct-obj.h" +#include "ov.h" +#include "ov-typeinfo.h" +#include "ops.h" + +#include "ov-perm.h" +#include "ov-re-sparse.h" + +// permutation matrix by sparse matrix ops + +DEFBINOP (mul_pm_sm, perm_matrix, sparse_matrix) +{ + CAST_BINOP_ARGS (const octave_perm_matrix&, const octave_sparse_matrix&); + + if (v2.rows() == 1 && v2.columns() == 1) + { + double d = v2.scalar_value (); + + return octave_value (v1.sparse_matrix_value () * d); + } + else if (v1.rows() == 1 && v1.columns() == 1) + return octave_value (v2.sparse_matrix_value ()); + else + return v1.perm_matrix_value () * v2.sparse_matrix_value (); +} + +DEFBINOP (ldiv_pm_sm, perm_matrix, sparse_matrix) +{ + CAST_BINOP_ARGS (const octave_perm_matrix&, const octave_sparse_matrix&); + + return v1.perm_matrix_value ().inverse () * v2.sparse_matrix_value (); +} + +// sparse matrix by diagonal matrix ops + +DEFBINOP (mul_sm_pm, sparse_matrix, perm_matrix) +{ + CAST_BINOP_ARGS (const octave_sparse_matrix&, const octave_perm_matrix&); + + if (v1.rows() == 1 && v1.columns() == 1) + { + double d = v1.scalar_value (); + + return octave_value (d * v2.sparse_matrix_value ()); + } + else if (v2.rows() == 1 && v2.columns() == 1) + return octave_value (v1.sparse_matrix_value ()); + else + return v1.sparse_matrix_value () * v2.perm_matrix_value (); +} + +DEFBINOP (div_sm_pm, sparse_matrix, perm_matrix) +{ + CAST_BINOP_ARGS (const octave_sparse_matrix&, const octave_perm_matrix&); + + return v1.sparse_matrix_value () * v2.perm_matrix_value ().inverse (); +} + +void +install_pm_sm_ops (void) +{ + INSTALL_BINOP (op_mul, octave_perm_matrix, octave_sparse_matrix, + mul_pm_sm); + INSTALL_BINOP (op_ldiv, octave_perm_matrix, octave_sparse_matrix, + ldiv_pm_sm); + INSTALL_BINOP (op_mul, octave_sparse_matrix, octave_perm_matrix, + mul_sm_pm); + INSTALL_BINOP (op_div, octave_sparse_matrix, octave_perm_matrix, + div_sm_pm); +}
--- a/test/ChangeLog +++ b/test/ChangeLog @@ -1,3 +1,9 @@ +2009-03-10 Jason Riedy <jason@acm.org> + + * test_diag_perm.m: Add tests for permuting sparse matrices and + for the correct types from interactions between + pseudo-scalars (1x1 matrices). + 2009-03-10 Jason Riedy <jason@acm.org> * build_sparse_tests.sh: Add LU tests to the rectangular tests.
--- a/test/test_diag_perm.m +++ b/test/test_diag_perm.m @@ -57,6 +57,52 @@ %! assert (typeinfo (rand (n) * Pc), "matrix"); %! assert (typeinfo (Pr * rand (n)), "matrix"); +## preserve sparse matrix structure +%!test +%! n = 7; +%! Pc = eye (n) (:, randperm (n)); +%! Ac = sprand (n-3, n, .5) + I () * sprand (n-3, n, .5); +%! Pr = eye (n) (randperm (n), :); +%! Ar = sprand (n, n+2, .5); +%! assert (typeinfo (Ac * Pc), "sparse complex matrix"); +%! assert (full (Ac * Pc), full (Ac) * Pc); +%! assert (full (Ac / Pc), full (Ac) / Pc); +%! assert (typeinfo (Pr * Ar), "sparse matrix"); +%! assert (full (Pr * Ar), Pr * full (Ar)); +%! assert (full (Pr \ Ar), Pr \ full (Ar)); + +## structure rules for 1x1 dense / scalar and 1x1 perm +%!test +%! n = 7; +%! P1 = eye (1) (:, [1]); +%! A1 = 1; +%! P = eye (n) (:, randperm (n)); +%! A = rand (n-3, n, .5); +%! assert (typeinfo (A * P1), "matrix"); +%! assert (full (A * P1), full (A) * P1); +%! assert (typeinfo (P1 * A), "matrix"); +%! assert (full (P1 * A), P1 * full (A)); +%! assert (typeinfo (A1 * P), "matrix"); +%! assert (full (A1 * P), full (A1) * P); +%! assert (typeinfo (P * A1), "matrix"); +%! assert (full (P * A1), P * full (A1)); + +## structure rules for 1x1 sparse and 1x1 perm +%!test +%! n = 7; +%! P1 = eye (1) (:, [1]); +%! A1 = sparse (1, 1, 2); +%! P = eye (n) (:, randperm (n)); +%! A = sprand (n-3, n, .5); +%! assert (typeinfo (A * P1), "sparse matrix"); +%! assert (full (A * P1), full (A) * P1); +%! assert (typeinfo (P1 * A), "sparse matrix"); +%! assert (full (P1 * A), P1 * full (A)); +%! assert (typeinfo (A1 * P), "sparse matrix"); +%! assert (full (A1 * P), full (A1) * P); +%! assert (typeinfo (P * A1), "sparse matrix"); +%! assert (full (P * A1), P * full (A1)); + ## permuting a matrix with exceptional values does not introduce new ones. %!test %! n = 5;