Mercurial > octave-libgccjit
comparison src/OPERATORS/op-dm-sm.cc @ 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 | |
children | 42aff15e059b |
comparison
equal
deleted
inserted
replaced
8963:d1eab3ddb02d | 8964:f4f4d65faaa0 |
---|---|
1 /* | |
2 | |
3 Copyright (C) 2009 Jason Riedy | |
4 | |
5 This file is part of Octave. | |
6 | |
7 Octave is free software; you can redistribute it and/or modify it | |
8 under the terms of the GNU General Public License as published by the | |
9 Free Software Foundation; either version 3 of the License, or (at your | |
10 option) any later version. | |
11 | |
12 Octave is distributed in the hope that it will be useful, but WITHOUT | |
13 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or | |
14 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License | |
15 for more details. | |
16 | |
17 You should have received a copy of the GNU General Public License | |
18 along with Octave; see the file COPYING. If not, see | |
19 <http://www.gnu.org/licenses/>. | |
20 | |
21 */ | |
22 | |
23 #ifdef HAVE_CONFIG_H | |
24 #include <config.h> | |
25 #endif | |
26 | |
27 #include "gripes.h" | |
28 #include "oct-obj.h" | |
29 #include "ov.h" | |
30 #include "ov-typeinfo.h" | |
31 #include "ops.h" | |
32 | |
33 #include "ov-re-diag.h" | |
34 #include "ov-re-sparse.h" | |
35 | |
36 // diagonal matrix by sparse matrix ops | |
37 | |
38 DEFBINOP (mul_dm_sm, diag_matrix, sparse_matrix) | |
39 { | |
40 CAST_BINOP_ARGS (const octave_diag_matrix&, const octave_sparse_matrix&); | |
41 | |
42 if (v2.rows() == 1 && v2.columns() == 1) | |
43 // If v2 is a scalar in disguise, return a diagonal matrix rather than | |
44 // a sparse matrix. | |
45 { | |
46 double d = v2.scalar_value (); | |
47 | |
48 return octave_value (v1.diag_matrix_value () * d); | |
49 } | |
50 else | |
51 { | |
52 MatrixType typ = v2.matrix_type (); | |
53 SparseMatrix ret = v1.diag_matrix_value () * v2.sparse_matrix_value (); | |
54 octave_value out = octave_value (ret); | |
55 typ.mark_as_unsymmetric (); | |
56 out.matrix_type (typ); | |
57 return out; | |
58 } | |
59 } | |
60 | |
61 // sparse matrix by diagonal matrix ops | |
62 | |
63 DEFBINOP (mul_sm_dm, sparse_matrix, diag_matrix) | |
64 { | |
65 CAST_BINOP_ARGS (const octave_sparse_matrix&, const octave_diag_matrix&); | |
66 | |
67 if (v1.rows() == 1 && v1.columns() == 1) | |
68 // If v1 is a scalar in disguise, return a diagonal matrix rather than | |
69 // a sparse matrix. | |
70 { | |
71 double d = v1.scalar_value (); | |
72 | |
73 return octave_value (d * v2.diag_matrix_value ()); | |
74 } | |
75 else | |
76 { | |
77 MatrixType typ = v1.matrix_type (); | |
78 SparseMatrix ret = v1.sparse_matrix_value () * v2.diag_matrix_value (); | |
79 octave_value out = octave_value (ret); | |
80 typ.mark_as_unsymmetric (); | |
81 out.matrix_type (typ); | |
82 return out; | |
83 } | |
84 } | |
85 | |
86 void | |
87 install_dm_sm_ops (void) | |
88 { | |
89 INSTALL_BINOP (op_mul, octave_diag_matrix, octave_sparse_matrix, | |
90 mul_dm_sm); | |
91 | |
92 INSTALL_BINOP (op_mul, octave_sparse_matrix, octave_diag_matrix, | |
93 mul_sm_dm); | |
94 } |