[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[PATCH 2 of 4] Implement sparse * diagonal and diagonal * sparse operati
|
From: |
Jason Riedy |
|
Subject: |
[PATCH 2 of 4] Implement sparse * diagonal and diagonal * sparse operations, double-prec only |
|
Date: |
Mon, 09 Mar 2009 17:52:20 -0400 |
# HG changeset patch
# User Jason Riedy <address@hidden>
# Date 1236635353 14400
# Node ID 6f9fc3415a747112998fbaa662bbbccef5c11631
# Parent 8ec69bd63c7360d9b186677968a5f93038c8c899
Implement sparse * diagonal and diagonal * sparse operations, double-prec only.
>From da16376d5519cf82cb3384330739784fb00fe57c Mon Sep 17 00:00:00 2001
Date: Sun, 8 Mar 2009 16:28:18 -0400
These preserve sparsity, so eye(5) * sprand (5, 5, .2) is *sparse*
and not dense. This may affect people who use multiplication by
eye() rather than full().
The liboctave routines do *not* check if arguments are scalars in
disguise. There is a type problem with checking at that level. I
suspect we want diag * "sparse scalar" to stay diagonal, but we have
to return a sparse matrix at the liboctave. Rather than worrying
about that in liboctave, we cope with it when binding to Octave and
return the correct higher-level type.
The implementation is in Sparse-diag-op-defs.h rather than
Sparse-op-defs.h to limit recompilation. And the implementations
are templates rather than macros to produce better compiler errors
and debugging information.
Signed-off-by: Jason Riedy <address@hidden>
---
liboctave/CSparse.cc | 102 ++++++++++
liboctave/CSparse.h | 23 +++
liboctave/ChangeLog | 15 ++
liboctave/Sparse-diag-op-defs.h | 117 +++++++++++
liboctave/dSparse.cc | 30 +++
liboctave/dSparse.h | 7 +-
src/ChangeLog | 9 +
src/Makefile.in | 2 +-
src/OPERATORS/op-dm-scm.cc | 403 +++++++++++++++++++++++++++++++++++++++
src/OPERATORS/op-dm-sm.cc | 118 ++++++++++++
test/ChangeLog | 4 +
test/test_diag_perm.m | 9 +-
12 files changed, 836 insertions(+), 3 deletions(-)
create mode 100644 liboctave/Sparse-diag-op-defs.h
create mode 100644 src/OPERATORS/op-dm-scm.cc
create mode 100644 src/OPERATORS/op-dm-sm.cc
diff --git a/liboctave/CSparse.cc b/liboctave/CSparse.cc
--- a/liboctave/CSparse.cc
+++ b/liboctave/CSparse.cc
@@ -37,6 +37,8 @@
#include "dRowVector.h"
#include "oct-locbuf.h"
+#include "dDiagMatrix.h"
+#include "CDiagMatrix.h"
#include "CSparse.h"
#include "boolSparse.h"
#include "dSparse.h"
@@ -48,6 +50,8 @@
#include "SparseCmplxCHOL.h"
#include "SparseCmplxQR.h"
+#include "Sparse-diag-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
@@ -7584,6 +7588,104 @@
SPARSE_FULL_TRANS_MUL (ComplexMatrix, Complex, Complex (0.,0.), conj);
}
+// diag * sparse and sparse * diag
+SparseComplexMatrix
+operator * (const DiagMatrix& d, const SparseComplexMatrix& a)
+{
+ return octave_impl::do_mul_dm_sm<SparseComplexMatrix> (d, a);
+}
+SparseComplexMatrix
+trans_mul (const DiagMatrix& d, const SparseComplexMatrix& a)
+{
+ return operator * (d, a);
+}
+SparseComplexMatrix
+herm_mul (const DiagMatrix& d, const SparseComplexMatrix& a)
+{
+ return operator * (d, a);
+}
+SparseComplexMatrix
+operator * (const SparseComplexMatrix& a, const DiagMatrix& d)
+{
+ return octave_impl::do_mul_sm_dm<SparseComplexMatrix> (a, d);
+}
+SparseComplexMatrix
+mul_trans (const SparseComplexMatrix& a, const DiagMatrix& d)
+{
+ return operator * (a, d);
+}
+SparseComplexMatrix
+mul_herm (const SparseComplexMatrix& a, const DiagMatrix& d)
+{
+ return operator * (a, d);
+}
+
+SparseComplexMatrix
+operator * (const ComplexDiagMatrix& d, const SparseMatrix& a)
+{
+ return octave_impl::do_mul_dm_sm<SparseComplexMatrix> (d, a);
+}
+SparseComplexMatrix
+trans_mul (const ComplexDiagMatrix& d, const SparseMatrix& a)
+{
+ return operator * (d, a);
+}
+SparseComplexMatrix
+herm_mul (const ComplexDiagMatrix& d, const SparseMatrix& a)
+{
+ return octave_impl::do_mul_dm_sm<SparseComplexMatrix> (d, a,
+ octave_impl::conj_val<ComplexDiagMatrix::element_type> ());
+}
+SparseComplexMatrix
+operator * (const SparseMatrix& a, const ComplexDiagMatrix& d)
+{
+ return octave_impl::do_mul_sm_dm<SparseComplexMatrix> (a, d);
+}
+SparseComplexMatrix
+mul_trans (const SparseMatrix& a, const ComplexDiagMatrix& d)
+{
+ return operator * (a, d);
+}
+SparseComplexMatrix
+mul_herm (const SparseMatrix& a, const ComplexDiagMatrix& d)
+{
+ return octave_impl::do_mul_sm_dm<SparseComplexMatrix> (a, d,
+ octave_impl::conj_val<ComplexDiagMatrix::element_type> ());
+}
+
+SparseComplexMatrix
+operator * (const ComplexDiagMatrix& d, const SparseComplexMatrix& a)
+{
+ return octave_impl::do_mul_dm_sm<SparseComplexMatrix> (d, a);
+}
+SparseComplexMatrix
+trans_mul (const ComplexDiagMatrix& d, const SparseComplexMatrix& a)
+{
+ return operator * (d, a);
+}
+SparseComplexMatrix
+herm_mul (const ComplexDiagMatrix& d, const SparseComplexMatrix& a)
+{
+ return octave_impl::do_mul_dm_sm<SparseComplexMatrix> (d, a,
+ octave_impl::conj_val<ComplexDiagMatrix::element_type> ());
+}
+SparseComplexMatrix
+operator * (const SparseComplexMatrix& a, const ComplexDiagMatrix& d)
+{
+ return octave_impl::do_mul_sm_dm<SparseComplexMatrix> (a, d);
+}
+SparseComplexMatrix
+mul_trans (const SparseComplexMatrix& a, const ComplexDiagMatrix& d)
+{
+ return operator * (a, d);
+}
+SparseComplexMatrix
+mul_herm (const SparseComplexMatrix& a, const ComplexDiagMatrix& d)
+{
+ return octave_impl::do_mul_sm_dm<SparseComplexMatrix> (a, d,
+ octave_impl::conj_val<ComplexDiagMatrix::element_type> ());
+}
+
// FIXME -- it would be nice to share code among the min/max
// functions below.
diff --git a/liboctave/CSparse.h b/liboctave/CSparse.h
--- a/liboctave/CSparse.h
+++ b/liboctave/CSparse.h
@@ -37,6 +37,8 @@
#include "Sparse-op-defs.h"
#include "MatrixType.h"
+class DiagMatrix;
+class ComplexDiagMatrix;
class SparseMatrix;
class SparseBoolMatrix;
@@ -473,6 +475,27 @@
extern OCTAVE_API ComplexMatrix herm_mul (const SparseComplexMatrix&,
const ComplexMatrix&);
+extern OCTAVE_API SparseComplexMatrix operator * (const DiagMatrix&, const
SparseComplexMatrix&);
+extern OCTAVE_API SparseComplexMatrix trans_mul (const DiagMatrix&, const
SparseComplexMatrix&);
+extern OCTAVE_API SparseComplexMatrix herm_mul (const DiagMatrix&, const
SparseComplexMatrix&);
+extern OCTAVE_API SparseComplexMatrix operator * (const SparseComplexMatrix&,
const DiagMatrix&);
+extern OCTAVE_API SparseComplexMatrix mul_trans (const SparseComplexMatrix&,
const DiagMatrix&);
+extern OCTAVE_API SparseComplexMatrix mul_herm (const SparseComplexMatrix&,
const DiagMatrix&);
+
+extern OCTAVE_API SparseComplexMatrix operator * (const ComplexDiagMatrix&,
const SparseMatrix&);
+extern OCTAVE_API SparseComplexMatrix trans_mul (const ComplexDiagMatrix&,
const SparseMatrix&);
+extern OCTAVE_API SparseComplexMatrix herm_mul (const ComplexDiagMatrix&,
const SparseMatrix&);
+extern OCTAVE_API SparseComplexMatrix operator * (const SparseMatrix&, const
ComplexDiagMatrix&);
+extern OCTAVE_API SparseComplexMatrix mul_trans (const SparseMatrix&, const
ComplexDiagMatrix&);
+extern OCTAVE_API SparseComplexMatrix mul_herm (const SparseMatrix&, const
ComplexDiagMatrix&);
+
+extern OCTAVE_API SparseComplexMatrix operator * (const ComplexDiagMatrix&,
const SparseComplexMatrix&);
+extern OCTAVE_API SparseComplexMatrix trans_mul (const ComplexDiagMatrix&,
const SparseComplexMatrix&);
+extern OCTAVE_API SparseComplexMatrix herm_mul (const ComplexDiagMatrix&,
const SparseComplexMatrix&);
+extern OCTAVE_API SparseComplexMatrix operator * (const SparseComplexMatrix&,
const ComplexDiagMatrix&);
+extern OCTAVE_API SparseComplexMatrix mul_trans (const SparseComplexMatrix&,
const ComplexDiagMatrix&);
+extern OCTAVE_API SparseComplexMatrix mul_herm (const SparseComplexMatrix&,
const ComplexDiagMatrix&);
+
extern OCTAVE_API SparseComplexMatrix min (const Complex& c,
const SparseComplexMatrix& m);
extern OCTAVE_API SparseComplexMatrix min (const SparseComplexMatrix& m,
diff --git a/liboctave/ChangeLog b/liboctave/ChangeLog
--- a/liboctave/ChangeLog
+++ b/liboctave/ChangeLog
@@ -1,3 +1,18 @@
+2009-03-08 Jason Riedy <address@hidden>
+
+ * Sparse-diag-op-defs.h (octave_impl::do_mul_dm_sm)
+ (octave_impl::do_mul_sm_dm): New template
+ functions. Implementations for sparse * diag and diag * sparse.
+
+ * CSparse.h (operator *, trans_mul, herm_mul): Add overloads for
+ DiagMatrix and ComplexDiagMatrix.
+ * CSparse.cc (operator *, trans_mul, herm_mul): Implement
+ operations by calling approprate functions in
+ Sparse-diag-op-defs.h.
+ * dSparse.h (operator *, trans_mul): Add overloads for DiagMatrix.
+ * dSparse.cc (operator *, trans_mul): Implement operations by
+ calling approprate functions in Sparse-diag-op-defs.h.
+
2009-03-08 Jason Riedy <address@hidden>
* functor.h (octave_impl::identity_val): New function object.
diff --git a/liboctave/Sparse-diag-op-defs.h b/liboctave/Sparse-diag-op-defs.h
new file mode 100644
--- /dev/null
+++ b/liboctave/Sparse-diag-op-defs.h
@@ -0,0 +1,117 @@
+/* -*- 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_diag_op_defs_h)
+#define octave_sparse_diag_op_defs_h 1
+
+namespace octave_impl {
+
+template <typename RT, typename DM, typename SM, typename UnaryOp>
+RT do_mul_dm_sm (const DM& d, const SM& a, UnaryOp conj_fn)
+{
+ const octave_idx_type nr = d.rows ();
+ const octave_idx_type nc = d.cols ();
+
+ const octave_idx_type a_nr = a.rows ();
+ const octave_idx_type a_nc = a.cols ();
+
+ if (nc != a_nr)
+ {
+ gripe_nonconformant ("operator *", nr, nc, a_nr, a_nc);
+ return RT ();
+ }
+ else
+ {
+ const octave_idx_type nz = a.nnz ();
+ RT r (a_nr, a_nc, nz);
+
+ for (octave_idx_type k = 0; k < nz; k++)
+ {
+ OCTAVE_QUIT;
+ const octave_idx_type i = a.ridx (k);
+ r.xdata (k) = conj_fn (d.dgelem (i)) * a.data (k);
+ r.xridx (k) = i;
+ }
+ for (octave_idx_type j = 0; j < a_nc + 1; j++)
+ {
+ OCTAVE_QUIT;
+ r.xcidx (j) = a.cidx (j);
+ }
+
+ r.maybe_compress (true);
+ return r;
+ }
+}
+
+template <typename RT, typename DM, typename SM>
+RT do_mul_dm_sm (const DM& d, const SM& a)
+{
+ return do_mul_dm_sm<RT> (d, a, octave_impl::identity_val<typename
DM::element_type>());
+}
+
+template <typename RT, typename SM, typename DM, typename UnaryOp>
+RT do_mul_sm_dm (const SM& a, const DM& d, UnaryOp conj_fn)
+{
+ const octave_idx_type nr = d.rows ();
+ const octave_idx_type nc = d.cols ();
+
+ const octave_idx_type a_nr = a.rows ();
+ const octave_idx_type a_nc = a.cols ();
+
+ if (nc != a_nr)
+ {
+ gripe_nonconformant ("operator *", nr, nc, a_nr, a_nc);
+ return RT ();
+ }
+ else
+ {
+ const octave_idx_type nz = a.nnz ();
+ RT r (a_nr, a_nc, nz);
+
+ for (octave_idx_type j = 0; j < nc; ++j)
+ {
+ OCTAVE_QUIT;
+ const typename DM::element_type s = conj_fn (d.dgelem (j));
+ const octave_idx_type colend = a.cidx (j+1);
+ r.xcidx (j) = a.cidx (j);
+ for (octave_idx_type k = a.cidx (j); k < colend; ++k)
+ {
+ r.xdata (k) = s * a.data (k);
+ r.xridx (k) = a.ridx (k);
+ }
+ }
+ r.xcidx (nc) = a.cidx (nc);
+
+ r.maybe_compress (true);
+ return r;
+ }
+}
+
+template <typename RT, typename SM, typename DM>
+RT do_mul_sm_dm (const SM& a, const DM& d)
+{
+ return do_mul_sm_dm<RT> (a, d, octave_impl::identity_val<typename
DM::element_type>());
+}
+
+} // namespace octave_impl
+
+#endif // octave_sparse_diag_op_defs_h
diff --git a/liboctave/dSparse.cc b/liboctave/dSparse.cc
--- a/liboctave/dSparse.cc
+++ b/liboctave/dSparse.cc
@@ -29,6 +29,7 @@
#include <iostream>
#include <vector>
+#include <functional>
#include "quit.h"
#include "lo-ieee.h"
@@ -37,6 +38,7 @@
#include "dRowVector.h"
#include "oct-locbuf.h"
+#include "dDiagMatrix.h"
#include "CSparse.h"
#include "boolSparse.h"
#include "dSparse.h"
@@ -49,6 +51,8 @@
#include "SparsedbleCHOL.h"
#include "SparseQR.h"
+#include "Sparse-diag-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
@@ -7695,6 +7699,32 @@
SPARSE_FULL_TRANS_MUL (Matrix, double, 0., );
}
+// diag * sparse and sparse * diag
+
+SparseMatrix
+operator * (const DiagMatrix& d, const SparseMatrix& a)
+{
+ return octave_impl::do_mul_dm_sm<SparseMatrix> (d, a);
+}
+
+SparseMatrix
+trans_mul (const DiagMatrix& d, const SparseMatrix& a)
+{
+ return operator * (d, a);
+}
+
+SparseMatrix
+operator * (const SparseMatrix& a, const DiagMatrix& d)
+{
+ return octave_impl::do_mul_sm_dm<SparseMatrix> (a, d);
+}
+
+SparseMatrix
+mul_trans (const SparseMatrix& a, const DiagMatrix& d)
+{
+ return operator * (a, d);
+}
+
// FIXME -- it would be nice to share code among the min/max
// functions below.
diff --git a/liboctave/dSparse.h b/liboctave/dSparse.h
--- a/liboctave/dSparse.h
+++ b/liboctave/dSparse.h
@@ -29,7 +29,6 @@
#include "CMatrix.h"
#include "dColVector.h"
#include "CColVector.h"
-#include "dDiagMatrix.h"
#include "DET.h"
#include "MSparse.h"
@@ -37,6 +36,7 @@
#include "Sparse-op-defs.h"
#include "MatrixType.h"
+class DiagMatrix;
class SparseComplexMatrix;
class SparseBoolMatrix;
@@ -451,6 +451,11 @@
extern OCTAVE_API Matrix trans_mul (const SparseMatrix& a,
const Matrix& b);
+extern OCTAVE_API SparseMatrix operator * (const DiagMatrix&, const
SparseMatrix&);
+extern OCTAVE_API SparseMatrix trans_mul (const DiagMatrix&, const
SparseMatrix&);
+extern OCTAVE_API SparseMatrix operator * (const SparseMatrix&, const
DiagMatrix&);
+extern OCTAVE_API SparseMatrix mul_trans (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 --git a/src/ChangeLog b/src/ChangeLog
--- a/src/ChangeLog
+++ b/src/ChangeLog
@@ -1,3 +1,12 @@
+2009-03-08 Jason Riedy <address@hidden>
+
+ * OPERATORS/op-dm-scm.cc: New file. Implement multiplication
+ between diagonal matrices (both real and complex) and complex
+ sparse matrices.
+ * OPERATORS/op-dm-sm.cc: New file. Implement multiplication
+ between diagonal matrices and sparse matrices, all real.
+ * Makefile.in (DIAG_OP_XSRC): Add op-dm-sm.cc and op-dm-scm.cc.
+
2009-03-09 Jaroslav Hajek <address@hidden>
* data.cc (F__accumarray_sum__): New function.
diff --git a/src/Makefile.in b/src/Makefile.in
--- a/src/Makefile.in
+++ b/src/Makefile.in
@@ -166,7 +166,7 @@
DIAG_OP_XSRC := op-cdm-cdm.cc op-cdm-cm.cc op-cdm-cs.cc op-cdm-dm.cc \
op-cdm-m.cc op-cdm-s.cc op-cm-cdm.cc op-cm-dm.cc op-dm-cdm.cc \
op-dm-cm.cc op-dm-cs.cc op-dm-dm.cc op-dm-m.cc op-dm-s.cc \
- op-m-cdm.cc op-m-dm.cc
+ op-m-cdm.cc op-m-dm.cc op-dm-sm.cc op-dm-scm.cc
FDIAG_OP_XSRC := op-fcdm-fcdm.cc op-fcdm-fcm.cc op-fcdm-fcs.cc op-fcdm-fdm.cc \
op-fcdm-fm.cc op-fcdm-fs.cc op-fcm-fcdm.cc op-fcm-fdm.cc \
diff --git a/src/OPERATORS/op-dm-scm.cc b/src/OPERATORS/op-dm-scm.cc
new file mode 100644
--- /dev/null
+++ b/src/OPERATORS/op-dm-scm.cc
@@ -0,0 +1,403 @@
+/*
+
+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-re-diag.h"
+#include "ov-cx-diag.h"
+#include "ov-re-sparse.h"
+#include "ov-cx-sparse.h"
+
+// diagonal matrix by sparse matrix ops
+
+DEFBINOP (mul_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<double> d = v2.complex_value ();
+
+ return octave_value (v1.diag_matrix_value () * d);
+ }
+ else if (v1.rows() == 1 && v1.columns() == 1)
+ // If v1 is a scalar in disguise, don't bother with further dispatching.
+ {
+ double d = v1.scalar_value ();
+
+ return octave_value (d * v1.sparse_complex_matrix_value ());
+ }
+ else
+ {
+
+ MatrixType typ = v2.matrix_type ();
+ SparseComplexMatrix ret = v1.diag_matrix_value () *
v2.sparse_complex_matrix_value ();
+ octave_value out = octave_value (ret);
+ typ.mark_as_unsymmetric ();
+ out.matrix_type (typ);
+ return out;
+ }
+}
+
+DEFBINOP (mul_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.
+ {
+ std::complex<double> d = v2.scalar_value ();
+
+ return octave_value (v1.complex_diag_matrix_value () * d);
+ }
+ else if (v1.rows() == 1 && v1.columns() == 1)
+ // If v1 is a scalar in disguise, don't bother with further dispatching.
+ {
+ std::complex<double> d = v1.complex_value ();
+
+ return octave_value (d * v1.sparse_complex_matrix_value ());
+ }
+ else
+ {
+
+ MatrixType typ = v2.matrix_type ();
+ SparseComplexMatrix ret = v1.complex_diag_matrix_value () *
v2.sparse_matrix_value ();
+ octave_value out = octave_value (ret);
+ typ.mark_as_unsymmetric ();
+ out.matrix_type (typ);
+ return out;
+ }
+}
+
+DEFBINOP (herm_mul_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.
+ {
+ std::complex<double> d = v2.scalar_value ();
+
+ return octave_value (v1.complex_diag_matrix_value () * d);
+ }
+ else if (v1.rows() == 1 && v1.columns() == 1)
+ // If v1 is a scalar in disguise, don't bother with further dispatching.
+ {
+ using std::conj;
+ std::complex<double> d = conj (v1.complex_value ());
+
+ return octave_value (d * v1.sparse_complex_matrix_value ());
+ }
+ else
+ {
+
+ MatrixType typ = v2.matrix_type ();
+ SparseComplexMatrix ret = herm_mul (v1.complex_diag_matrix_value (),
+ v2.sparse_matrix_value ());
+ octave_value out = octave_value (ret);
+ typ.mark_as_unsymmetric ();
+ out.matrix_type (typ);
+ return out;
+ }
+}
+
+DEFBINOP (mul_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<double> d = v2.complex_value ();
+
+ return octave_value (v1.complex_diag_matrix_value () * d);
+ }
+ else if (v1.rows() == 1 && v1.columns() == 1)
+ // If v1 is a scalar in disguise, don't bother with further dispatching.
+ {
+ std::complex<double> d = v1.complex_value ();
+
+ return octave_value (d * v1.sparse_complex_matrix_value ());
+ }
+ else
+ {
+
+ MatrixType typ = v2.matrix_type ();
+ SparseComplexMatrix ret = v1.complex_diag_matrix_value () *
v2.sparse_complex_matrix_value ();
+ octave_value out = octave_value (ret);
+ typ.mark_as_unsymmetric ();
+ out.matrix_type (typ);
+ return out;
+ }
+}
+
+DEFBINOP (herm_mul_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<double> d = v2.complex_value ();
+
+ return octave_value (v1.complex_diag_matrix_value () * d);
+ }
+ else if (v1.rows() == 1 && v1.columns() == 1)
+ // If v1 is a scalar in disguise, don't bother with further dispatching.
+ {
+ using std::conj;
+ std::complex<double> d = conj (v1.complex_value ());
+
+ return octave_value (d * v1.sparse_complex_matrix_value ());
+ }
+ else
+ {
+
+ MatrixType typ = v2.matrix_type ();
+ SparseComplexMatrix ret = herm_mul (v1.complex_diag_matrix_value (),
+ v2.sparse_complex_matrix_value ());
+ octave_value out = octave_value (ret);
+ typ.mark_as_unsymmetric ();
+ out.matrix_type (typ);
+ return out;
+ }
+}
+
+// sparse matrix by diagonal matrix ops
+
+DEFBINOP (mul_scm_dm, sparse_complex_matrix, diag_matrix)
+{
+ CAST_BINOP_ARGS (const octave_sparse_complex_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.
+ {
+ std::complex<double> d = v1.complex_value ();
+
+ return octave_value (d * v2.diag_matrix_value ());
+ }
+ else if (v2.rows() == 1 && v2.columns() == 1)
+ // If v2 is a scalar in disguise, don't bother with further dispatching.
+ {
+ double d = v2.scalar_value ();
+
+ return octave_value (v1.sparse_complex_matrix_value () * d);
+ }
+ else
+ {
+
+ MatrixType typ = v1.matrix_type ();
+ SparseComplexMatrix ret = v1.sparse_complex_matrix_value () *
v2.diag_matrix_value ();
+ octave_value out = octave_value (ret);
+ typ.mark_as_unsymmetric ();
+ out.matrix_type (typ);
+ return out;
+ }
+}
+
+DEFBINOP (mul_sm_cdm, sparse_matrix, complex_diag_matrix)
+{
+ CAST_BINOP_ARGS (const octave_sparse_matrix&, const
octave_complex_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.
+ {
+ std::complex<double> d = v1.complex_value ();
+
+ return octave_value (d * v2.complex_diag_matrix_value ());
+ }
+ else if (v2.rows() == 1 && v2.columns() == 1)
+ // If v2 is a scalar in disguise, don't bother with further dispatching.
+ {
+ std::complex<double> d = v2.complex_value ();
+
+ return octave_value (v1.sparse_complex_matrix_value () * d);
+ }
+ else
+ {
+
+ MatrixType typ = v1.matrix_type ();
+ SparseComplexMatrix ret = v1.sparse_matrix_value () *
v2.complex_diag_matrix_value ();
+ octave_value out = octave_value (ret);
+ typ.mark_as_unsymmetric ();
+ out.matrix_type (typ);
+ return out;
+ }
+}
+
+DEFBINOP (mul_herm_sm_cdm, sparse_matrix, complex_diag_matrix)
+{
+ CAST_BINOP_ARGS (const octave_sparse_matrix&, const
octave_complex_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.
+ {
+ std::complex<double> d = v1.complex_value ();
+
+ return octave_value (d * v2.complex_diag_matrix_value ());
+ }
+ else if (v2.rows() == 1 && v2.columns() == 1)
+ // If v2 is a scalar in disguise, don't bother with further dispatching.
+ {
+ using std::conj;
+ std::complex<double> d = conj (v2.complex_value ());
+
+ return octave_value (v1.sparse_complex_matrix_value () * d);
+ }
+ else
+ {
+
+ MatrixType typ = v1.matrix_type ();
+ SparseComplexMatrix ret = mul_herm (v1.sparse_matrix_value (),
+ v2.complex_diag_matrix_value ());
+ octave_value out = octave_value (ret);
+ typ.mark_as_unsymmetric ();
+ out.matrix_type (typ);
+ return out;
+ }
+}
+
+DEFBINOP (mul_scm_cdm, sparse_complex_matrix, complex_diag_matrix)
+{
+ CAST_BINOP_ARGS (const octave_sparse_complex_matrix&, const
octave_complex_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.
+ {
+ std::complex<double> d = v1.complex_value ();
+
+ return octave_value (d * v2.complex_diag_matrix_value ());
+ }
+ else if (v2.rows() == 1 && v2.columns() == 1)
+ // If v2 is a scalar in disguise, don't bother with further dispatching.
+ {
+ std::complex<double> d = v2.complex_value ();
+
+ return octave_value (v1.sparse_complex_matrix_value () * d);
+ }
+ else
+ {
+
+ MatrixType typ = v1.matrix_type ();
+ SparseComplexMatrix ret = v1.sparse_complex_matrix_value () *
v2.complex_diag_matrix_value ();
+ octave_value out = octave_value (ret);
+ typ.mark_as_unsymmetric ();
+ out.matrix_type (typ);
+ return out;
+ }
+}
+
+DEFBINOP (mul_herm_scm_cdm, sparse_complex_matrix, complex_diag_matrix)
+{
+ CAST_BINOP_ARGS (const octave_sparse_complex_matrix&, const
octave_complex_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.
+ {
+ std::complex<double> d = v1.complex_value ();
+
+ return octave_value (d * v2.complex_diag_matrix_value ());
+ }
+ else if (v2.rows() == 1 && v2.columns() == 1)
+ // If v2 is a scalar in disguise, don't bother with further dispatching.
+ {
+ std::complex<double> d = std::conj (v2.complex_value ());
+
+ return octave_value (v1.sparse_complex_matrix_value () * d);
+ }
+ else
+ {
+
+ MatrixType typ = v1.matrix_type ();
+ SparseComplexMatrix ret = mul_herm (v1.sparse_complex_matrix_value (),
+ v2.complex_diag_matrix_value ());
+ octave_value out = octave_value (ret);
+ typ.mark_as_unsymmetric ();
+ out.matrix_type (typ);
+ return out;
+ }
+}
+
+void
+install_dm_scm_ops (void)
+{
+ INSTALL_BINOP (op_mul, octave_diag_matrix, octave_sparse_complex_matrix,
+ mul_dm_scm);
+ INSTALL_BINOP (op_trans_mul, octave_diag_matrix,
octave_sparse_complex_matrix,
+ mul_dm_scm);
+ INSTALL_BINOP (op_herm_mul, octave_diag_matrix, octave_sparse_complex_matrix,
+ mul_dm_scm);
+
+ INSTALL_BINOP (op_mul, octave_complex_diag_matrix, octave_sparse_matrix,
+ mul_cdm_sm);
+ INSTALL_BINOP (op_trans_mul, octave_complex_diag_matrix,
octave_sparse_matrix,
+ mul_cdm_sm);
+ INSTALL_BINOP (op_herm_mul, octave_complex_diag_matrix, octave_sparse_matrix,
+ herm_mul_cdm_sm);
+
+ INSTALL_BINOP (op_mul, octave_complex_diag_matrix,
octave_sparse_complex_matrix,
+ mul_cdm_scm);
+ INSTALL_BINOP (op_trans_mul, octave_complex_diag_matrix,
octave_sparse_complex_matrix,
+ mul_cdm_scm);
+ INSTALL_BINOP (op_herm_mul, octave_complex_diag_matrix,
octave_sparse_complex_matrix,
+ herm_mul_cdm_scm);
+
+ INSTALL_BINOP (op_mul, octave_sparse_complex_matrix, octave_diag_matrix,
+ mul_scm_dm);
+ INSTALL_BINOP (op_mul_trans, octave_sparse_complex_matrix,
octave_diag_matrix,
+ mul_scm_dm);
+ INSTALL_BINOP (op_mul_herm, 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_trans, octave_sparse_matrix,
octave_complex_diag_matrix,
+ mul_sm_cdm);
+ INSTALL_BINOP (op_mul_herm, octave_sparse_matrix, octave_complex_diag_matrix,
+ mul_herm_sm_cdm);
+
+ INSTALL_BINOP (op_mul, octave_sparse_complex_matrix,
octave_complex_diag_matrix,
+ mul_scm_cdm);
+ INSTALL_BINOP (op_mul_trans, octave_sparse_complex_matrix,
octave_complex_diag_matrix,
+ mul_scm_cdm);
+ INSTALL_BINOP (op_mul_herm, octave_sparse_complex_matrix,
octave_complex_diag_matrix,
+ mul_herm_scm_cdm);
+}
diff --git a/src/OPERATORS/op-dm-sm.cc b/src/OPERATORS/op-dm-sm.cc
new file mode 100644
--- /dev/null
+++ b/src/OPERATORS/op-dm-sm.cc
@@ -0,0 +1,118 @@
+/*
+
+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-re-diag.h"
+#include "ov-re-sparse.h"
+
+// diagonal matrix by sparse matrix ops
+
+DEFBINOP (mul_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 if (v1.rows() == 1 && v1.columns() == 1)
+ // If v1 is a scalar in disguise, don't bother with further dispatching.
+ {
+ double d = v1.scalar_value ();
+
+ return octave_value (d * v1.sparse_matrix_value ());
+ }
+ else
+ {
+
+ MatrixType typ = v2.matrix_type ();
+ SparseMatrix ret = v1.diag_matrix_value () * v2.sparse_matrix_value ();
+ octave_value out = octave_value (ret);
+ typ.mark_as_unsymmetric ();
+ out.matrix_type (typ);
+ return out;
+ }
+}
+
+// sparse matrix by diagonal matrix ops
+
+DEFBINOP (mul_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 if (v2.rows() == 1 && v2.columns() == 1)
+ // If v2 is a scalar in disguise, don't bother with further dispatching.
+ {
+ double d = v2.scalar_value ();
+
+ return octave_value (v1.sparse_matrix_value () * d);
+ }
+ else
+ {
+
+ MatrixType typ = v1.matrix_type ();
+ SparseMatrix ret = v1.sparse_matrix_value () * v2.diag_matrix_value ();
+ octave_value out = octave_value (ret);
+ typ.mark_as_unsymmetric ();
+ out.matrix_type (typ);
+ return out;
+ }
+}
+
+void
+install_dm_sm_ops (void)
+{
+ INSTALL_BINOP (op_mul, octave_diag_matrix, octave_sparse_matrix,
+ mul_dm_sm);
+ INSTALL_BINOP (op_trans_mul, octave_diag_matrix, octave_sparse_matrix,
+ mul_dm_sm);
+ INSTALL_BINOP (op_herm_mul, octave_diag_matrix, octave_sparse_matrix,
+ mul_dm_sm);
+
+ INSTALL_BINOP (op_mul, octave_sparse_matrix, octave_diag_matrix,
+ mul_sm_dm);
+ INSTALL_BINOP (op_trans_mul, octave_sparse_matrix, octave_diag_matrix,
+ mul_sm_dm);
+ INSTALL_BINOP (op_herm_mul, octave_sparse_matrix, octave_diag_matrix,
+ mul_sm_dm);
+}
diff --git a/test/ChangeLog b/test/ChangeLog
--- a/test/ChangeLog
+++ b/test/ChangeLog
@@ -1,3 +1,7 @@
+2009-03-08 Jason Riedy <address@hidden>
+
+ * test_diag_perm.m: Add tests for preserving sparse structure when
scaling.
+
2009-02-25 John W. Eaton <address@hidden>
* build_sparse_tests.sh: Note that saving sparse matrices to MAT
diff --git a/test/test_diag_perm.m b/test/test_diag_perm.m
--- a/test/test_diag_perm.m
+++ b/test/test_diag_perm.m
@@ -124,10 +124,17 @@
%! Dslice = D (1:(m-3), 1:(n-2));
%! assert (typeinfo (Dslice), "diagonal matrix");
-## preserve dense matrix structure
+## preserve dense matrix structure when scaling
%!assert (typeinfo (rand (8) * (3 * eye (8))), "matrix");
%!assert (typeinfo ((3 * eye (8)) * rand (8)), "matrix");
+## preserve sparse matrix structure when scaling
+%!assert (typeinfo (sprand (8, 8, .5) * (3 * eye (8))), "sparse matrix");
+%!assert (typeinfo (sprand (8, 8, .5) * (3 * eye (8))'), "sparse matrix");
+%!assert (typeinfo (((3 + 2 * I ()) * eye (8)) * sprand (8, 8, .5)), "sparse
complex matrix");
+%!assert (typeinfo (((3 + 2 * I ()) * eye (8))' * sprand (8, 8, .5)), "sparse
complex matrix");
+%!assert (typeinfo (sprand (8, 8, .5) * ((3 + 2 * I ()) * eye (8)).'), "sparse
complex matrix");
+
## scaling a matrix with exceptional values does not introduce new ones.
%!test
%! n = 6;
- [PATCH 0 of 4] Implement basic sparse op diag and diag op sparse support., Jason Riedy, 2009/03/09
- [PATCH 1 of 4] Add identity_val and conj_val functor objects to functor.h, Jason Riedy, 2009/03/09
- [PATCH 2 of 4] Implement sparse * diagonal and diagonal * sparse operations, double-prec only,
Jason Riedy <=
- [PATCH 3 of 4] Implement diag \ sparse and sparse / diag, Jason Riedy, 2009/03/09
- [PATCH 4 of 4] Implement diag + sparse, diag - sparse, sparse + diag, sparse - diag, Jason Riedy, 2009/03/09
- [PATCH 4 of 4] Implement diag + sparse, diag - sparse, sparse + diag, sparse - diag, John W. Eaton, 2009/03/09
- Re: [PATCH 4 of 4] Implement diag + sparse, diag - sparse, sparse + diag, sparse - diag, Jason Riedy, 2009/03/10
- Re: [PATCH 4 of 4] Implement diag + sparse, diag - sparse, sparse + diag, sparse - diag, John W. Eaton, 2009/03/10
- Re: [PATCH 4 of 4] Implement diag + sparse, diag - sparse, sparse + diag, sparse - diag, Jason Riedy, 2009/03/10
- Re: [PATCH 4 of 4] Implement diag + sparse, diag - sparse, sparse + diag, sparse - diag, Jaroslav Hajek, 2009/03/11
[PATCH 0 of 4] Implement basic sparse op diag and diag op sparse support., John W. Eaton, 2009/03/09