Mercurial > octave-libgccjit
comparison src/OPERATORS/op-dm-scm.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-cx-diag.h" | |
35 #include "ov-re-sparse.h" | |
36 #include "ov-cx-sparse.h" | |
37 | |
38 // diagonal matrix by sparse matrix ops | |
39 | |
40 DEFBINOP (mul_dm_scm, diag_matrix, sparse_complex_matrix) | |
41 { | |
42 CAST_BINOP_ARGS (const octave_diag_matrix&, const octave_sparse_complex_matrix&); | |
43 | |
44 if (v2.rows() == 1 && v2.columns() == 1) | |
45 // If v2 is a scalar in disguise, return a diagonal matrix rather than | |
46 // a sparse matrix. | |
47 { | |
48 std::complex<double> d = v2.complex_value (); | |
49 | |
50 return octave_value (v1.diag_matrix_value () * d); | |
51 } | |
52 else | |
53 { | |
54 MatrixType typ = v2.matrix_type (); | |
55 SparseComplexMatrix ret = v1.diag_matrix_value () * v2.sparse_complex_matrix_value (); | |
56 octave_value out = octave_value (ret); | |
57 typ.mark_as_unsymmetric (); | |
58 out.matrix_type (typ); | |
59 return out; | |
60 } | |
61 } | |
62 | |
63 DEFBINOP (mul_cdm_sm, complex_diag_matrix, sparse_matrix) | |
64 { | |
65 CAST_BINOP_ARGS (const octave_complex_diag_matrix&, const octave_sparse_matrix&); | |
66 | |
67 if (v2.rows() == 1 && v2.columns() == 1) | |
68 // If v2 is a scalar in disguise, return a diagonal matrix rather than | |
69 // a sparse matrix. | |
70 { | |
71 std::complex<double> d = v2.scalar_value (); | |
72 | |
73 return octave_value (v1.complex_diag_matrix_value () * d); | |
74 } | |
75 else | |
76 { | |
77 MatrixType typ = v2.matrix_type (); | |
78 SparseComplexMatrix ret = v1.complex_diag_matrix_value () * v2.sparse_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 DEFBINOP (mul_cdm_scm, complex_diag_matrix, sparse_complex_matrix) | |
87 { | |
88 CAST_BINOP_ARGS (const octave_complex_diag_matrix&, const octave_sparse_complex_matrix&); | |
89 | |
90 if (v2.rows() == 1 && v2.columns() == 1) | |
91 // If v2 is a scalar in disguise, return a diagonal matrix rather than | |
92 // a sparse matrix. | |
93 { | |
94 std::complex<double> d = v2.complex_value (); | |
95 | |
96 return octave_value (v1.complex_diag_matrix_value () * d); | |
97 } | |
98 else | |
99 { | |
100 MatrixType typ = v2.matrix_type (); | |
101 SparseComplexMatrix ret = v1.complex_diag_matrix_value () * v2.sparse_complex_matrix_value (); | |
102 octave_value out = octave_value (ret); | |
103 typ.mark_as_unsymmetric (); | |
104 out.matrix_type (typ); | |
105 return out; | |
106 } | |
107 } | |
108 | |
109 // sparse matrix by diagonal matrix ops | |
110 | |
111 DEFBINOP (mul_scm_dm, sparse_complex_matrix, diag_matrix) | |
112 { | |
113 CAST_BINOP_ARGS (const octave_sparse_complex_matrix&, const octave_diag_matrix&); | |
114 | |
115 if (v1.rows() == 1 && v1.columns() == 1) | |
116 // If v1 is a scalar in disguise, return a diagonal matrix rather than | |
117 // a sparse matrix. | |
118 { | |
119 std::complex<double> d = v1.complex_value (); | |
120 | |
121 return octave_value (d * v2.diag_matrix_value ()); | |
122 } | |
123 else | |
124 { | |
125 MatrixType typ = v1.matrix_type (); | |
126 SparseComplexMatrix ret = v1.sparse_complex_matrix_value () * v2.diag_matrix_value (); | |
127 octave_value out = octave_value (ret); | |
128 typ.mark_as_unsymmetric (); | |
129 out.matrix_type (typ); | |
130 return out; | |
131 } | |
132 } | |
133 | |
134 DEFBINOP (mul_sm_cdm, sparse_matrix, complex_diag_matrix) | |
135 { | |
136 CAST_BINOP_ARGS (const octave_sparse_matrix&, const octave_complex_diag_matrix&); | |
137 | |
138 if (v1.rows() == 1 && v1.columns() == 1) | |
139 // If v1 is a scalar in disguise, return a diagonal matrix rather than | |
140 // a sparse matrix. | |
141 { | |
142 std::complex<double> d = v1.complex_value (); | |
143 | |
144 return octave_value (d * v2.complex_diag_matrix_value ()); | |
145 } | |
146 else | |
147 { | |
148 MatrixType typ = v1.matrix_type (); | |
149 SparseComplexMatrix ret = v1.sparse_matrix_value () * v2.complex_diag_matrix_value (); | |
150 octave_value out = octave_value (ret); | |
151 typ.mark_as_unsymmetric (); | |
152 out.matrix_type (typ); | |
153 return out; | |
154 } | |
155 } | |
156 | |
157 DEFBINOP (mul_scm_cdm, sparse_complex_matrix, complex_diag_matrix) | |
158 { | |
159 CAST_BINOP_ARGS (const octave_sparse_complex_matrix&, const octave_complex_diag_matrix&); | |
160 | |
161 if (v1.rows() == 1 && v1.columns() == 1) | |
162 // If v1 is a scalar in disguise, return a diagonal matrix rather than | |
163 // a sparse matrix. | |
164 { | |
165 std::complex<double> d = v1.complex_value (); | |
166 | |
167 return octave_value (d * v2.complex_diag_matrix_value ()); | |
168 } | |
169 else if (v2.rows() == 1 && v2.columns() == 1) | |
170 // If v2 is a scalar in disguise, don't bother with further dispatching. | |
171 { | |
172 std::complex<double> d = v2.complex_value (); | |
173 | |
174 return octave_value (v1.sparse_complex_matrix_value () * d); | |
175 } | |
176 else | |
177 { | |
178 MatrixType typ = v1.matrix_type (); | |
179 SparseComplexMatrix ret = v1.sparse_complex_matrix_value () * v2.complex_diag_matrix_value (); | |
180 octave_value out = octave_value (ret); | |
181 typ.mark_as_unsymmetric (); | |
182 out.matrix_type (typ); | |
183 return out; | |
184 } | |
185 } | |
186 | |
187 void | |
188 install_dm_scm_ops (void) | |
189 { | |
190 INSTALL_BINOP (op_mul, octave_diag_matrix, octave_sparse_complex_matrix, | |
191 mul_dm_scm); | |
192 INSTALL_BINOP (op_mul, octave_complex_diag_matrix, octave_sparse_matrix, | |
193 mul_cdm_sm); | |
194 INSTALL_BINOP (op_mul, octave_complex_diag_matrix, octave_sparse_complex_matrix, | |
195 mul_cdm_scm); | |
196 INSTALL_BINOP (op_mul, octave_sparse_complex_matrix, octave_diag_matrix, | |
197 mul_scm_dm); | |
198 INSTALL_BINOP (op_mul, octave_sparse_matrix, octave_complex_diag_matrix, | |
199 mul_sm_cdm); | |
200 INSTALL_BINOP (op_mul, octave_sparse_complex_matrix, octave_complex_diag_matrix, | |
201 mul_scm_cdm); | |
202 } |