Mercurial > octave-nkf
changeset 8964:f4f4d65faaa0
Implement sparse * diagonal and diagonal * sparse operations, double-prec only.
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.
| author | Jason Riedy <jason@acm.org> |
|---|---|
| date | Mon, 09 Mar 2009 17:49:13 -0400 |
| parents | d1eab3ddb02d |
| children | 42aff15e059b |
| files | liboctave/CSparse.cc liboctave/CSparse.h liboctave/ChangeLog liboctave/Sparse-diag-op-defs.h liboctave/dSparse.cc liboctave/dSparse.h src/ChangeLog src/Makefile.in src/OPERATORS/op-dm-scm.cc src/OPERATORS/op-dm-sm.cc test/ChangeLog test/test_diag_perm.m |
| diffstat | 12 files changed, 510 insertions(+), 3 deletions(-) [+] |
line wrap: on
line diff
--- a/liboctave/CSparse.cc Thu Mar 12 17:08:14 2009 -0400 +++ b/liboctave/CSparse.cc Mon Mar 09 17:49:13 2009 -0400 @@ -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,40 @@ 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 +operator * (const SparseComplexMatrix& a, const DiagMatrix& d) +{ + return octave_impl::do_mul_sm_dm<SparseComplexMatrix> (a, d); +} + +SparseComplexMatrix +operator * (const ComplexDiagMatrix& d, const SparseMatrix& a) +{ + return octave_impl::do_mul_dm_sm<SparseComplexMatrix> (d, a); +} +SparseComplexMatrix +operator * (const SparseMatrix& a, const ComplexDiagMatrix& d) +{ + return octave_impl::do_mul_sm_dm<SparseComplexMatrix> (a, d); +} + +SparseComplexMatrix +operator * (const ComplexDiagMatrix& d, const SparseComplexMatrix& a) +{ + return octave_impl::do_mul_dm_sm<SparseComplexMatrix> (d, a); +} +SparseComplexMatrix +operator * (const SparseComplexMatrix& a, const ComplexDiagMatrix& d) +{ + return octave_impl::do_mul_sm_dm<SparseComplexMatrix> (a, d); +} + // FIXME -- it would be nice to share code among the min/max // functions below.
--- a/liboctave/CSparse.h Thu Mar 12 17:08:14 2009 -0400 +++ b/liboctave/CSparse.h Mon Mar 09 17:49:13 2009 -0400 @@ -37,6 +37,8 @@ #include "Sparse-op-defs.h" #include "MatrixType.h" +class DiagMatrix; +class ComplexDiagMatrix; class SparseMatrix; class SparseBoolMatrix; @@ -473,6 +475,15 @@ extern OCTAVE_API ComplexMatrix herm_mul (const SparseComplexMatrix&, const ComplexMatrix&); +extern OCTAVE_API SparseComplexMatrix operator * (const DiagMatrix&, const SparseComplexMatrix&); +extern OCTAVE_API SparseComplexMatrix operator * (const SparseComplexMatrix&, const DiagMatrix&); + +extern OCTAVE_API SparseComplexMatrix operator * (const ComplexDiagMatrix&, const SparseMatrix&); +extern OCTAVE_API SparseComplexMatrix operator * (const SparseMatrix&, const ComplexDiagMatrix&); + +extern OCTAVE_API SparseComplexMatrix operator * (const ComplexDiagMatrix&, const SparseComplexMatrix&); +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,
--- a/liboctave/ChangeLog Thu Mar 12 17:08:14 2009 -0400 +++ b/liboctave/ChangeLog Mon Mar 09 17:49:13 2009 -0400 @@ -1,3 +1,18 @@ +2009-03-08 Jason Riedy <jason@acm.org> + + * 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-12 John W. Eaton <jwe@octave.org> * oct-inttypes.h (bitshift): Apply mask even if not shifting.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/Sparse-diag-op-defs.h Mon Mar 09 17:49:13 2009 -0400 @@ -0,0 +1,106 @@ +/* -*- C++ -*- + +Copyright (C) 2009 Jason Riedy, Jaroslav Hajek + +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 + +template <typename RT, typename DM, typename SM> +RT do_mul_dm_sm (const DM& d, const SM& a) +{ + 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 + { + RT r (nr, a_nc, a.nnz ()); + + octave_idx_type l = 0; + + for (octave_idx_type j = 0; j < a_nc; j++) + { + r.xcidx (j) = l; + const octave_idx_type colend = a.cidx (j+1); + for (octave_idx_type k = a.cidx (j); k < colend; k++) + { + const octave_idx_type i = a.ridx (k); + if (i >= nr) break; + r.xdata (l) = d.dgelem (i) * a.data (k); + r.xridx (l) = i; + l++; + } + } + + r.xcidx (a_nc) = l; + + r.maybe_compress (true); + return r; + } +} + +template <typename RT, typename SM, typename DM> +RT do_mul_sm_dm (const SM& a, const DM& d) +{ + 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 (nr != a_nc) + { + gripe_nonconformant ("operator *", a_nr, a_nc, nr, nc); + return RT (); + } + else + { + + const octave_idx_type mnc = nc < a_nc ? nc: a_nc; + RT r (a_nr, nc, a.cidx (mnc)); + + for (octave_idx_type j = 0; j < mnc; ++j) + { + const typename DM::element_type s = 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); + } + } + for (octave_idx_type j = mnc; j <= nc; ++j) + r.xcidx (j) = a.cidx (mnc); + + r.maybe_compress (true); + return r; + } +} + +#endif // octave_sparse_diag_op_defs_h
--- a/liboctave/dSparse.cc Thu Mar 12 17:08:14 2009 -0400 +++ b/liboctave/dSparse.cc Mon Mar 09 17:49:13 2009 -0400 @@ -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 @@ -7698,6 +7702,20 @@ 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 +operator * (const SparseMatrix& a, const DiagMatrix& d) +{ + return octave_impl::do_mul_sm_dm<SparseMatrix> (a, d); +} + // FIXME -- it would be nice to share code among the min/max // functions below.
--- a/liboctave/dSparse.h Thu Mar 12 17:08:14 2009 -0400 +++ b/liboctave/dSparse.h Mon Mar 09 17:49:13 2009 -0400 @@ -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,9 @@ 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 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);
--- a/src/ChangeLog Thu Mar 12 17:08:14 2009 -0400 +++ b/src/ChangeLog Mon Mar 09 17:49:13 2009 -0400 @@ -99,6 +99,15 @@ * Makefile.in (maintainer-clean): Remove y.tab.h here. (distclean): Not here. +2009-03-08 Jason Riedy <jason@acm.org> + + * 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 <highegg@gmail.com> * data.cc (F__accumarray_sum__): New function.
--- a/src/Makefile.in Thu Mar 12 17:08:14 2009 -0400 +++ b/src/Makefile.in Mon Mar 09 17:49:13 2009 -0400 @@ -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 \
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/OPERATORS/op-dm-scm.cc Mon Mar 09 17:49:13 2009 -0400 @@ -0,0 +1,202 @@ +/* + +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 + { + 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 + { + 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 (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 + { + 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; + } +} + +// 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 + { + 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 + { + 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_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; + } +} + +void +install_dm_scm_ops (void) +{ + INSTALL_BINOP (op_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_mul, octave_complex_diag_matrix, octave_sparse_complex_matrix, + mul_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); +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/OPERATORS/op-dm-sm.cc Mon Mar 09 17:49:13 2009 -0400 @@ -0,0 +1,94 @@ +/* + +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 + { + 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 + { + 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_mul, octave_sparse_matrix, octave_diag_matrix, + mul_sm_dm); +}
--- a/test/ChangeLog Thu Mar 12 17:08:14 2009 -0400 +++ b/test/ChangeLog Mon Mar 09 17:49:13 2009 -0400 @@ -6,6 +6,10 @@ * test_diag_perm.m: Add a test for conversion to sparse form. +2009-03-08 Jason Riedy <jason@acm.org> + + * test_diag_perm.m: Add tests for preserving sparse structure when scaling. + 2009-02-25 John W. Eaton <jwe@octave.org> * build_sparse_tests.sh: Note that saving sparse matrices to MAT
--- a/test/test_diag_perm.m Thu Mar 12 17:08:14 2009 -0400 +++ b/test/test_diag_perm.m Mon Mar 09 17:49:13 2009 -0400 @@ -134,10 +134,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;