Mercurial > octave

comparison liboctave/operators/Sparse-op-defs.h @ 22197:e43d83253e28

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
refill multi-line macro definitions Use the Emacs C++ mode style for line continuation markers in multi-line macro definitions. * make_int.cc, __dsearchn__.cc, __magick_read__.cc, besselj.cc, bitfcns.cc, bsxfun.cc, cellfun.cc, data.cc, defun-dld.h, defun-int.h, defun.h, det.cc, error.h, find.cc, gcd.cc, graphics.cc, interpreter.h, jit-ir.h, jit-typeinfo.h, lookup.cc, ls-mat5.cc, max.cc, mexproto.h, mxarray.in.h, oct-stream.cc, ordschur.cc, pr-output.cc, profiler.h, psi.cc, regexp.cc, sparse-xdiv.cc, sparse-xpow.cc, tril.cc, txt-eng.h, utils.cc, variables.cc, variables.h, xdiv.cc, xpow.cc, __glpk__.cc, ov-base.cc, ov-base.h, ov-cell.cc, ov-ch-mat.cc, ov-classdef.cc, ov-complex.cc, ov-cx-mat.cc, ov-cx-sparse.cc, ov-float.cc, ov-float.h, ov-flt-complex.cc, ov-flt-cx-mat.cc, ov-flt-re-mat.cc, ov-int-traits.h, ov-lazy-idx.h, ov-perm.cc, ov-re-mat.cc, ov-re-sparse.cc, ov-scalar.cc, ov-scalar.h, ov-str-mat.cc, ov-type-conv.h, ov.cc, ov.h, op-class.cc, op-int-conv.cc, op-int.h, op-str-str.cc, ops.h, lex.ll, Array.cc, CMatrix.cc, CSparse.cc, MArray.cc, MArray.h, MDiagArray2.cc, MDiagArray2.h, MSparse.h, Sparse.cc, dMatrix.cc, dSparse.cc, fCMatrix.cc, fMatrix.cc, idx-vector.cc, f77-fcn.h, quit.h, bsxfun-decl.h, bsxfun-defs.cc, lo-specfun.cc, oct-convn.cc, oct-convn.h, oct-norm.cc, oct-norm.h, oct-rand.cc, Sparse-op-decls.h, Sparse-op-defs.h, mx-inlines.cc, mx-op-decl.h, mx-op-defs.h, mach-info.cc, oct-group.cc, oct-passwd.cc, oct-syscalls.cc, oct-time.cc, data-conv.cc, kpse.cc, lo-ieee.h, lo-macros.h, oct-cmplx.h, oct-glob.cc, oct-inttypes.cc, oct-inttypes.h, oct-locbuf.h, oct-sparse.h, url-transfer.cc, oct-conf-post.in.h, shared-fcns.h: Refill macro definitions.
author John W. Eaton <jwe@octave.org>
date Mon, 01 Aug 2016 12:40:18 -0400
parents 278fc29b69ca
children bac0d6f07a3e
comparison
equal deleted inserted replaced
22196:dd992fd74fce 22197:e43d83253e28
31 #include "oct-locbuf.h" 31 #include "oct-locbuf.h"
32 #include "mx-inlines.cc" 32 #include "mx-inlines.cc"
33 33
34 // sparse matrix by scalar operations. 34 // sparse matrix by scalar operations.
35 35
36 #define SPARSE_SMS_BIN_OP_1(R, F, OP, M, S) \ 36 #define SPARSE_SMS_BIN_OP_1(R, F, OP, M, S) \
37 R \ 37 R \
38 F (const M& m, const S& s) \ 38 F (const M& m, const S& s) \
39 { \ 39 { \
40 octave_idx_type nr = m.rows (); \ 40 octave_idx_type nr = m.rows (); \
41 octave_idx_type nc = m.cols (); \ 41 octave_idx_type nc = m.cols (); \
42 \ 42 \
43 R r (nr, nc, (0.0 OP s)); \ 43 R r (nr, nc, (0.0 OP s)); \
44 \ 44 \
45 for (octave_idx_type j = 0; j < nc; j++) \ 45 for (octave_idx_type j = 0; j < nc; j++) \
46 for (octave_idx_type i = m.cidx (j); i < m.cidx (j+1); i++) \ 46 for (octave_idx_type i = m.cidx (j); i < m.cidx (j+1); i++) \
47 r.xelem (m.ridx (i), j) = m.data (i) OP s; \ 47 r.xelem (m.ridx (i), j) = m.data (i) OP s; \
48 return r; \ 48 return r; \
49 } 49 }
50 50
51 #define SPARSE_SMS_BIN_OP_2(R, F, OP, M, S) \ 51 #define SPARSE_SMS_BIN_OP_2(R, F, OP, M, S) \
52 R \ 52 R \
53 F (const M& m, const S& s) \ 53 F (const M& m, const S& s) \
54 { \ 54 { \
55 octave_idx_type nr = m.rows (); \ 55 octave_idx_type nr = m.rows (); \
56 octave_idx_type nc = m.cols (); \ 56 octave_idx_type nc = m.cols (); \
57 octave_idx_type nz = m.nnz (); \ 57 octave_idx_type nz = m.nnz (); \
58 \ 58 \
59 R r (nr, nc, nz); \ 59 R r (nr, nc, nz); \
60 \ 60 \
61 for (octave_idx_type i = 0; i < nz; i++) \ 61 for (octave_idx_type i = 0; i < nz; i++) \
62 { \ 62 { \
63 r.xdata (i) = m.data (i) OP s; \ 63 r.xdata (i) = m.data (i) OP s; \
64 r.xridx (i) = m.ridx (i); \ 64 r.xridx (i) = m.ridx (i); \
65 } \ 65 } \
66 for (octave_idx_type i = 0; i < nc + 1; i++) \ 66 for (octave_idx_type i = 0; i < nc + 1; i++) \
67 r.xcidx (i) = m.cidx (i); \ 67 r.xcidx (i) = m.cidx (i); \
68 \ 68 \
69 r.maybe_compress (true); \ 69 r.maybe_compress (true); \
70 return r; \ 70 return r; \
71 } 71 }
72 72
73 #define SPARSE_SMS_BIN_OPS(R1, R2, M, S) \ 73 #define SPARSE_SMS_BIN_OPS(R1, R2, M, S) \
74 SPARSE_SMS_BIN_OP_1 (R1, operator +, +, M, S) \ 74 SPARSE_SMS_BIN_OP_1 (R1, operator +, +, M, S) \
75 SPARSE_SMS_BIN_OP_1 (R1, operator -, -, M, S) \ 75 SPARSE_SMS_BIN_OP_1 (R1, operator -, -, M, S) \
76 SPARSE_SMS_BIN_OP_2 (R2, operator *, *, M, S) \ 76 SPARSE_SMS_BIN_OP_2 (R2, operator *, *, M, S) \
77 SPARSE_SMS_BIN_OP_2 (R2, operator /, /, M, S) 77 SPARSE_SMS_BIN_OP_2 (R2, operator /, /, M, S)
78 78
79 #define SPARSE_SMS_CMP_OP(F, OP, M, MZ, MC, S, SZ, SC) \ 79 #define SPARSE_SMS_CMP_OP(F, OP, M, MZ, MC, S, SZ, SC) \
80 SparseBoolMatrix \ 80 SparseBoolMatrix \
81 F (const M& m, const S& s) \ 81 F (const M& m, const S& s) \
82 { \ 82 { \
83 octave_idx_type nr = m.rows (); \ 83 octave_idx_type nr = m.rows (); \
84 octave_idx_type nc = m.cols (); \ 84 octave_idx_type nc = m.cols (); \
85 SparseBoolMatrix r; \ 85 SparseBoolMatrix r; \
86 \ 86 \
87 if (MC (MZ) OP SC (s)) \ 87 if (MC (MZ) OP SC (s)) \
88 { \ 88 { \
89 r = SparseBoolMatrix (nr, nc, true); \ 89 r = SparseBoolMatrix (nr, nc, true); \
90 for (octave_idx_type j = 0; j < nc; j++) \ 90 for (octave_idx_type j = 0; j < nc; j++) \
91 for (octave_idx_type i = m.cidx (j); i < m.cidx (j+1); i++) \ 91 for (octave_idx_type i = m.cidx (j); i < m.cidx (j+1); i++) \
92 if (! (MC (m.data (i)) OP SC (s))) \ 92 if (! (MC (m.data (i)) OP SC (s))) \
93 r.data (m.ridx (i) + j * nr) = false; \ 93 r.data (m.ridx (i) + j * nr) = false; \
94 r.maybe_compress (true); \ 94 r.maybe_compress (true); \
95 } \ 95 } \
96 else \ 96 else \
97 { \ 97 { \
98 r = SparseBoolMatrix (nr, nc, m.nnz ()); \ 98 r = SparseBoolMatrix (nr, nc, m.nnz ()); \
99 r.cidx (0) = static_cast<octave_idx_type> (0); \ 99 r.cidx (0) = static_cast<octave_idx_type> (0); \
100 octave_idx_type nel = 0; \ 100 octave_idx_type nel = 0; \
101 for (octave_idx_type j = 0; j < nc; j++) \ 101 for (octave_idx_type j = 0; j < nc; j++) \
102 { \ 102 { \
103 for (octave_idx_type i = m.cidx (j); i < m.cidx (j+1); i++) \ 103 for (octave_idx_type i = m.cidx (j); i < m.cidx (j+1); i++) \
104 if (MC (m.data (i)) OP SC (s)) \ 104 if (MC (m.data (i)) OP SC (s)) \
105 { \ 105 { \
106 r.ridx (nel) = m.ridx (i); \ 106 r.ridx (nel) = m.ridx (i); \
107 r.data (nel++) = true; \ 107 r.data (nel++) = true; \
108 } \ 108 } \
109 r.cidx (j + 1) = nel; \ 109 r.cidx (j + 1) = nel; \
110 } \ 110 } \
111 r.maybe_compress (false); \ 111 r.maybe_compress (false); \
112 } \ 112 } \
113 return r; \ 113 return r; \
114 } 114 }
115 115
116 #define SPARSE_SMS_CMP_OPS(M, MZ, CM, S, SZ, CS) \ 116 #define SPARSE_SMS_CMP_OPS(M, MZ, CM, S, SZ, CS) \
117 SPARSE_SMS_CMP_OP (mx_el_lt, <, M, MZ, , S, SZ, ) \ 117 SPARSE_SMS_CMP_OP (mx_el_lt, <, M, MZ, , S, SZ, ) \
118 SPARSE_SMS_CMP_OP (mx_el_le, <=, M, MZ, , S, SZ, ) \ 118 SPARSE_SMS_CMP_OP (mx_el_le, <=, M, MZ, , S, SZ, ) \
119 SPARSE_SMS_CMP_OP (mx_el_ge, >=, M, MZ, , S, SZ, ) \ 119 SPARSE_SMS_CMP_OP (mx_el_ge, >=, M, MZ, , S, SZ, ) \
120 SPARSE_SMS_CMP_OP (mx_el_gt, >, M, MZ, , S, SZ, ) \ 120 SPARSE_SMS_CMP_OP (mx_el_gt, >, M, MZ, , S, SZ, ) \
121 SPARSE_SMS_CMP_OP (mx_el_eq, ==, M, MZ, , S, SZ, ) \ 121 SPARSE_SMS_CMP_OP (mx_el_eq, ==, M, MZ, , S, SZ, ) \
122 SPARSE_SMS_CMP_OP (mx_el_ne, !=, M, MZ, , S, SZ, ) 122 SPARSE_SMS_CMP_OP (mx_el_ne, !=, M, MZ, , S, SZ, )
123 123
124 #define SPARSE_SMS_EQNE_OPS(M, MZ, CM, S, SZ, CS) \ 124 #define SPARSE_SMS_EQNE_OPS(M, MZ, CM, S, SZ, CS) \
125 SPARSE_SMS_CMP_OP (mx_el_eq, ==, M, MZ, , S, SZ, ) \ 125 SPARSE_SMS_CMP_OP (mx_el_eq, ==, M, MZ, , S, SZ, ) \
126 SPARSE_SMS_CMP_OP (mx_el_ne, !=, M, MZ, , S, SZ, ) 126 SPARSE_SMS_CMP_OP (mx_el_ne, !=, M, MZ, , S, SZ, )
127 127
128 #define SPARSE_SMS_BOOL_OP(F, OP, M, S, LHS_ZERO, RHS_ZERO) \ 128 #define SPARSE_SMS_BOOL_OP(F, OP, M, S, LHS_ZERO, RHS_ZERO) \
129 SparseBoolMatrix \ 129 SparseBoolMatrix \
130 F (const M& m, const S& s) \ 130 F (const M& m, const S& s) \
131 { \ 131 { \
132 octave_idx_type nr = m.rows (); \ 132 octave_idx_type nr = m.rows (); \
133 octave_idx_type nc = m.cols (); \ 133 octave_idx_type nc = m.cols (); \
134 SparseBoolMatrix r; \ 134 SparseBoolMatrix r; \
135 \ 135 \
136 if (nr > 0 && nc > 0) \ 136 if (nr > 0 && nc > 0) \
137 { \ 137 { \
138 if (LHS_ZERO OP (s != RHS_ZERO)) \ 138 if (LHS_ZERO OP (s != RHS_ZERO)) \
139 { \ 139 { \
140 r = SparseBoolMatrix (nr, nc, true); \ 140 r = SparseBoolMatrix (nr, nc, true); \
141 for (octave_idx_type j = 0; j < nc; j++) \ 141 for (octave_idx_type j = 0; j < nc; j++) \
142 for (octave_idx_type i = m.cidx (j); i < m.cidx (j+1); i++) \ 142 for (octave_idx_type i = m.cidx (j); i < m.cidx (j+1); i++) \
143 if (! ((m.data (i) != LHS_ZERO) OP (s != RHS_ZERO))) \ 143 if (! ((m.data (i) != LHS_ZERO) OP (s != RHS_ZERO))) \
144 r.data (m.ridx (i) + j * nr) = false; \ 144 r.data (m.ridx (i) + j * nr) = false; \
145 r.maybe_compress (true); \ 145 r.maybe_compress (true); \
146 } \ 146 } \
147 else \ 147 else \
148 { \ 148 { \
149 r = SparseBoolMatrix (nr, nc, m.nnz ()); \ 149 r = SparseBoolMatrix (nr, nc, m.nnz ()); \
150 r.cidx (0) = static_cast<octave_idx_type> (0); \ 150 r.cidx (0) = static_cast<octave_idx_type> (0); \
151 octave_idx_type nel = 0; \ 151 octave_idx_type nel = 0; \
152 for (octave_idx_type j = 0; j < nc; j++) \ 152 for (octave_idx_type j = 0; j < nc; j++) \
153 { \ 153 { \
154 for (octave_idx_type i = m.cidx (j); i < m.cidx (j+1); i++) \ 154 for (octave_idx_type i = m.cidx (j); i < m.cidx (j+1); i++) \
155 if ((m.data (i) != LHS_ZERO) OP (s != RHS_ZERO)) \ 155 if ((m.data (i) != LHS_ZERO) OP (s != RHS_ZERO)) \
156 { \ 156 { \
157 r.ridx (nel) = m.ridx (i); \ 157 r.ridx (nel) = m.ridx (i); \
158 r.data (nel++) = true; \ 158 r.data (nel++) = true; \
159 } \ 159 } \
160 r.cidx (j + 1) = nel; \ 160 r.cidx (j + 1) = nel; \
161 } \ 161 } \
162 r.maybe_compress (false); \ 162 r.maybe_compress (false); \
163 } \ 163 } \
164 } \ 164 } \
165 return r; \ 165 return r; \
166 } 166 }
167 167
168 #define SPARSE_SMS_BOOL_OPS2(M, S, LHS_ZERO, RHS_ZERO) \ 168 #define SPARSE_SMS_BOOL_OPS2(M, S, LHS_ZERO, RHS_ZERO) \
169 SPARSE_SMS_BOOL_OP (mx_el_and, &&, M, S, LHS_ZERO, RHS_ZERO) \ 169 SPARSE_SMS_BOOL_OP (mx_el_and, &&, M, S, LHS_ZERO, RHS_ZERO) \
170 SPARSE_SMS_BOOL_OP (mx_el_or, ||, M, S, LHS_ZERO, RHS_ZERO) 170 SPARSE_SMS_BOOL_OP (mx_el_or, ||, M, S, LHS_ZERO, RHS_ZERO)
171 171
172 #define SPARSE_SMS_BOOL_OPS(M, S, ZERO) \ 172 #define SPARSE_SMS_BOOL_OPS(M, S, ZERO) \
173 SPARSE_SMS_BOOL_OPS2(M, S, ZERO, ZERO) 173 SPARSE_SMS_BOOL_OPS2(M, S, ZERO, ZERO)
174 174
175 // scalar by sparse matrix operations. 175 // scalar by sparse matrix operations.
176 176
177 #define SPARSE_SSM_BIN_OP_1(R, F, OP, S, M) \ 177 #define SPARSE_SSM_BIN_OP_1(R, F, OP, S, M) \
178 R \ 178 R \
179 F (const S& s, const M& m) \ 179 F (const S& s, const M& m) \
180 { \ 180 { \
181 octave_idx_type nr = m.rows (); \ 181 octave_idx_type nr = m.rows (); \
182 octave_idx_type nc = m.cols (); \ 182 octave_idx_type nc = m.cols (); \
183 \ 183 \
184 R r (nr, nc, (s OP 0.0)); \ 184 R r (nr, nc, (s OP 0.0)); \
185 \ 185 \
186 for (octave_idx_type j = 0; j < nc; j++) \ 186 for (octave_idx_type j = 0; j < nc; j++) \
187 for (octave_idx_type i = m.cidx (j); i < m.cidx (j+1); i++) \ 187 for (octave_idx_type i = m.cidx (j); i < m.cidx (j+1); i++) \
188 r.xelem (m.ridx (i), j) = s OP m.data (i); \ 188 r.xelem (m.ridx (i), j) = s OP m.data (i); \
189 \ 189 \
190 return r; \ 190 return r; \
191 } 191 }
192 192
193 #define SPARSE_SSM_BIN_OP_2(R, F, OP, S, M) \ 193 #define SPARSE_SSM_BIN_OP_2(R, F, OP, S, M) \
194 R \ 194 R \
195 F (const S& s, const M& m) \ 195 F (const S& s, const M& m) \
196 { \ 196 { \
197 octave_idx_type nr = m.rows (); \ 197 octave_idx_type nr = m.rows (); \
198 octave_idx_type nc = m.cols (); \ 198 octave_idx_type nc = m.cols (); \
199 octave_idx_type nz = m.nnz (); \ 199 octave_idx_type nz = m.nnz (); \
200 \ 200 \
201 R r (nr, nc, nz); \ 201 R r (nr, nc, nz); \
202 \ 202 \
203 for (octave_idx_type i = 0; i < nz; i++) \ 203 for (octave_idx_type i = 0; i < nz; i++) \
204 { \ 204 { \
205 r.xdata (i) = s OP m.data (i); \ 205 r.xdata (i) = s OP m.data (i); \
206 r.xridx (i) = m.ridx (i); \ 206 r.xridx (i) = m.ridx (i); \
207 } \ 207 } \
208 for (octave_idx_type i = 0; i < nc + 1; i++) \ 208 for (octave_idx_type i = 0; i < nc + 1; i++) \
209 r.xcidx (i) = m.cidx (i); \ 209 r.xcidx (i) = m.cidx (i); \
210 \ 210 \
211 r.maybe_compress(true); \ 211 r.maybe_compress(true); \
212 return r; \ 212 return r; \
213 } 213 }
214 214
215 #define SPARSE_SSM_BIN_OPS(R1, R2, S, M) \ 215 #define SPARSE_SSM_BIN_OPS(R1, R2, S, M) \
216 SPARSE_SSM_BIN_OP_1 (R1, operator +, +, S, M) \ 216 SPARSE_SSM_BIN_OP_1 (R1, operator +, +, S, M) \
217 SPARSE_SSM_BIN_OP_1 (R1, operator -, -, S, M) \ 217 SPARSE_SSM_BIN_OP_1 (R1, operator -, -, S, M) \
218 SPARSE_SSM_BIN_OP_2 (R2, operator *, *, S, M) \ 218 SPARSE_SSM_BIN_OP_2 (R2, operator *, *, S, M) \
219 SPARSE_SSM_BIN_OP_2 (R2, operator /, /, S, M) 219 SPARSE_SSM_BIN_OP_2 (R2, operator /, /, S, M)
220 220
221 #define SPARSE_SSM_CMP_OP(F, OP, S, SZ, SC, M, MZ, MC) \ 221 #define SPARSE_SSM_CMP_OP(F, OP, S, SZ, SC, M, MZ, MC) \
222 SparseBoolMatrix \ 222 SparseBoolMatrix \
223 F (const S& s, const M& m) \ 223 F (const S& s, const M& m) \
224 { \ 224 { \
225 octave_idx_type nr = m.rows (); \ 225 octave_idx_type nr = m.rows (); \
226 octave_idx_type nc = m.cols (); \ 226 octave_idx_type nc = m.cols (); \
227 SparseBoolMatrix r; \ 227 SparseBoolMatrix r; \
228 \ 228 \
229 if (SC (s) OP SC (MZ)) \ 229 if (SC (s) OP SC (MZ)) \
230 { \ 230 { \
231 r = SparseBoolMatrix (nr, nc, true); \ 231 r = SparseBoolMatrix (nr, nc, true); \
232 for (octave_idx_type j = 0; j < nc; j++) \ 232 for (octave_idx_type j = 0; j < nc; j++) \
233 for (octave_idx_type i = m.cidx (j); i < m.cidx (j+1); i++) \ 233 for (octave_idx_type i = m.cidx (j); i < m.cidx (j+1); i++) \
234 if (! (SC (s) OP MC (m.data (i)))) \ 234 if (! (SC (s) OP MC (m.data (i)))) \
235 r.data (m.ridx (i) + j * nr) = false; \ 235 r.data (m.ridx (i) + j * nr) = false; \
236 r.maybe_compress (true); \ 236 r.maybe_compress (true); \
237 } \ 237 } \
238 else \ 238 else \
239 { \ 239 { \
240 r = SparseBoolMatrix (nr, nc, m.nnz ()); \ 240 r = SparseBoolMatrix (nr, nc, m.nnz ()); \
241 r.cidx (0) = static_cast<octave_idx_type> (0); \ 241 r.cidx (0) = static_cast<octave_idx_type> (0); \
242 octave_idx_type nel = 0; \ 242 octave_idx_type nel = 0; \
243 for (octave_idx_type j = 0; j < nc; j++) \ 243 for (octave_idx_type j = 0; j < nc; j++) \
244 { \ 244 { \
245 for (octave_idx_type i = m.cidx (j); i < m.cidx (j+1); i++) \ 245 for (octave_idx_type i = m.cidx (j); i < m.cidx (j+1); i++) \
246 if (SC (s) OP MC (m.data (i))) \ 246 if (SC (s) OP MC (m.data (i))) \
247 { \ 247 { \
248 r.ridx (nel) = m.ridx (i); \ 248 r.ridx (nel) = m.ridx (i); \
249 r.data (nel++) = true; \ 249 r.data (nel++) = true; \
250 } \ 250 } \
251 r.cidx (j + 1) = nel; \ 251 r.cidx (j + 1) = nel; \
252 } \ 252 } \
253 r.maybe_compress (false); \ 253 r.maybe_compress (false); \
254 } \ 254 } \
255 return r; \ 255 return r; \
256 } 256 }
257 257
258 #define SPARSE_SSM_CMP_OPS(S, SZ, SC, M, MZ, MC) \ 258 #define SPARSE_SSM_CMP_OPS(S, SZ, SC, M, MZ, MC) \
259 SPARSE_SSM_CMP_OP (mx_el_lt, <, S, SZ, , M, MZ, ) \ 259 SPARSE_SSM_CMP_OP (mx_el_lt, <, S, SZ, , M, MZ, ) \
260 SPARSE_SSM_CMP_OP (mx_el_le, <=, S, SZ, , M, MZ, ) \ 260 SPARSE_SSM_CMP_OP (mx_el_le, <=, S, SZ, , M, MZ, ) \
261 SPARSE_SSM_CMP_OP (mx_el_ge, >=, S, SZ, , M, MZ, ) \ 261 SPARSE_SSM_CMP_OP (mx_el_ge, >=, S, SZ, , M, MZ, ) \
262 SPARSE_SSM_CMP_OP (mx_el_gt, >, S, SZ, , M, MZ, ) \ 262 SPARSE_SSM_CMP_OP (mx_el_gt, >, S, SZ, , M, MZ, ) \
263 SPARSE_SSM_CMP_OP (mx_el_eq, ==, S, SZ, , M, MZ, ) \ 263 SPARSE_SSM_CMP_OP (mx_el_eq, ==, S, SZ, , M, MZ, ) \
264 SPARSE_SSM_CMP_OP (mx_el_ne, !=, S, SZ, , M, MZ, ) 264 SPARSE_SSM_CMP_OP (mx_el_ne, !=, S, SZ, , M, MZ, )
265 265
266 #define SPARSE_SSM_EQNE_OPS(S, SZ, SC, M, MZ, MC) \ 266 #define SPARSE_SSM_EQNE_OPS(S, SZ, SC, M, MZ, MC) \
267 SPARSE_SSM_CMP_OP (mx_el_eq, ==, S, SZ, , M, MZ, ) \ 267 SPARSE_SSM_CMP_OP (mx_el_eq, ==, S, SZ, , M, MZ, ) \
268 SPARSE_SSM_CMP_OP (mx_el_ne, !=, S, SZ, , M, MZ, ) 268 SPARSE_SSM_CMP_OP (mx_el_ne, !=, S, SZ, , M, MZ, )
269 269
270 #define SPARSE_SSM_BOOL_OP(F, OP, S, M, LHS_ZERO, RHS_ZERO) \ 270 #define SPARSE_SSM_BOOL_OP(F, OP, S, M, LHS_ZERO, RHS_ZERO) \
271 SparseBoolMatrix \ 271 SparseBoolMatrix \
272 F (const S& s, const M& m) \ 272 F (const S& s, const M& m) \
273 { \ 273 { \
274 octave_idx_type nr = m.rows (); \ 274 octave_idx_type nr = m.rows (); \
275 octave_idx_type nc = m.cols (); \ 275 octave_idx_type nc = m.cols (); \
276 SparseBoolMatrix r; \ 276 SparseBoolMatrix r; \
277 \ 277 \
278 if (nr > 0 && nc > 0) \ 278 if (nr > 0 && nc > 0) \
279 { \ 279 { \
280 if ((s != LHS_ZERO) OP RHS_ZERO) \ 280 if ((s != LHS_ZERO) OP RHS_ZERO) \
281 { \ 281 { \
282 r = SparseBoolMatrix (nr, nc, true); \ 282 r = SparseBoolMatrix (nr, nc, true); \
283 for (octave_idx_type j = 0; j < nc; j++) \ 283 for (octave_idx_type j = 0; j < nc; j++) \
284 for (octave_idx_type i = m.cidx (j); i < m.cidx (j+1); i++) \ 284 for (octave_idx_type i = m.cidx (j); i < m.cidx (j+1); i++) \
285 if (! ((s != LHS_ZERO) OP (m.data (i) != RHS_ZERO))) \ 285 if (! ((s != LHS_ZERO) OP (m.data (i) != RHS_ZERO))) \
286 r.data (m.ridx (i) + j * nr) = false; \ 286 r.data (m.ridx (i) + j * nr) = false; \
287 r.maybe_compress (true); \ 287 r.maybe_compress (true); \
288 } \ 288 } \
289 else \ 289 else \
290 { \ 290 { \
291 r = SparseBoolMatrix (nr, nc, m.nnz ()); \ 291 r = SparseBoolMatrix (nr, nc, m.nnz ()); \
292 r.cidx (0) = static_cast<octave_idx_type> (0); \ 292 r.cidx (0) = static_cast<octave_idx_type> (0); \
293 octave_idx_type nel = 0; \ 293 octave_idx_type nel = 0; \
294 for (octave_idx_type j = 0; j < nc; j++) \ 294 for (octave_idx_type j = 0; j < nc; j++) \
295 { \ 295 { \
296 for (octave_idx_type i = m.cidx (j); i < m.cidx (j+1); i++) \ 296 for (octave_idx_type i = m.cidx (j); i < m.cidx (j+1); i++) \
297 if ((s != LHS_ZERO) OP (m.data (i) != RHS_ZERO)) \ 297 if ((s != LHS_ZERO) OP (m.data (i) != RHS_ZERO)) \
298 { \ 298 { \
299 r.ridx (nel) = m.ridx (i); \ 299 r.ridx (nel) = m.ridx (i); \
300 r.data (nel++) = true; \ 300 r.data (nel++) = true; \
301 } \ 301 } \
302 r.cidx (j + 1) = nel; \ 302 r.cidx (j + 1) = nel; \
303 } \ 303 } \
304 r.maybe_compress (false); \ 304 r.maybe_compress (false); \
305 } \ 305 } \
306 } \ 306 } \
307 return r; \ 307 return r; \
308 } 308 }
309 309
310 #define SPARSE_SSM_BOOL_OPS2(S, M, LHS_ZERO, RHS_ZERO) \ 310 #define SPARSE_SSM_BOOL_OPS2(S, M, LHS_ZERO, RHS_ZERO) \
311 SPARSE_SSM_BOOL_OP (mx_el_and, &&, S, M, LHS_ZERO, RHS_ZERO) \ 311 SPARSE_SSM_BOOL_OP (mx_el_and, &&, S, M, LHS_ZERO, RHS_ZERO) \
312 SPARSE_SSM_BOOL_OP (mx_el_or, ||, S, M, LHS_ZERO, RHS_ZERO) 312 SPARSE_SSM_BOOL_OP (mx_el_or, ||, S, M, LHS_ZERO, RHS_ZERO)
313 313
314 #define SPARSE_SSM_BOOL_OPS(S, M, ZERO) \ 314 #define SPARSE_SSM_BOOL_OPS(S, M, ZERO) \
315 SPARSE_SSM_BOOL_OPS2(S, M, ZERO, ZERO) 315 SPARSE_SSM_BOOL_OPS2(S, M, ZERO, ZERO)
316 316
317 // sparse matrix by sparse matrix operations. 317 // sparse matrix by sparse matrix operations.
318 318
319 #define SPARSE_SMSM_BIN_OP_1(R, F, OP, M1, M2) \ 319 #define SPARSE_SMSM_BIN_OP_1(R, F, OP, M1, M2) \
320 R \ 320 R \
321 F (const M1& m1, const M2& m2) \ 321 F (const M1& m1, const M2& m2) \
322 { \ 322 { \
323 R r; \ 323 R r; \
324 \ 324 \
325 octave_idx_type m1_nr = m1.rows (); \ 325 octave_idx_type m1_nr = m1.rows (); \
326 octave_idx_type m1_nc = m1.cols (); \ 326 octave_idx_type m1_nc = m1.cols (); \
327 \ 327 \
328 octave_idx_type m2_nr = m2.rows (); \ 328 octave_idx_type m2_nr = m2.rows (); \
329 octave_idx_type m2_nc = m2.cols (); \ 329 octave_idx_type m2_nc = m2.cols (); \
330 \ 330 \
331 if (m1_nr == 1 && m1_nc == 1) \ 331 if (m1_nr == 1 && m1_nc == 1) \
332 { \ 332 { \
333 if (m1.elem (0,0) == 0.) \ 333 if (m1.elem (0,0) == 0.) \
334 r = OP R (m2); \ 334 r = OP R (m2); \
335 else \ 335 else \
336 { \ 336 { \
337 r = R (m2_nr, m2_nc, m1.data (0) OP 0.); \ 337 r = R (m2_nr, m2_nc, m1.data (0) OP 0.); \
338 \ 338 \
339 for (octave_idx_type j = 0 ; j < m2_nc ; j++) \ 339 for (octave_idx_type j = 0 ; j < m2_nc ; j++) \
340 { \ 340 { \
341 octave_quit (); \ 341 octave_quit (); \
342 octave_idx_type idxj = j * m2_nr; \ 342 octave_idx_type idxj = j * m2_nr; \
343 for (octave_idx_type i = m2.cidx (j) ; i < m2.cidx (j+1) ; i++) \ 343 for (octave_idx_type i = m2.cidx (j) ; i < m2.cidx (j+1) ; i++) \
344 { \ 344 { \
345 octave_quit (); \ 345 octave_quit (); \
346 r.data (idxj + m2.ridx (i)) = m1.data (0) OP m2.data (i); \ 346 r.data (idxj + m2.ridx (i)) = m1.data (0) OP m2.data (i); \
347 } \ 347 } \
348 } \ 348 } \
349 r.maybe_compress (); \ 349 r.maybe_compress (); \
350 } \ 350 } \
351 } \ 351 } \
352 else if (m2_nr == 1 && m2_nc == 1) \ 352 else if (m2_nr == 1 && m2_nc == 1) \
353 { \ 353 { \
354 if (m2.elem (0,0) == 0.) \ 354 if (m2.elem (0,0) == 0.) \
355 r = R (m1); \ 355 r = R (m1); \
356 else \ 356 else \
357 { \ 357 { \
358 r = R (m1_nr, m1_nc, 0. OP m2.data (0)); \ 358 r = R (m1_nr, m1_nc, 0. OP m2.data (0)); \
359 \ 359 \
360 for (octave_idx_type j = 0 ; j < m1_nc ; j++) \ 360 for (octave_idx_type j = 0 ; j < m1_nc ; j++) \
361 { \ 361 { \
362 octave_quit (); \ 362 octave_quit (); \
363 octave_idx_type idxj = j * m1_nr; \ 363 octave_idx_type idxj = j * m1_nr; \
364 for (octave_idx_type i = m1.cidx (j) ; i < m1.cidx (j+1) ; i++) \ 364 for (octave_idx_type i = m1.cidx (j) ; i < m1.cidx (j+1) ; i++) \
365 { \ 365 { \
366 octave_quit (); \ 366 octave_quit (); \
367 r.data (idxj + m1.ridx (i)) = m1.data (i) OP m2.data (0); \ 367 r.data (idxj + m1.ridx (i)) = m1.data (i) OP m2.data (0); \
368 } \ 368 } \
369 } \ 369 } \
370 r.maybe_compress (); \ 370 r.maybe_compress (); \
371 } \ 371 } \
372 } \ 372 } \
373 else if (m1_nr != m2_nr || m1_nc != m2_nc) \ 373 else if (m1_nr != m2_nr || m1_nc != m2_nc) \
374 err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \ 374 err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \
375 else \ 375 else \
376 { \ 376 { \
377 r = R (m1_nr, m1_nc, (m1.nnz () + m2.nnz ())); \ 377 r = R (m1_nr, m1_nc, (m1.nnz () + m2.nnz ())); \
378 \ 378 \
379 octave_idx_type jx = 0; \ 379 octave_idx_type jx = 0; \
380 r.cidx (0) = 0; \ 380 r.cidx (0) = 0; \
381 for (octave_idx_type i = 0 ; i < m1_nc ; i++) \ 381 for (octave_idx_type i = 0 ; i < m1_nc ; i++) \
382 { \ 382 { \
383 octave_idx_type ja = m1.cidx (i); \ 383 octave_idx_type ja = m1.cidx (i); \
384 octave_idx_type ja_max = m1.cidx (i+1); \ 384 octave_idx_type ja_max = m1.cidx (i+1); \
385 bool ja_lt_max = ja < ja_max; \ 385 bool ja_lt_max = ja < ja_max; \
386 \ 386 \
387 octave_idx_type jb = m2.cidx (i); \ 387 octave_idx_type jb = m2.cidx (i); \
388 octave_idx_type jb_max = m2.cidx (i+1); \ 388 octave_idx_type jb_max = m2.cidx (i+1); \
389 bool jb_lt_max = jb < jb_max; \ 389 bool jb_lt_max = jb < jb_max; \
390 \ 390 \
391 while (ja_lt_max || jb_lt_max) \ 391 while (ja_lt_max || jb_lt_max) \
392 { \ 392 { \
393 octave_quit (); \ 393 octave_quit (); \
394 if ((! jb_lt_max) || \ 394 if ((! jb_lt_max) || \
395 (ja_lt_max && (m1.ridx (ja) < m2.ridx (jb)))) \ 395 (ja_lt_max && (m1.ridx (ja) < m2.ridx (jb)))) \
396 { \ 396 { \
397 r.ridx (jx) = m1.ridx (ja); \ 397 r.ridx (jx) = m1.ridx (ja); \
398 r.data (jx) = m1.data (ja) OP 0.; \ 398 r.data (jx) = m1.data (ja) OP 0.; \
399 jx++; \ 399 jx++; \
400 ja++; \ 400 ja++; \
401 ja_lt_max= ja < ja_max; \ 401 ja_lt_max= ja < ja_max; \
402 } \ 402 } \
403 else if ((! ja_lt_max) || \ 403 else if ((! ja_lt_max) || \
404 (jb_lt_max && (m2.ridx (jb) < m1.ridx (ja)))) \ 404 (jb_lt_max && (m2.ridx (jb) < m1.ridx (ja)))) \
405 { \ 405 { \
406 r.ridx (jx) = m2.ridx (jb); \ 406 r.ridx (jx) = m2.ridx (jb); \
407 r.data (jx) = 0. OP m2.data (jb); \ 407 r.data (jx) = 0. OP m2.data (jb); \
408 jx++; \ 408 jx++; \
409 jb++; \ 409 jb++; \
410 jb_lt_max= jb < jb_max; \ 410 jb_lt_max= jb < jb_max; \
411 } \ 411 } \
412 else \ 412 else \
413 { \ 413 { \
414 if ((m1.data (ja) OP m2.data (jb)) != 0.) \ 414 if ((m1.data (ja) OP m2.data (jb)) != 0.) \
415 { \ 415 { \
416 r.data (jx) = m1.data (ja) OP m2.data (jb); \ 416 r.data (jx) = m1.data (ja) OP m2.data (jb); \
417 r.ridx (jx) = m1.ridx (ja); \ 417 r.ridx (jx) = m1.ridx (ja); \
418 jx++; \ 418 jx++; \
419 } \ 419 } \
420 ja++; \ 420 ja++; \
421 ja_lt_max= ja < ja_max; \ 421 ja_lt_max= ja < ja_max; \
422 jb++; \ 422 jb++; \
423 jb_lt_max= jb < jb_max; \ 423 jb_lt_max= jb < jb_max; \
424 } \ 424 } \
425 } \ 425 } \
426 r.cidx (i+1) = jx; \ 426 r.cidx (i+1) = jx; \
427 } \ 427 } \
428 \ 428 \
429 r.maybe_compress (); \ 429 r.maybe_compress (); \
430 } \ 430 } \
431 \ 431 \
432 return r; \ 432 return r; \
433 } 433 }
434 434
435 #define SPARSE_SMSM_BIN_OP_2(R, F, OP, M1, M2) \ 435 #define SPARSE_SMSM_BIN_OP_2(R, F, OP, M1, M2) \
436 R \ 436 R \
437 F (const M1& m1, const M2& m2) \ 437 F (const M1& m1, const M2& m2) \
438 { \ 438 { \
439 R r; \ 439 R r; \
440 \ 440 \
441 octave_idx_type m1_nr = m1.rows (); \ 441 octave_idx_type m1_nr = m1.rows (); \
442 octave_idx_type m1_nc = m1.cols (); \ 442 octave_idx_type m1_nc = m1.cols (); \
443 \ 443 \
444 octave_idx_type m2_nr = m2.rows (); \ 444 octave_idx_type m2_nr = m2.rows (); \
445 octave_idx_type m2_nc = m2.cols (); \ 445 octave_idx_type m2_nc = m2.cols (); \
446 \ 446 \
447 if (m1_nr == 1 && m1_nc == 1) \ 447 if (m1_nr == 1 && m1_nc == 1) \
448 { \ 448 { \
449 if (m1.elem (0,0) == 0.) \ 449 if (m1.elem (0,0) == 0.) \
450 r = R (m2_nr, m2_nc); \ 450 r = R (m2_nr, m2_nc); \
451 else \ 451 else \
452 { \ 452 { \
453 r = R (m2); \ 453 r = R (m2); \
454 octave_idx_type m2_nnz = m2.nnz (); \ 454 octave_idx_type m2_nnz = m2.nnz (); \
455 \ 455 \
456 for (octave_idx_type i = 0 ; i < m2_nnz ; i++) \ 456 for (octave_idx_type i = 0 ; i < m2_nnz ; i++) \
457 { \ 457 { \
458 octave_quit (); \ 458 octave_quit (); \
459 r.data (i) = m1.data (0) OP r.data (i); \ 459 r.data (i) = m1.data (0) OP r.data (i); \
460 } \ 460 } \
461 r.maybe_compress (); \ 461 r.maybe_compress (); \
462 } \ 462 } \
463 } \ 463 } \
464 else if (m2_nr == 1 && m2_nc == 1) \ 464 else if (m2_nr == 1 && m2_nc == 1) \
465 { \ 465 { \
466 if (m2.elem (0,0) == 0.) \ 466 if (m2.elem (0,0) == 0.) \
467 r = R (m1_nr, m1_nc); \ 467 r = R (m1_nr, m1_nc); \
468 else \ 468 else \
469 { \ 469 { \
470 r = R (m1); \ 470 r = R (m1); \
471 octave_idx_type m1_nnz = m1.nnz (); \ 471 octave_idx_type m1_nnz = m1.nnz (); \
472 \ 472 \
473 for (octave_idx_type i = 0 ; i < m1_nnz ; i++) \ 473 for (octave_idx_type i = 0 ; i < m1_nnz ; i++) \
474 { \ 474 { \
475 octave_quit (); \ 475 octave_quit (); \
476 r.data (i) = r.data (i) OP m2.data (0); \ 476 r.data (i) = r.data (i) OP m2.data (0); \
477 } \ 477 } \
478 r.maybe_compress (); \ 478 r.maybe_compress (); \
479 } \ 479 } \
480 } \ 480 } \
481 else if (m1_nr != m2_nr || m1_nc != m2_nc) \ 481 else if (m1_nr != m2_nr || m1_nc != m2_nc) \
482 err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \ 482 err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \
483 else \ 483 else \
484 { \ 484 { \
485 r = R (m1_nr, m1_nc, (m1.nnz () > m2.nnz () ? m1.nnz () : m2.nnz ())); \ 485 r = R (m1_nr, m1_nc, (m1.nnz () > m2.nnz () ? m1.nnz () : m2.nnz ())); \
486 \ 486 \
487 octave_idx_type jx = 0; \ 487 octave_idx_type jx = 0; \
488 r.cidx (0) = 0; \ 488 r.cidx (0) = 0; \
489 for (octave_idx_type i = 0 ; i < m1_nc ; i++) \ 489 for (octave_idx_type i = 0 ; i < m1_nc ; i++) \
490 { \ 490 { \
491 octave_idx_type ja = m1.cidx (i); \ 491 octave_idx_type ja = m1.cidx (i); \
492 octave_idx_type ja_max = m1.cidx (i+1); \ 492 octave_idx_type ja_max = m1.cidx (i+1); \
493 bool ja_lt_max = ja < ja_max; \ 493 bool ja_lt_max = ja < ja_max; \
494 \ 494 \
495 octave_idx_type jb = m2.cidx (i); \ 495 octave_idx_type jb = m2.cidx (i); \
496 octave_idx_type jb_max = m2.cidx (i+1); \ 496 octave_idx_type jb_max = m2.cidx (i+1); \
497 bool jb_lt_max = jb < jb_max; \ 497 bool jb_lt_max = jb < jb_max; \
498 \ 498 \
499 while (ja_lt_max || jb_lt_max) \ 499 while (ja_lt_max || jb_lt_max) \
500 { \ 500 { \
501 octave_quit (); \ 501 octave_quit (); \
502 if ((! jb_lt_max) || \ 502 if ((! jb_lt_max) || \
503 (ja_lt_max && (m1.ridx (ja) < m2.ridx (jb)))) \ 503 (ja_lt_max && (m1.ridx (ja) < m2.ridx (jb)))) \
504 { \ 504 { \
505 ja++; ja_lt_max= ja < ja_max; \ 505 ja++; ja_lt_max= ja < ja_max; \
506 } \ 506 } \
507 else if ((! ja_lt_max) || \ 507 else if ((! ja_lt_max) || \
508 (jb_lt_max && (m2.ridx (jb) < m1.ridx (ja)))) \ 508 (jb_lt_max && (m2.ridx (jb) < m1.ridx (ja)))) \
509 { \ 509 { \
510 jb++; jb_lt_max= jb < jb_max; \ 510 jb++; jb_lt_max= jb < jb_max; \
511 } \ 511 } \
512 else \ 512 else \
513 { \ 513 { \
514 if ((m1.data (ja) OP m2.data (jb)) != 0.) \ 514 if ((m1.data (ja) OP m2.data (jb)) != 0.) \
515 { \ 515 { \
516 r.data (jx) = m1.data (ja) OP m2.data (jb); \ 516 r.data (jx) = m1.data (ja) OP m2.data (jb); \
517 r.ridx (jx) = m1.ridx (ja); \ 517 r.ridx (jx) = m1.ridx (ja); \
518 jx++; \ 518 jx++; \
519 } \ 519 } \
520 ja++; ja_lt_max= ja < ja_max; \ 520 ja++; ja_lt_max= ja < ja_max; \
521 jb++; jb_lt_max= jb < jb_max; \ 521 jb++; jb_lt_max= jb < jb_max; \
522 } \ 522 } \
523 } \ 523 } \
524 r.cidx (i+1) = jx; \ 524 r.cidx (i+1) = jx; \
525 } \ 525 } \
526 \ 526 \
527 r.maybe_compress (); \ 527 r.maybe_compress (); \
528 } \ 528 } \
529 \ 529 \
530 return r; \ 530 return r; \
531 } 531 }
532 532
533 #define SPARSE_SMSM_BIN_OP_3(R, F, OP, M1, M2) \ 533 #define SPARSE_SMSM_BIN_OP_3(R, F, OP, M1, M2) \
534 R \ 534 R \
535 F (const M1& m1, const M2& m2) \ 535 F (const M1& m1, const M2& m2) \
536 { \ 536 { \
537 R r; \ 537 R r; \
538 \ 538 \
539 octave_idx_type m1_nr = m1.rows (); \ 539 octave_idx_type m1_nr = m1.rows (); \
540 octave_idx_type m1_nc = m1.cols (); \ 540 octave_idx_type m1_nc = m1.cols (); \
541 \ 541 \
542 octave_idx_type m2_nr = m2.rows (); \ 542 octave_idx_type m2_nr = m2.rows (); \
543 octave_idx_type m2_nc = m2.cols (); \ 543 octave_idx_type m2_nc = m2.cols (); \
544 \ 544 \
545 if (m1_nr == 1 && m1_nc == 1) \ 545 if (m1_nr == 1 && m1_nc == 1) \
546 { \ 546 { \
547 if ((m1.elem (0,0) OP Complex ()) == Complex ()) \ 547 if ((m1.elem (0,0) OP Complex ()) == Complex ()) \
548 { \ 548 { \
549 octave_idx_type m2_nnz = m2.nnz (); \ 549 octave_idx_type m2_nnz = m2.nnz (); \
550 r = R (m2); \ 550 r = R (m2); \
551 for (octave_idx_type i = 0 ; i < m2_nnz ; i++) \ 551 for (octave_idx_type i = 0 ; i < m2_nnz ; i++) \
552 r.data (i) = m1.elem (0,0) OP r.data (i); \ 552 r.data (i) = m1.elem (0,0) OP r.data (i); \
553 r.maybe_compress (); \ 553 r.maybe_compress (); \
554 } \ 554 } \
555 else \ 555 else \
556 { \ 556 { \
557 r = R (m2_nr, m2_nc, m1.elem (0,0) OP Complex ()); \ 557 r = R (m2_nr, m2_nc, m1.elem (0,0) OP Complex ()); \
558 for (octave_idx_type j = 0 ; j < m2_nc ; j++) \ 558 for (octave_idx_type j = 0 ; j < m2_nc ; j++) \
559 { \ 559 { \
560 octave_quit (); \ 560 octave_quit (); \
561 octave_idx_type idxj = j * m2_nr; \ 561 octave_idx_type idxj = j * m2_nr; \
562 for (octave_idx_type i = m2.cidx (j) ; i < m2.cidx (j+1) ; i++) \ 562 for (octave_idx_type i = m2.cidx (j) ; i < m2.cidx (j+1) ; i++) \
563 { \ 563 { \
564 octave_quit (); \ 564 octave_quit (); \
565 r.data (idxj + m2.ridx (i)) = m1.elem (0,0) OP m2.data (i); \ 565 r.data (idxj + m2.ridx (i)) = m1.elem (0,0) OP m2.data (i); \
566 } \ 566 } \
567 } \ 567 } \
568 r.maybe_compress (); \ 568 r.maybe_compress (); \
569 } \ 569 } \
570 } \ 570 } \
571 else if (m2_nr == 1 && m2_nc == 1) \ 571 else if (m2_nr == 1 && m2_nc == 1) \
572 { \ 572 { \
573 if ((Complex () OP m1.elem (0,0)) == Complex ()) \ 573 if ((Complex () OP m1.elem (0,0)) == Complex ()) \
574 { \ 574 { \
575 octave_idx_type m1_nnz = m1.nnz (); \ 575 octave_idx_type m1_nnz = m1.nnz (); \
576 r = R (m1); \ 576 r = R (m1); \
577 for (octave_idx_type i = 0 ; i < m1_nnz ; i++) \ 577 for (octave_idx_type i = 0 ; i < m1_nnz ; i++) \
578 r.data (i) = r.data (i) OP m2.elem (0,0); \ 578 r.data (i) = r.data (i) OP m2.elem (0,0); \
579 r.maybe_compress (); \ 579 r.maybe_compress (); \
580 } \ 580 } \
581 else \ 581 else \
582 { \ 582 { \
583 r = R (m1_nr, m1_nc, Complex () OP m2.elem (0,0)); \ 583 r = R (m1_nr, m1_nc, Complex () OP m2.elem (0,0)); \
584 for (octave_idx_type j = 0 ; j < m1_nc ; j++) \ 584 for (octave_idx_type j = 0 ; j < m1_nc ; j++) \
585 { \ 585 { \
586 octave_quit (); \ 586 octave_quit (); \
587 octave_idx_type idxj = j * m1_nr; \ 587 octave_idx_type idxj = j * m1_nr; \
588 for (octave_idx_type i = m1.cidx (j) ; i < m1.cidx (j+1) ; i++) \ 588 for (octave_idx_type i = m1.cidx (j) ; i < m1.cidx (j+1) ; i++) \
589 { \ 589 { \
590 octave_quit (); \ 590 octave_quit (); \
591 r.data (idxj + m1.ridx (i)) = m1.data (i) OP m2.elem (0,0); \ 591 r.data (idxj + m1.ridx (i)) = m1.data (i) OP m2.elem (0,0); \
592 } \ 592 } \
593 } \ 593 } \
594 r.maybe_compress (); \ 594 r.maybe_compress (); \
595 } \ 595 } \
596 } \ 596 } \
597 else if (m1_nr != m2_nr || m1_nc != m2_nc) \ 597 else if (m1_nr != m2_nr || m1_nc != m2_nc) \
598 err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \ 598 err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \
599 else \ 599 else \
600 { \ 600 { \
601 \ 601 \
602 /* FIXME: Kludge... Always double/Complex, so Complex () */ \ 602 /* FIXME: Kludge... Always double/Complex, so Complex () */ \
603 r = R (m1_nr, m1_nc, (Complex () OP Complex ())); \ 603 r = R (m1_nr, m1_nc, (Complex () OP Complex ())); \
604 \ 604 \
605 for (octave_idx_type i = 0 ; i < m1_nc ; i++) \ 605 for (octave_idx_type i = 0 ; i < m1_nc ; i++) \
606 { \ 606 { \
607 octave_idx_type ja = m1.cidx (i); \ 607 octave_idx_type ja = m1.cidx (i); \
608 octave_idx_type ja_max = m1.cidx (i+1); \ 608 octave_idx_type ja_max = m1.cidx (i+1); \
609 bool ja_lt_max = ja < ja_max; \ 609 bool ja_lt_max = ja < ja_max; \
610 \ 610 \
611 octave_idx_type jb = m2.cidx (i); \ 611 octave_idx_type jb = m2.cidx (i); \
612 octave_idx_type jb_max = m2.cidx (i+1); \ 612 octave_idx_type jb_max = m2.cidx (i+1); \
613 bool jb_lt_max = jb < jb_max; \ 613 bool jb_lt_max = jb < jb_max; \
614 \ 614 \
615 while (ja_lt_max || jb_lt_max) \ 615 while (ja_lt_max || jb_lt_max) \
616 { \ 616 { \
617 octave_quit (); \ 617 octave_quit (); \
618 if ((! jb_lt_max) || \ 618 if ((! jb_lt_max) || \
619 (ja_lt_max && (m1.ridx (ja) < m2.ridx (jb)))) \ 619 (ja_lt_max && (m1.ridx (ja) < m2.ridx (jb)))) \
620 { \ 620 { \
621 /* keep those kludges coming */ \ 621 /* keep those kludges coming */ \
622 r.elem (m1.ridx (ja),i) = m1.data (ja) OP Complex (); \ 622 r.elem (m1.ridx (ja),i) = m1.data (ja) OP Complex (); \
623 ja++; \ 623 ja++; \
624 ja_lt_max= ja < ja_max; \ 624 ja_lt_max= ja < ja_max; \
625 } \ 625 } \
626 else if ((! ja_lt_max) || \ 626 else if ((! ja_lt_max) || \
627 (jb_lt_max && (m2.ridx (jb) < m1.ridx (ja)))) \ 627 (jb_lt_max && (m2.ridx (jb) < m1.ridx (ja)))) \
628 { \ 628 { \
629 /* keep those kludges coming */ \ 629 /* keep those kludges coming */ \
630 r.elem (m2.ridx (jb),i) = Complex () OP m2.data (jb); \ 630 r.elem (m2.ridx (jb),i) = Complex () OP m2.data (jb); \
631 jb++; \ 631 jb++; \
632 jb_lt_max= jb < jb_max; \ 632 jb_lt_max= jb < jb_max; \
633 } \ 633 } \
634 else \ 634 else \
635 { \ 635 { \
636 r.elem (m1.ridx (ja),i) = m1.data (ja) OP m2.data (jb); \ 636 r.elem (m1.ridx (ja),i) = m1.data (ja) OP m2.data (jb); \
637 ja++; \ 637 ja++; \
638 ja_lt_max= ja < ja_max; \ 638 ja_lt_max= ja < ja_max; \
639 jb++; \ 639 jb++; \
640 jb_lt_max= jb < jb_max; \ 640 jb_lt_max= jb < jb_max; \
641 } \ 641 } \
642 } \ 642 } \
643 } \ 643 } \
644 r.maybe_compress (true); \ 644 r.maybe_compress (true); \
645 } \ 645 } \
646 \ 646 \
647 return r; \ 647 return r; \
648 } 648 }
649 649
650 // Note that SM ./ SM needs to take into account the NaN and Inf values 650 // Note that SM ./ SM needs to take into account the NaN and Inf values
651 // implied by the division by zero. 651 // implied by the division by zero.
652 // FIXME: Are the NaNs double(NaN) or Complex(NaN,Nan) in the complex case? 652 // FIXME: Are the NaNs double(NaN) or Complex(NaN,Nan) in the complex case?
653 #define SPARSE_SMSM_BIN_OPS(R1, R2, M1, M2) \ 653 #define SPARSE_SMSM_BIN_OPS(R1, R2, M1, M2) \
654 SPARSE_SMSM_BIN_OP_1 (R1, operator +, +, M1, M2) \ 654 SPARSE_SMSM_BIN_OP_1 (R1, operator +, +, M1, M2) \
655 SPARSE_SMSM_BIN_OP_1 (R1, operator -, -, M1, M2) \ 655 SPARSE_SMSM_BIN_OP_1 (R1, operator -, -, M1, M2) \
656 SPARSE_SMSM_BIN_OP_2 (R2, product, *, M1, M2) \ 656 SPARSE_SMSM_BIN_OP_2 (R2, product, *, M1, M2) \
657 SPARSE_SMSM_BIN_OP_3 (R2, quotient, /, M1, M2) 657 SPARSE_SMSM_BIN_OP_3 (R2, quotient, /, M1, M2)
658 658
659 // FIXME: this macro duplicates the bodies of the template functions 659 // FIXME: this macro duplicates the bodies of the template functions
660 // defined in the SPARSE_SSM_CMP_OP and SPARSE_SMS_CMP_OP macros. 660 // defined in the SPARSE_SSM_CMP_OP and SPARSE_SMS_CMP_OP macros.
661 661
662 #define SPARSE_SMSM_CMP_OP(F, OP, M1, Z1, C1, M2, Z2, C2) \ 662 #define SPARSE_SMSM_CMP_OP(F, OP, M1, Z1, C1, M2, Z2, C2) \
663 SparseBoolMatrix \ 663 SparseBoolMatrix \
664 F (const M1& m1, const M2& m2) \ 664 F (const M1& m1, const M2& m2) \
665 { \ 665 { \
666 SparseBoolMatrix r; \ 666 SparseBoolMatrix r; \
667 \ 667 \
668 octave_idx_type m1_nr = m1.rows (); \ 668 octave_idx_type m1_nr = m1.rows (); \
669 octave_idx_type m1_nc = m1.cols (); \ 669 octave_idx_type m1_nc = m1.cols (); \
670 \ 670 \
671 octave_idx_type m2_nr = m2.rows (); \ 671 octave_idx_type m2_nr = m2.rows (); \
672 octave_idx_type m2_nc = m2.cols (); \ 672 octave_idx_type m2_nc = m2.cols (); \
673 \ 673 \
674 if (m1_nr == 1 && m1_nc == 1) \ 674 if (m1_nr == 1 && m1_nc == 1) \
675 { \ 675 { \
676 if (C1 (m1.elem (0,0)) OP C2 (Z2)) \ 676 if (C1 (m1.elem (0,0)) OP C2 (Z2)) \
677 { \ 677 { \
678 r = SparseBoolMatrix (m2_nr, m2_nc, true); \ 678 r = SparseBoolMatrix (m2_nr, m2_nc, true); \
679 for (octave_idx_type j = 0; j < m2_nc; j++) \ 679 for (octave_idx_type j = 0; j < m2_nc; j++) \
680 for (octave_idx_type i = m2.cidx (j); i < m2.cidx (j+1); i++) \ 680 for (octave_idx_type i = m2.cidx (j); i < m2.cidx (j+1); i++) \
681 if (! (C1 (m1.elem (0,0)) OP C2 (m2.data (i)))) \ 681 if (! (C1 (m1.elem (0,0)) OP C2 (m2.data (i)))) \
682 r.data (m2.ridx (i) + j * m2_nr) = false; \ 682 r.data (m2.ridx (i) + j * m2_nr) = false; \
683 r.maybe_compress (true); \ 683 r.maybe_compress (true); \
684 } \ 684 } \
685 else \ 685 else \
686 { \ 686 { \
687 r = SparseBoolMatrix (m2_nr, m2_nc, m2.nnz ()); \ 687 r = SparseBoolMatrix (m2_nr, m2_nc, m2.nnz ()); \
688 r.cidx (0) = static_cast<octave_idx_type> (0); \ 688 r.cidx (0) = static_cast<octave_idx_type> (0); \
689 octave_idx_type nel = 0; \ 689 octave_idx_type nel = 0; \
690 for (octave_idx_type j = 0; j < m2_nc; j++) \ 690 for (octave_idx_type j = 0; j < m2_nc; j++) \
691 { \ 691 { \
692 for (octave_idx_type i = m2.cidx (j); i < m2.cidx (j+1); i++) \ 692 for (octave_idx_type i = m2.cidx (j); i < m2.cidx (j+1); i++) \
693 if (C1 (m1.elem (0,0)) OP C2 (m2.data (i))) \ 693 if (C1 (m1.elem (0,0)) OP C2 (m2.data (i))) \
694 { \ 694 { \
695 r.ridx (nel) = m2.ridx (i); \ 695 r.ridx (nel) = m2.ridx (i); \
696 r.data (nel++) = true; \ 696 r.data (nel++) = true; \
697 } \ 697 } \
698 r.cidx (j + 1) = nel; \ 698 r.cidx (j + 1) = nel; \
699 } \ 699 } \
700 r.maybe_compress (false); \ 700 r.maybe_compress (false); \
701 } \ 701 } \
702 } \ 702 } \
703 else if (m2_nr == 1 && m2_nc == 1) \ 703 else if (m2_nr == 1 && m2_nc == 1) \
704 { \ 704 { \
705 if (C1 (Z1) OP C2 (m2.elem (0,0))) \ 705 if (C1 (Z1) OP C2 (m2.elem (0,0))) \
706 { \ 706 { \
707 r = SparseBoolMatrix (m1_nr, m1_nc, true); \ 707 r = SparseBoolMatrix (m1_nr, m1_nc, true); \
708 for (octave_idx_type j = 0; j < m1_nc; j++) \ 708 for (octave_idx_type j = 0; j < m1_nc; j++) \
709 for (octave_idx_type i = m1.cidx (j); i < m1.cidx (j+1); i++) \ 709 for (octave_idx_type i = m1.cidx (j); i < m1.cidx (j+1); i++) \
710 if (! (C1 (m1.data (i)) OP C2 (m2.elem (0,0)))) \ 710 if (! (C1 (m1.data (i)) OP C2 (m2.elem (0,0)))) \
711 r.data (m1.ridx (i) + j * m1_nr) = false; \ 711 r.data (m1.ridx (i) + j * m1_nr) = false; \
712 r.maybe_compress (true); \ 712 r.maybe_compress (true); \
713 } \ 713 } \
714 else \ 714 else \
715 { \ 715 { \
716 r = SparseBoolMatrix (m1_nr, m1_nc, m1.nnz ()); \ 716 r = SparseBoolMatrix (m1_nr, m1_nc, m1.nnz ()); \
717 r.cidx (0) = static_cast<octave_idx_type> (0); \ 717 r.cidx (0) = static_cast<octave_idx_type> (0); \
718 octave_idx_type nel = 0; \ 718 octave_idx_type nel = 0; \
719 for (octave_idx_type j = 0; j < m1_nc; j++) \ 719 for (octave_idx_type j = 0; j < m1_nc; j++) \
720 { \ 720 { \
721 for (octave_idx_type i = m1.cidx (j); i < m1.cidx (j+1); i++) \ 721 for (octave_idx_type i = m1.cidx (j); i < m1.cidx (j+1); i++) \
722 if (C1 (m1.data (i)) OP C2 (m2.elem (0,0))) \ 722 if (C1 (m1.data (i)) OP C2 (m2.elem (0,0))) \
723 { \ 723 { \
724 r.ridx (nel) = m1.ridx (i); \ 724 r.ridx (nel) = m1.ridx (i); \
725 r.data (nel++) = true; \ 725 r.data (nel++) = true; \
726 } \ 726 } \
727 r.cidx (j + 1) = nel; \ 727 r.cidx (j + 1) = nel; \
728 } \ 728 } \
729 r.maybe_compress (false); \ 729 r.maybe_compress (false); \
730 } \ 730 } \
731 } \ 731 } \
732 else if (m1_nr == m2_nr && m1_nc == m2_nc) \ 732 else if (m1_nr == m2_nr && m1_nc == m2_nc) \
733 { \ 733 { \
734 if (m1_nr != 0 || m1_nc != 0) \ 734 if (m1_nr != 0 || m1_nc != 0) \
735 { \ 735 { \
736 if (C1 (Z1) OP C2 (Z2)) \ 736 if (C1 (Z1) OP C2 (Z2)) \
737 { \ 737 { \
738 r = SparseBoolMatrix (m1_nr, m1_nc, true); \ 738 r = SparseBoolMatrix (m1_nr, m1_nc, true); \
739 for (octave_idx_type j = 0; j < m1_nc; j++) \ 739 for (octave_idx_type j = 0; j < m1_nc; j++) \
740 { \ 740 { \
741 octave_idx_type i1 = m1.cidx (j); \ 741 octave_idx_type i1 = m1.cidx (j); \
742 octave_idx_type e1 = m1.cidx (j+1); \ 742 octave_idx_type e1 = m1.cidx (j+1); \
743 octave_idx_type i2 = m2.cidx (j); \ 743 octave_idx_type i2 = m2.cidx (j); \
744 octave_idx_type e2 = m2.cidx (j+1); \ 744 octave_idx_type e2 = m2.cidx (j+1); \
745 while (i1 < e1 || i2 < e2) \ 745 while (i1 < e1 || i2 < e2) \
746 { \ 746 { \
747 if (i1 == e1 || (i2 < e2 && m1.ridx (i1) > m2.ridx (i2))) \ 747 if (i1 == e1 || (i2 < e2 && m1.ridx (i1) > m2.ridx (i2))) \
748 { \ 748 { \
749 if (! (C1 (Z1) OP C2 (m2.data (i2)))) \ 749 if (! (C1 (Z1) OP C2 (m2.data (i2)))) \
750 r.data (m2.ridx (i2) + j * m1_nr) = false; \ 750 r.data (m2.ridx (i2) + j * m1_nr) = false; \
751 i2++; \ 751 i2++; \
752 } \ 752 } \
753 else if (i2 == e2 || m1.ridx (i1) < m2.ridx (i2)) \ 753 else if (i2 == e2 || m1.ridx (i1) < m2.ridx (i2)) \
754 { \ 754 { \
755 if (! (C1 (m1.data (i1)) OP C2 (Z2))) \ 755 if (! (C1 (m1.data (i1)) OP C2 (Z2))) \
756 r.data (m1.ridx (i1) + j * m1_nr) = false; \ 756 r.data (m1.ridx (i1) + j * m1_nr) = false; \
757 i1++; \ 757 i1++; \
758 } \ 758 } \
759 else \ 759 else \
760 { \ 760 { \
761 if (! (C1 (m1.data (i1)) OP C2 (m2.data (i2)))) \ 761 if (! (C1 (m1.data (i1)) OP C2 (m2.data (i2)))) \
762 r.data (m1.ridx (i1) + j * m1_nr) = false; \ 762 r.data (m1.ridx (i1) + j * m1_nr) = false; \
763 i1++; \ 763 i1++; \
764 i2++; \ 764 i2++; \
765 } \ 765 } \
766 } \ 766 } \
767 } \ 767 } \
768 r.maybe_compress (true); \ 768 r.maybe_compress (true); \
769 } \ 769 } \
770 else \ 770 else \
771 { \ 771 { \
772 r = SparseBoolMatrix (m1_nr, m1_nc, m1.nnz () + m2.nnz ()); \ 772 r = SparseBoolMatrix (m1_nr, m1_nc, m1.nnz () + m2.nnz ()); \
773 r.cidx (0) = static_cast<octave_idx_type> (0); \ 773 r.cidx (0) = static_cast<octave_idx_type> (0); \
774 octave_idx_type nel = 0; \ 774 octave_idx_type nel = 0; \
775 for (octave_idx_type j = 0; j < m1_nc; j++) \ 775 for (octave_idx_type j = 0; j < m1_nc; j++) \
776 { \ 776 { \
777 octave_idx_type i1 = m1.cidx (j); \ 777 octave_idx_type i1 = m1.cidx (j); \
778 octave_idx_type e1 = m1.cidx (j+1); \ 778 octave_idx_type e1 = m1.cidx (j+1); \
779 octave_idx_type i2 = m2.cidx (j); \ 779 octave_idx_type i2 = m2.cidx (j); \
780 octave_idx_type e2 = m2.cidx (j+1); \ 780 octave_idx_type e2 = m2.cidx (j+1); \
781 while (i1 < e1 || i2 < e2) \ 781 while (i1 < e1 || i2 < e2) \
782 { \ 782 { \
783 if (i1 == e1 || (i2 < e2 && m1.ridx (i1) > m2.ridx (i2))) \ 783 if (i1 == e1 || (i2 < e2 && m1.ridx (i1) > m2.ridx (i2))) \
784 { \ 784 { \
785 if (C1 (Z1) OP C2 (m2.data (i2))) \ 785 if (C1 (Z1) OP C2 (m2.data (i2))) \
786 { \ 786 { \
787 r.ridx (nel) = m2.ridx (i2); \ 787 r.ridx (nel) = m2.ridx (i2); \
788 r.data (nel++) = true; \ 788 r.data (nel++) = true; \
789 } \ 789 } \
790 i2++; \ 790 i2++; \
791 } \ 791 } \
792 else if (i2 == e2 || m1.ridx (i1) < m2.ridx (i2)) \ 792 else if (i2 == e2 || m1.ridx (i1) < m2.ridx (i2)) \
793 { \ 793 { \
794 if (C1 (m1.data (i1)) OP C2 (Z2)) \ 794 if (C1 (m1.data (i1)) OP C2 (Z2)) \
795 { \ 795 { \
796 r.ridx (nel) = m1.ridx (i1); \ 796 r.ridx (nel) = m1.ridx (i1); \
797 r.data (nel++) = true; \ 797 r.data (nel++) = true; \
798 } \ 798 } \
799 i1++; \ 799 i1++; \
800 } \ 800 } \
801 else \ 801 else \
802 { \ 802 { \
803 if (C1 (m1.data (i1)) OP C2 (m2.data (i2))) \ 803 if (C1 (m1.data (i1)) OP C2 (m2.data (i2))) \
804 { \ 804 { \
805 r.ridx (nel) = m1.ridx (i1); \ 805 r.ridx (nel) = m1.ridx (i1); \
806 r.data (nel++) = true; \ 806 r.data (nel++) = true; \
807 } \ 807 } \
808 i1++; \ 808 i1++; \
809 i2++; \ 809 i2++; \
810 } \ 810 } \
811 } \ 811 } \
812 r.cidx (j + 1) = nel; \ 812 r.cidx (j + 1) = nel; \
813 } \ 813 } \
814 r.maybe_compress (false); \ 814 r.maybe_compress (false); \
815 } \ 815 } \
816 } \ 816 } \
817 } \ 817 } \
818 else \ 818 else \
819 { \ 819 { \
820 if ((m1_nr != 0 || m1_nc != 0) && (m2_nr != 0 || m2_nc != 0)) \ 820 if ((m1_nr != 0 || m1_nc != 0) && (m2_nr != 0 || m2_nc != 0)) \
821 err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \ 821 err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \
822 } \ 822 } \
823 return r; \ 823 return r; \
824 } 824 }
825 825
826 #define SPARSE_SMSM_CMP_OPS(M1, Z1, C1, M2, Z2, C2) \ 826 #define SPARSE_SMSM_CMP_OPS(M1, Z1, C1, M2, Z2, C2) \
827 SPARSE_SMSM_CMP_OP (mx_el_lt, <, M1, Z1, , M2, Z2, ) \ 827 SPARSE_SMSM_CMP_OP (mx_el_lt, <, M1, Z1, , M2, Z2, ) \
828 SPARSE_SMSM_CMP_OP (mx_el_le, <=, M1, Z1, , M2, Z2, ) \ 828 SPARSE_SMSM_CMP_OP (mx_el_le, <=, M1, Z1, , M2, Z2, ) \
829 SPARSE_SMSM_CMP_OP (mx_el_ge, >=, M1, Z1, , M2, Z2, ) \ 829 SPARSE_SMSM_CMP_OP (mx_el_ge, >=, M1, Z1, , M2, Z2, ) \
830 SPARSE_SMSM_CMP_OP (mx_el_gt, >, M1, Z1, , M2, Z2, ) \ 830 SPARSE_SMSM_CMP_OP (mx_el_gt, >, M1, Z1, , M2, Z2, ) \
831 SPARSE_SMSM_CMP_OP (mx_el_eq, ==, M1, Z1, , M2, Z2, ) \ 831 SPARSE_SMSM_CMP_OP (mx_el_eq, ==, M1, Z1, , M2, Z2, ) \
832 SPARSE_SMSM_CMP_OP (mx_el_ne, !=, M1, Z1, , M2, Z2, ) 832 SPARSE_SMSM_CMP_OP (mx_el_ne, !=, M1, Z1, , M2, Z2, )
833 833
834 #define SPARSE_SMSM_EQNE_OPS(M1, Z1, C1, M2, Z2, C2) \ 834 #define SPARSE_SMSM_EQNE_OPS(M1, Z1, C1, M2, Z2, C2) \
835 SPARSE_SMSM_CMP_OP (mx_el_eq, ==, M1, Z1, , M2, Z2, ) \ 835 SPARSE_SMSM_CMP_OP (mx_el_eq, ==, M1, Z1, , M2, Z2, ) \
836 SPARSE_SMSM_CMP_OP (mx_el_ne, !=, M1, Z1, , M2, Z2, ) 836 SPARSE_SMSM_CMP_OP (mx_el_ne, !=, M1, Z1, , M2, Z2, )
837 837
838 // FIXME: this macro duplicates the bodies of the template functions 838 // FIXME: this macro duplicates the bodies of the template functions
839 // defined in the SPARSE_SSM_BOOL_OP and SPARSE_SMS_BOOL_OP macros. 839 // defined in the SPARSE_SSM_BOOL_OP and SPARSE_SMS_BOOL_OP macros.
840 840
841 #define SPARSE_SMSM_BOOL_OP(F, OP, M1, M2, LHS_ZERO, RHS_ZERO) \ 841 #define SPARSE_SMSM_BOOL_OP(F, OP, M1, M2, LHS_ZERO, RHS_ZERO) \
842 SparseBoolMatrix \ 842 SparseBoolMatrix \
843 F (const M1& m1, const M2& m2) \ 843 F (const M1& m1, const M2& m2) \
844 { \ 844 { \
845 SparseBoolMatrix r; \ 845 SparseBoolMatrix r; \
846 \ 846 \
847 octave_idx_type m1_nr = m1.rows (); \ 847 octave_idx_type m1_nr = m1.rows (); \
848 octave_idx_type m1_nc = m1.cols (); \ 848 octave_idx_type m1_nc = m1.cols (); \
849 \ 849 \
850 octave_idx_type m2_nr = m2.rows (); \ 850 octave_idx_type m2_nr = m2.rows (); \
851 octave_idx_type m2_nc = m2.cols (); \ 851 octave_idx_type m2_nc = m2.cols (); \
852 \ 852 \
853 if (m1_nr == 1 && m1_nc == 1) \ 853 if (m1_nr == 1 && m1_nc == 1) \
854 { \ 854 { \
855 if (m2_nr > 0 && m2_nc > 0) \ 855 if (m2_nr > 0 && m2_nc > 0) \
856 { \ 856 { \
857 if ((m1.elem (0,0) != LHS_ZERO) OP RHS_ZERO) \ 857 if ((m1.elem (0,0) != LHS_ZERO) OP RHS_ZERO) \
858 { \ 858 { \
859 r = SparseBoolMatrix (m2_nr, m2_nc, true); \ 859 r = SparseBoolMatrix (m2_nr, m2_nc, true); \
860 for (octave_idx_type j = 0; j < m2_nc; j++) \ 860 for (octave_idx_type j = 0; j < m2_nc; j++) \
861 for (octave_idx_type i = m2.cidx (j); i < m2.cidx (j+1); i++) \ 861 for (octave_idx_type i = m2.cidx (j); i < m2.cidx (j+1); i++) \
862 if (! ((m1.elem (0,0) != LHS_ZERO) OP (m2.data (i) != RHS_ZERO))) \ 862 if (! ((m1.elem (0,0) != LHS_ZERO) OP (m2.data (i) != RHS_ZERO))) \
863 r.data (m2.ridx (i) + j * m2_nr) = false; \ 863 r.data (m2.ridx (i) + j * m2_nr) = false; \
864 r.maybe_compress (true); \ 864 r.maybe_compress (true); \
865 } \ 865 } \
866 else \ 866 else \
867 { \ 867 { \
868 r = SparseBoolMatrix (m2_nr, m2_nc, m2.nnz ()); \ 868 r = SparseBoolMatrix (m2_nr, m2_nc, m2.nnz ()); \
869 r.cidx (0) = static_cast<octave_idx_type> (0); \ 869 r.cidx (0) = static_cast<octave_idx_type> (0); \
870 octave_idx_type nel = 0; \ 870 octave_idx_type nel = 0; \
871 for (octave_idx_type j = 0; j < m2_nc; j++) \ 871 for (octave_idx_type j = 0; j < m2_nc; j++) \
872 { \ 872 { \
873 for (octave_idx_type i = m2.cidx (j); i < m2.cidx (j+1); i++) \ 873 for (octave_idx_type i = m2.cidx (j); i < m2.cidx (j+1); i++) \
874 if ((m1.elem (0,0) != LHS_ZERO) OP (m2.data (i) != RHS_ZERO)) \ 874 if ((m1.elem (0,0) != LHS_ZERO) OP (m2.data (i) != RHS_ZERO)) \
875 { \ 875 { \
876 r.ridx (nel) = m2.ridx (i); \ 876 r.ridx (nel) = m2.ridx (i); \
877 r.data (nel++) = true; \ 877 r.data (nel++) = true; \
878 } \ 878 } \
879 r.cidx (j + 1) = nel; \ 879 r.cidx (j + 1) = nel; \
880 } \ 880 } \
881 r.maybe_compress (false); \ 881 r.maybe_compress (false); \
882 } \ 882 } \
883 } \ 883 } \
884 } \ 884 } \
885 else if (m2_nr == 1 && m2_nc == 1) \ 885 else if (m2_nr == 1 && m2_nc == 1) \
886 { \ 886 { \
887 if (m1_nr > 0 && m1_nc > 0) \ 887 if (m1_nr > 0 && m1_nc > 0) \
888 { \ 888 { \
889 if (LHS_ZERO OP (m2.elem (0,0) != RHS_ZERO)) \ 889 if (LHS_ZERO OP (m2.elem (0,0) != RHS_ZERO)) \
890 { \ 890 { \
891 r = SparseBoolMatrix (m1_nr, m1_nc, true); \ 891 r = SparseBoolMatrix (m1_nr, m1_nc, true); \
892 for (octave_idx_type j = 0; j < m1_nc; j++) \ 892 for (octave_idx_type j = 0; j < m1_nc; j++) \
893 for (octave_idx_type i = m1.cidx (j); i < m1.cidx (j+1); i++) \ 893 for (octave_idx_type i = m1.cidx (j); i < m1.cidx (j+1); i++) \
894 if (! ((m1.data (i) != LHS_ZERO) OP (m2.elem (0,0) != RHS_ZERO))) \ 894 if (! ((m1.data (i) != LHS_ZERO) OP (m2.elem (0,0) != RHS_ZERO))) \
895 r.data (m1.ridx (i) + j * m1_nr) = false; \ 895 r.data (m1.ridx (i) + j * m1_nr) = false; \
896 r.maybe_compress (true); \ 896 r.maybe_compress (true); \
897 } \ 897 } \
898 else \ 898 else \
899 { \ 899 { \
900 r = SparseBoolMatrix (m1_nr, m1_nc, m1.nnz ()); \ 900 r = SparseBoolMatrix (m1_nr, m1_nc, m1.nnz ()); \
901 r.cidx (0) = static_cast<octave_idx_type> (0); \ 901 r.cidx (0) = static_cast<octave_idx_type> (0); \
902 octave_idx_type nel = 0; \ 902 octave_idx_type nel = 0; \
903 for (octave_idx_type j = 0; j < m1_nc; j++) \ 903 for (octave_idx_type j = 0; j < m1_nc; j++) \
904 { \ 904 { \
905 for (octave_idx_type i = m1.cidx (j); i < m1.cidx (j+1); i++) \ 905 for (octave_idx_type i = m1.cidx (j); i < m1.cidx (j+1); i++) \
906 if ((m1.data (i) != LHS_ZERO) OP (m2.elem (0,0) != RHS_ZERO)) \ 906 if ((m1.data (i) != LHS_ZERO) OP (m2.elem (0,0) != RHS_ZERO)) \
907 { \ 907 { \
908 r.ridx (nel) = m1.ridx (i); \ 908 r.ridx (nel) = m1.ridx (i); \
909 r.data (nel++) = true; \ 909 r.data (nel++) = true; \
910 } \ 910 } \
911 r.cidx (j + 1) = nel; \ 911 r.cidx (j + 1) = nel; \
912 } \ 912 } \
913 r.maybe_compress (false); \ 913 r.maybe_compress (false); \
914 } \ 914 } \
915 } \ 915 } \
916 } \ 916 } \
917 else if (m1_nr == m2_nr && m1_nc == m2_nc) \ 917 else if (m1_nr == m2_nr && m1_nc == m2_nc) \
918 { \ 918 { \
919 if (m1_nr != 0 || m1_nc != 0) \ 919 if (m1_nr != 0 || m1_nc != 0) \
920 { \ 920 { \
921 r = SparseBoolMatrix (m1_nr, m1_nc, m1.nnz () + m2.nnz ()); \ 921 r = SparseBoolMatrix (m1_nr, m1_nc, m1.nnz () + m2.nnz ()); \
922 r.cidx (0) = static_cast<octave_idx_type> (0); \ 922 r.cidx (0) = static_cast<octave_idx_type> (0); \
923 octave_idx_type nel = 0; \ 923 octave_idx_type nel = 0; \
924 for (octave_idx_type j = 0; j < m1_nc; j++) \ 924 for (octave_idx_type j = 0; j < m1_nc; j++) \
925 { \ 925 { \
926 octave_idx_type i1 = m1.cidx (j); \ 926 octave_idx_type i1 = m1.cidx (j); \
927 octave_idx_type e1 = m1.cidx (j+1); \ 927 octave_idx_type e1 = m1.cidx (j+1); \
928 octave_idx_type i2 = m2.cidx (j); \ 928 octave_idx_type i2 = m2.cidx (j); \
929 octave_idx_type e2 = m2.cidx (j+1); \ 929 octave_idx_type e2 = m2.cidx (j+1); \
930 while (i1 < e1 || i2 < e2) \ 930 while (i1 < e1 || i2 < e2) \
931 { \ 931 { \
932 if (i1 == e1 || (i2 < e2 && m1.ridx (i1) > m2.ridx (i2))) \ 932 if (i1 == e1 || (i2 < e2 && m1.ridx (i1) > m2.ridx (i2))) \
933 { \ 933 { \
934 if (LHS_ZERO OP m2.data (i2) != RHS_ZERO) \ 934 if (LHS_ZERO OP m2.data (i2) != RHS_ZERO) \
935 { \ 935 { \
936 r.ridx (nel) = m2.ridx (i2); \ 936 r.ridx (nel) = m2.ridx (i2); \
937 r.data (nel++) = true; \ 937 r.data (nel++) = true; \
938 } \ 938 } \
939 i2++; \ 939 i2++; \
940 } \ 940 } \
941 else if (i2 == e2 || m1.ridx (i1) < m2.ridx (i2)) \ 941 else if (i2 == e2 || m1.ridx (i1) < m2.ridx (i2)) \
942 { \ 942 { \
943 if (m1.data (i1) != LHS_ZERO OP RHS_ZERO) \ 943 if (m1.data (i1) != LHS_ZERO OP RHS_ZERO) \
944 { \ 944 { \
945 r.ridx (nel) = m1.ridx (i1); \ 945 r.ridx (nel) = m1.ridx (i1); \
946 r.data (nel++) = true; \ 946 r.data (nel++) = true; \
947 } \ 947 } \
948 i1++; \ 948 i1++; \
949 } \ 949 } \
950 else \ 950 else \
951 { \ 951 { \
952 if (m1.data (i1) != LHS_ZERO OP m2.data (i2) != RHS_ZERO) \ 952 if (m1.data (i1) != LHS_ZERO OP m2.data (i2) != RHS_ZERO) \
953 { \ 953 { \
954 r.ridx (nel) = m1.ridx (i1); \ 954 r.ridx (nel) = m1.ridx (i1); \
955 r.data (nel++) = true; \ 955 r.data (nel++) = true; \
956 } \ 956 } \
957 i1++; \ 957 i1++; \
958 i2++; \ 958 i2++; \
959 } \ 959 } \
960 } \ 960 } \
961 r.cidx (j + 1) = nel; \ 961 r.cidx (j + 1) = nel; \
962 } \ 962 } \
963 r.maybe_compress (false); \ 963 r.maybe_compress (false); \
964 } \ 964 } \
965 } \ 965 } \
966 else \ 966 else \
967 { \ 967 { \
968 if ((m1_nr != 0 || m1_nc != 0) && (m2_nr != 0 || m2_nc != 0)) \ 968 if ((m1_nr != 0 || m1_nc != 0) && (m2_nr != 0 || m2_nc != 0)) \
969 err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \ 969 err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \
970 } \ 970 } \
971 return r; \ 971 return r; \
972 } 972 }
973 973
974 #define SPARSE_SMSM_BOOL_OPS2(M1, M2, LHS_ZERO, RHS_ZERO) \ 974 #define SPARSE_SMSM_BOOL_OPS2(M1, M2, LHS_ZERO, RHS_ZERO) \
975 SPARSE_SMSM_BOOL_OP (mx_el_and, &&, M1, M2, LHS_ZERO, RHS_ZERO) \ 975 SPARSE_SMSM_BOOL_OP (mx_el_and, &&, M1, M2, LHS_ZERO, RHS_ZERO) \
976 SPARSE_SMSM_BOOL_OP (mx_el_or, ||, M1, M2, LHS_ZERO, RHS_ZERO) 976 SPARSE_SMSM_BOOL_OP (mx_el_or, ||, M1, M2, LHS_ZERO, RHS_ZERO)
977 977
978 #define SPARSE_SMSM_BOOL_OPS(M1, M2, ZERO) \ 978 #define SPARSE_SMSM_BOOL_OPS(M1, M2, ZERO) \
979 SPARSE_SMSM_BOOL_OPS2(M1, M2, ZERO, ZERO) 979 SPARSE_SMSM_BOOL_OPS2(M1, M2, ZERO, ZERO)
980 980
981 // matrix by sparse matrix operations. 981 // matrix by sparse matrix operations.
982 982
983 #define SPARSE_MSM_BIN_OP_1(R, F, OP, M1, M2) \ 983 #define SPARSE_MSM_BIN_OP_1(R, F, OP, M1, M2) \
984 R \ 984 R \
985 F (const M1& m1, const M2& m2) \ 985 F (const M1& m1, const M2& m2) \
986 { \ 986 { \
987 R r; \ 987 R r; \
988 \ 988 \
989 octave_idx_type m1_nr = m1.rows (); \ 989 octave_idx_type m1_nr = m1.rows (); \
990 octave_idx_type m1_nc = m1.cols (); \ 990 octave_idx_type m1_nc = m1.cols (); \
991 \ 991 \
992 octave_idx_type m2_nr = m2.rows (); \ 992 octave_idx_type m2_nr = m2.rows (); \
993 octave_idx_type m2_nc = m2.cols (); \ 993 octave_idx_type m2_nc = m2.cols (); \
994 \ 994 \
995 if (m2_nr == 1 && m2_nc == 1) \ 995 if (m2_nr == 1 && m2_nc == 1) \
996 r = R (m1 OP m2.elem (0,0)); \ 996 r = R (m1 OP m2.elem (0,0)); \
997 else if (m1_nr != m2_nr || m1_nc != m2_nc) \ 997 else if (m1_nr != m2_nr || m1_nc != m2_nc) \
998 err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \ 998 err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \
999 else \ 999 else \
1000 { \ 1000 { \
1001 r = R (F (m1, m2.matrix_value ())); \ 1001 r = R (F (m1, m2.matrix_value ())); \
1002 } \ 1002 } \
1003 return r; \ 1003 return r; \
1004 } 1004 }
1005 1005
1006 #define SPARSE_MSM_BIN_OP_2(R, F, OP, M1, M2) \ 1006 #define SPARSE_MSM_BIN_OP_2(R, F, OP, M1, M2) \
1007 R \ 1007 R \
1008 F (const M1& m1, const M2& m2) \ 1008 F (const M1& m1, const M2& m2) \
1009 { \ 1009 { \
1010 R r; \ 1010 R r; \
1011 \ 1011 \
1012 octave_idx_type m1_nr = m1.rows (); \ 1012 octave_idx_type m1_nr = m1.rows (); \
1013 octave_idx_type m1_nc = m1.cols (); \ 1013 octave_idx_type m1_nc = m1.cols (); \
1014 \ 1014 \
1015 octave_idx_type m2_nr = m2.rows (); \ 1015 octave_idx_type m2_nr = m2.rows (); \
1016 octave_idx_type m2_nc = m2.cols (); \ 1016 octave_idx_type m2_nc = m2.cols (); \
1017 \ 1017 \
1018 if (m2_nr == 1 && m2_nc == 1) \ 1018 if (m2_nr == 1 && m2_nc == 1) \
1019 r = R (m1 OP m2.elem (0,0)); \ 1019 r = R (m1 OP m2.elem (0,0)); \
1020 else if (m1_nr != m2_nr || m1_nc != m2_nc) \ 1020 else if (m1_nr != m2_nr || m1_nc != m2_nc) \
1021 err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \ 1021 err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \
1022 else \ 1022 else \
1023 { \ 1023 { \
1024 if (do_mx_check (m1, mx_inline_all_finite<M1::element_type>)) \ 1024 if (do_mx_check (m1, mx_inline_all_finite<M1::element_type>)) \
1025 { \ 1025 { \
1026 /* Sparsity pattern is preserved. */ \ 1026 /* Sparsity pattern is preserved. */ \
1027 octave_idx_type m2_nz = m2.nnz (); \ 1027 octave_idx_type m2_nz = m2.nnz (); \
1028 r = R (m2_nr, m2_nc, m2_nz); \ 1028 r = R (m2_nr, m2_nc, m2_nz); \
1029 for (octave_idx_type j = 0, k = 0; j < m2_nc; j++) \ 1029 for (octave_idx_type j = 0, k = 0; j < m2_nc; j++) \
1030 { \ 1030 { \
1031 octave_quit (); \ 1031 octave_quit (); \
1032 for (octave_idx_type i = m2.cidx (j); i < m2.cidx (j+1); i++) \ 1032 for (octave_idx_type i = m2.cidx (j); i < m2.cidx (j+1); i++) \
1033 { \ 1033 { \
1034 octave_idx_type mri = m2.ridx (i); \ 1034 octave_idx_type mri = m2.ridx (i); \
1035 R::element_type x = m1(mri, j) OP m2.data (i); \ 1035 R::element_type x = m1(mri, j) OP m2.data (i); \
1036 if (x != 0.0) \ 1036 if (x != 0.0) \
1037 { \ 1037 { \
1038 r.xdata (k) = x; \ 1038 r.xdata (k) = x; \
1039 r.xridx (k) = m2.ridx (i); \ 1039 r.xridx (k) = m2.ridx (i); \
1040 k++; \ 1040 k++; \
1041 } \ 1041 } \
1042 } \ 1042 } \
1043 r.xcidx (j+1) = k; \ 1043 r.xcidx (j+1) = k; \
1044 } \ 1044 } \
1045 r.maybe_compress (false); \ 1045 r.maybe_compress (false); \
1046 return r; \ 1046 return r; \
1047 } \ 1047 } \
1048 else \ 1048 else \
1049 r = R (F (m1, m2.matrix_value ())); \ 1049 r = R (F (m1, m2.matrix_value ())); \
1050 } \ 1050 } \
1051 \ 1051 \
1052 return r; \ 1052 return r; \
1053 } 1053 }
1054 1054
1055 // FIXME: Pass a specific ZERO value 1055 // FIXME: Pass a specific ZERO value
1056 #define SPARSE_MSM_BIN_OPS(R1, R2, M1, M2) \ 1056 #define SPARSE_MSM_BIN_OPS(R1, R2, M1, M2) \
1057 SPARSE_MSM_BIN_OP_1 (R1, operator +, +, M1, M2) \ 1057 SPARSE_MSM_BIN_OP_1 (R1, operator +, +, M1, M2) \
1058 SPARSE_MSM_BIN_OP_1 (R1, operator -, -, M1, M2) \ 1058 SPARSE_MSM_BIN_OP_1 (R1, operator -, -, M1, M2) \
1059 SPARSE_MSM_BIN_OP_2 (R2, product, *, M1, M2) \ 1059 SPARSE_MSM_BIN_OP_2 (R2, product, *, M1, M2) \
1060 SPARSE_MSM_BIN_OP_1 (R2, quotient, /, M1, M2) 1060 SPARSE_MSM_BIN_OP_1 (R2, quotient, /, M1, M2)
1061 1061
1062 #define SPARSE_MSM_CMP_OP(F, OP, M1, C1, M2, C2) \ 1062 #define SPARSE_MSM_CMP_OP(F, OP, M1, C1, M2, C2) \
1063 SparseBoolMatrix \ 1063 SparseBoolMatrix \
1064 F (const M1& m1, const M2& m2) \ 1064 F (const M1& m1, const M2& m2) \
1065 { \ 1065 { \
1066 SparseBoolMatrix r; \ 1066 SparseBoolMatrix r; \
1067 \ 1067 \
1068 octave_idx_type m1_nr = m1.rows (); \ 1068 octave_idx_type m1_nr = m1.rows (); \
1069 octave_idx_type m1_nc = m1.cols (); \ 1069 octave_idx_type m1_nc = m1.cols (); \
1070 \ 1070 \
1071 octave_idx_type m2_nr = m2.rows (); \ 1071 octave_idx_type m2_nr = m2.rows (); \
1072 octave_idx_type m2_nc = m2.cols (); \ 1072 octave_idx_type m2_nc = m2.cols (); \
1073 \ 1073 \
1074 if (m2_nr == 1 && m2_nc == 1) \ 1074 if (m2_nr == 1 && m2_nc == 1) \
1075 r = SparseBoolMatrix (F (m1, m2.elem (0,0))); \ 1075 r = SparseBoolMatrix (F (m1, m2.elem (0,0))); \
1076 else if (m1_nr == m2_nr && m1_nc == m2_nc) \ 1076 else if (m1_nr == m2_nr && m1_nc == m2_nc) \
1077 { \ 1077 { \
1078 if (m1_nr != 0 || m1_nc != 0) \ 1078 if (m1_nr != 0 || m1_nc != 0) \
1079 { \ 1079 { \
1080 /* Count num of nonzero elements */ \ 1080 /* Count num of nonzero elements */ \
1081 octave_idx_type nel = 0; \ 1081 octave_idx_type nel = 0; \
1082 for (octave_idx_type j = 0; j < m1_nc; j++) \ 1082 for (octave_idx_type j = 0; j < m1_nc; j++) \
1083 for (octave_idx_type i = 0; i < m1_nr; i++) \ 1083 for (octave_idx_type i = 0; i < m1_nr; i++) \
1084 if (C1 (m1.elem (i, j)) OP C2 (m2.elem (i, j))) \ 1084 if (C1 (m1.elem (i, j)) OP C2 (m2.elem (i, j))) \
1085 nel++; \ 1085 nel++; \
1086 \ 1086 \
1087 r = SparseBoolMatrix (m1_nr, m1_nc, nel); \ 1087 r = SparseBoolMatrix (m1_nr, m1_nc, nel); \
1088 \ 1088 \
1089 octave_idx_type ii = 0; \ 1089 octave_idx_type ii = 0; \
1090 r.cidx (0) = 0; \ 1090 r.cidx (0) = 0; \
1091 for (octave_idx_type j = 0; j < m1_nc; j++) \ 1091 for (octave_idx_type j = 0; j < m1_nc; j++) \
1092 { \ 1092 { \
1093 for (octave_idx_type i = 0; i < m1_nr; i++) \ 1093 for (octave_idx_type i = 0; i < m1_nr; i++) \
1094 { \ 1094 { \
1095 bool el = C1 (m1.elem (i, j)) OP C2 (m2.elem (i, j)); \ 1095 bool el = C1 (m1.elem (i, j)) OP C2 (m2.elem (i, j)); \
1096 if (el) \ 1096 if (el) \
1097 { \ 1097 { \
1098 r.data (ii) = el; \ 1098 r.data (ii) = el; \
1099 r.ridx (ii++) = i; \ 1099 r.ridx (ii++) = i; \
1100 } \ 1100 } \
1101 } \ 1101 } \
1102 r.cidx (j+1) = ii; \ 1102 r.cidx (j+1) = ii; \
1103 } \ 1103 } \
1104 } \ 1104 } \
1105 } \ 1105 } \
1106 else \ 1106 else \
1107 { \ 1107 { \
1108 if ((m1_nr != 0 || m1_nc != 0) && (m2_nr != 0 || m2_nc != 0)) \ 1108 if ((m1_nr != 0 || m1_nc != 0) && (m2_nr != 0 || m2_nc != 0)) \
1109 err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \ 1109 err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \
1110 } \ 1110 } \
1111 return r; \ 1111 return r; \
1112 } 1112 }
1113 1113
1114 #define SPARSE_MSM_CMP_OPS(M1, Z1, C1, M2, Z2, C2) \ 1114 #define SPARSE_MSM_CMP_OPS(M1, Z1, C1, M2, Z2, C2) \
1115 SPARSE_MSM_CMP_OP (mx_el_lt, <, M1, , M2, ) \ 1115 SPARSE_MSM_CMP_OP (mx_el_lt, <, M1, , M2, ) \
1116 SPARSE_MSM_CMP_OP (mx_el_le, <=, M1, , M2, ) \ 1116 SPARSE_MSM_CMP_OP (mx_el_le, <=, M1, , M2, ) \
1117 SPARSE_MSM_CMP_OP (mx_el_ge, >=, M1, , M2, ) \ 1117 SPARSE_MSM_CMP_OP (mx_el_ge, >=, M1, , M2, ) \
1118 SPARSE_MSM_CMP_OP (mx_el_gt, >, M1, , M2, ) \ 1118 SPARSE_MSM_CMP_OP (mx_el_gt, >, M1, , M2, ) \
1119 SPARSE_MSM_CMP_OP (mx_el_eq, ==, M1, , M2, ) \ 1119 SPARSE_MSM_CMP_OP (mx_el_eq, ==, M1, , M2, ) \
1120 SPARSE_MSM_CMP_OP (mx_el_ne, !=, M1, , M2, ) 1120 SPARSE_MSM_CMP_OP (mx_el_ne, !=, M1, , M2, )
1121 1121
1122 #define SPARSE_MSM_EQNE_OPS(M1, Z1, C1, M2, Z2, C2) \ 1122 #define SPARSE_MSM_EQNE_OPS(M1, Z1, C1, M2, Z2, C2) \
1123 SPARSE_MSM_CMP_OP (mx_el_eq, ==, M1, , M2, ) \ 1123 SPARSE_MSM_CMP_OP (mx_el_eq, ==, M1, , M2, ) \
1124 SPARSE_MSM_CMP_OP (mx_el_ne, !=, M1, , M2, ) 1124 SPARSE_MSM_CMP_OP (mx_el_ne, !=, M1, , M2, )
1125 1125
1126 #define SPARSE_MSM_BOOL_OP(F, OP, M1, M2, LHS_ZERO, RHS_ZERO) \ 1126 #define SPARSE_MSM_BOOL_OP(F, OP, M1, M2, LHS_ZERO, RHS_ZERO) \
1127 SparseBoolMatrix \ 1127 SparseBoolMatrix \
1128 F (const M1& m1, const M2& m2) \ 1128 F (const M1& m1, const M2& m2) \
1129 { \ 1129 { \
1130 SparseBoolMatrix r; \ 1130 SparseBoolMatrix r; \
1131 \ 1131 \
1132 octave_idx_type m1_nr = m1.rows (); \ 1132 octave_idx_type m1_nr = m1.rows (); \
1133 octave_idx_type m1_nc = m1.cols (); \ 1133 octave_idx_type m1_nc = m1.cols (); \
1134 \ 1134 \
1135 octave_idx_type m2_nr = m2.rows (); \ 1135 octave_idx_type m2_nr = m2.rows (); \
1136 octave_idx_type m2_nc = m2.cols (); \ 1136 octave_idx_type m2_nc = m2.cols (); \
1137 \ 1137 \
1138 if (m2_nr == 1 && m2_nc == 1) \ 1138 if (m2_nr == 1 && m2_nc == 1) \
1139 r = SparseBoolMatrix (F (m1, m2.elem (0,0))); \ 1139 r = SparseBoolMatrix (F (m1, m2.elem (0,0))); \
1140 else if (m1_nr == m2_nr && m1_nc == m2_nc) \ 1140 else if (m1_nr == m2_nr && m1_nc == m2_nc) \
1141 { \ 1141 { \
1142 if (m1_nr != 0 || m1_nc != 0) \ 1142 if (m1_nr != 0 || m1_nc != 0) \
1143 { \ 1143 { \
1144 /* Count num of nonzero elements */ \ 1144 /* Count num of nonzero elements */ \
1145 octave_idx_type nel = 0; \ 1145 octave_idx_type nel = 0; \
1146 for (octave_idx_type j = 0; j < m1_nc; j++) \ 1146 for (octave_idx_type j = 0; j < m1_nc; j++) \
1147 for (octave_idx_type i = 0; i < m1_nr; i++) \ 1147 for (octave_idx_type i = 0; i < m1_nr; i++) \
1148 if ((m1.elem (i, j) != LHS_ZERO) \ 1148 if ((m1.elem (i, j) != LHS_ZERO) \
1149 OP (m2.elem (i, j) != RHS_ZERO)) \ 1149 OP (m2.elem (i, j) != RHS_ZERO)) \
1150 nel++; \ 1150 nel++; \
1151 \ 1151 \
1152 r = SparseBoolMatrix (m1_nr, m1_nc, nel); \ 1152 r = SparseBoolMatrix (m1_nr, m1_nc, nel); \
1153 \ 1153 \
1154 octave_idx_type ii = 0; \ 1154 octave_idx_type ii = 0; \
1155 r.cidx (0) = 0; \ 1155 r.cidx (0) = 0; \
1156 for (octave_idx_type j = 0; j < m1_nc; j++) \ 1156 for (octave_idx_type j = 0; j < m1_nc; j++) \
1157 { \ 1157 { \
1158 for (octave_idx_type i = 0; i < m1_nr; i++) \ 1158 for (octave_idx_type i = 0; i < m1_nr; i++) \
1159 { \ 1159 { \
1160 bool el = (m1.elem (i, j) != LHS_ZERO) \ 1160 bool el = (m1.elem (i, j) != LHS_ZERO) \
1161 OP (m2.elem (i, j) != RHS_ZERO); \ 1161 OP (m2.elem (i, j) != RHS_ZERO); \
1162 if (el) \ 1162 if (el) \
1163 { \ 1163 { \
1164 r.data (ii) = el; \ 1164 r.data (ii) = el; \
1165 r.ridx (ii++) = i; \ 1165 r.ridx (ii++) = i; \
1166 } \ 1166 } \
1167 } \ 1167 } \
1168 r.cidx (j+1) = ii; \ 1168 r.cidx (j+1) = ii; \
1169 } \ 1169 } \
1170 } \ 1170 } \
1171 } \ 1171 } \
1172 else \ 1172 else \
1173 { \ 1173 { \
1174 if ((m1_nr != 0 || m1_nc != 0) && (m2_nr != 0 || m2_nc != 0)) \ 1174 if ((m1_nr != 0 || m1_nc != 0) && (m2_nr != 0 || m2_nc != 0)) \
1175 err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \ 1175 err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \
1176 } \ 1176 } \
1177 return r; \ 1177 return r; \
1178 } 1178 }
1179 1179
1180 #define SPARSE_MSM_BOOL_OPS2(M1, M2, LHS_ZERO, RHS_ZERO) \ 1180 #define SPARSE_MSM_BOOL_OPS2(M1, M2, LHS_ZERO, RHS_ZERO) \
1181 SPARSE_MSM_BOOL_OP (mx_el_and, &&, M1, M2, LHS_ZERO, RHS_ZERO) \ 1181 SPARSE_MSM_BOOL_OP (mx_el_and, &&, M1, M2, LHS_ZERO, RHS_ZERO) \
1182 SPARSE_MSM_BOOL_OP (mx_el_or, ||, M1, M2, LHS_ZERO, RHS_ZERO) 1182 SPARSE_MSM_BOOL_OP (mx_el_or, ||, M1, M2, LHS_ZERO, RHS_ZERO)
1183 1183
1184 #define SPARSE_MSM_BOOL_OPS(M1, M2, ZERO) \ 1184 #define SPARSE_MSM_BOOL_OPS(M1, M2, ZERO) \
1185 SPARSE_MSM_BOOL_OPS2(M1, M2, ZERO, ZERO) 1185 SPARSE_MSM_BOOL_OPS2(M1, M2, ZERO, ZERO)
1186 1186
1187 // sparse matrix by matrix operations. 1187 // sparse matrix by matrix operations.
1188 1188
1189 #define SPARSE_SMM_BIN_OP_1(R, F, OP, M1, M2) \ 1189 #define SPARSE_SMM_BIN_OP_1(R, F, OP, M1, M2) \
1190 R \ 1190 R \
1191 F (const M1& m1, const M2& m2) \ 1191 F (const M1& m1, const M2& m2) \
1192 { \ 1192 { \
1193 R r; \ 1193 R r; \
1194 \ 1194 \
1195 octave_idx_type m1_nr = m1.rows (); \ 1195 octave_idx_type m1_nr = m1.rows (); \
1196 octave_idx_type m1_nc = m1.cols (); \ 1196 octave_idx_type m1_nc = m1.cols (); \
1197 \ 1197 \
1198 octave_idx_type m2_nr = m2.rows (); \ 1198 octave_idx_type m2_nr = m2.rows (); \
1199 octave_idx_type m2_nc = m2.cols (); \ 1199 octave_idx_type m2_nc = m2.cols (); \
1200 \ 1200 \
1201 if (m1_nr == 1 && m1_nc == 1) \ 1201 if (m1_nr == 1 && m1_nc == 1) \
1202 r = R (m1.elem (0,0) OP m2); \ 1202 r = R (m1.elem (0,0) OP m2); \
1203 else if (m1_nr != m2_nr || m1_nc != m2_nc) \ 1203 else if (m1_nr != m2_nr || m1_nc != m2_nc) \
1204 err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \ 1204 err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \
1205 else \ 1205 else \
1206 { \ 1206 { \
1207 r = R (m1.matrix_value () OP m2); \ 1207 r = R (m1.matrix_value () OP m2); \
1208 } \ 1208 } \
1209 return r; \ 1209 return r; \
1210 } 1210 }
1211 1211
1212 // sm .* m preserves sparsity if m contains no Infs nor Nans. 1212 // sm .* m preserves sparsity if m contains no Infs nor Nans.
1213 #define SPARSE_SMM_BIN_OP_2_CHECK_product(ET) \ 1213 #define SPARSE_SMM_BIN_OP_2_CHECK_product(ET) \
1214 do_mx_check (m2, mx_inline_all_finite<ET>) 1214 do_mx_check (m2, mx_inline_all_finite<ET>)
1215 1215
1216 // sm ./ m preserves sparsity if m contains no NaNs or zeros. 1216 // sm ./ m preserves sparsity if m contains no NaNs or zeros.
1217 #define SPARSE_SMM_BIN_OP_2_CHECK_quotient(ET) \ 1217 #define SPARSE_SMM_BIN_OP_2_CHECK_quotient(ET) \
1218 ! do_mx_check (m2, mx_inline_any_nan<ET>) && m2.nnz () == m2.numel () 1218 ! do_mx_check (m2, mx_inline_any_nan<ET>) && m2.nnz () == m2.numel ()
1219 1219
1220 #define SPARSE_SMM_BIN_OP_2(R, F, OP, M1, M2) \ 1220 #define SPARSE_SMM_BIN_OP_2(R, F, OP, M1, M2) \
1221 R \ 1221 R \
1222 F (const M1& m1, const M2& m2) \ 1222 F (const M1& m1, const M2& m2) \
1223 { \ 1223 { \
1224 R r; \ 1224 R r; \
1225 \ 1225 \
1226 octave_idx_type m1_nr = m1.rows (); \ 1226 octave_idx_type m1_nr = m1.rows (); \
1227 octave_idx_type m1_nc = m1.cols (); \ 1227 octave_idx_type m1_nc = m1.cols (); \
1228 \ 1228 \
1229 octave_idx_type m2_nr = m2.rows (); \ 1229 octave_idx_type m2_nr = m2.rows (); \
1230 octave_idx_type m2_nc = m2.cols (); \ 1230 octave_idx_type m2_nc = m2.cols (); \
1231 \ 1231 \
1232 if (m1_nr == 1 && m1_nc == 1) \ 1232 if (m1_nr == 1 && m1_nc == 1) \
1233 r = R (m1.elem (0,0) OP m2); \ 1233 r = R (m1.elem (0,0) OP m2); \
1234 else if (m1_nr != m2_nr || m1_nc != m2_nc) \ 1234 else if (m1_nr != m2_nr || m1_nc != m2_nc) \
1235 err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \ 1235 err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \
1236 else \ 1236 else \
1237 { \ 1237 { \
1238 if (SPARSE_SMM_BIN_OP_2_CHECK_ ## F(M2::element_type)) \ 1238 if (SPARSE_SMM_BIN_OP_2_CHECK_ ## F(M2::element_type)) \
1239 { \ 1239 { \
1240 /* Sparsity pattern is preserved. */ \ 1240 /* Sparsity pattern is preserved. */ \
1241 octave_idx_type m1_nz = m1.nnz (); \ 1241 octave_idx_type m1_nz = m1.nnz (); \
1242 r = R (m1_nr, m1_nc, m1_nz); \ 1242 r = R (m1_nr, m1_nc, m1_nz); \
1243 for (octave_idx_type j = 0, k = 0; j < m1_nc; j++) \ 1243 for (octave_idx_type j = 0, k = 0; j < m1_nc; j++) \
1244 { \ 1244 { \
1245 octave_quit (); \ 1245 octave_quit (); \
1246 for (octave_idx_type i = m1.cidx (j); i < m1.cidx (j+1); i++) \ 1246 for (octave_idx_type i = m1.cidx (j); i < m1.cidx (j+1); i++) \
1247 { \ 1247 { \
1248 octave_idx_type mri = m1.ridx (i); \ 1248 octave_idx_type mri = m1.ridx (i); \
1249 R::element_type x = m1.data (i) OP m2 (mri, j); \ 1249 R::element_type x = m1.data (i) OP m2 (mri, j); \
1250 if (x != 0.0) \ 1250 if (x != 0.0) \
1251 { \ 1251 { \
1252 r.xdata (k) = x; \ 1252 r.xdata (k) = x; \
1253 r.xridx (k) = m1.ridx (i); \ 1253 r.xridx (k) = m1.ridx (i); \
1254 k++; \ 1254 k++; \
1255 } \ 1255 } \
1256 } \ 1256 } \
1257 r.xcidx (j+1) = k; \ 1257 r.xcidx (j+1) = k; \
1258 } \ 1258 } \
1259 r.maybe_compress (false); \ 1259 r.maybe_compress (false); \
1260 return r; \ 1260 return r; \
1261 } \ 1261 } \
1262 else \ 1262 else \
1263 r = R (F (m1.matrix_value (), m2)); \ 1263 r = R (F (m1.matrix_value (), m2)); \
1264 } \ 1264 } \
1265 \ 1265 \
1266 return r; \ 1266 return r; \
1267 } 1267 }
1268 1268
1269 #define SPARSE_SMM_BIN_OPS(R1, R2, M1, M2) \ 1269 #define SPARSE_SMM_BIN_OPS(R1, R2, M1, M2) \
1270 SPARSE_SMM_BIN_OP_1 (R1, operator +, +, M1, M2) \ 1270 SPARSE_SMM_BIN_OP_1 (R1, operator +, +, M1, M2) \
1271 SPARSE_SMM_BIN_OP_1 (R1, operator -, -, M1, M2) \ 1271 SPARSE_SMM_BIN_OP_1 (R1, operator -, -, M1, M2) \
1272 SPARSE_SMM_BIN_OP_2 (R2, product, *, M1, M2) \ 1272 SPARSE_SMM_BIN_OP_2 (R2, product, *, M1, M2) \
1273 SPARSE_SMM_BIN_OP_2 (R2, quotient, /, M1, M2) 1273 SPARSE_SMM_BIN_OP_2 (R2, quotient, /, M1, M2)
1274 1274
1275 #define SPARSE_SMM_CMP_OP(F, OP, M1, C1, M2, C2) \ 1275 #define SPARSE_SMM_CMP_OP(F, OP, M1, C1, M2, C2) \
1276 SparseBoolMatrix \ 1276 SparseBoolMatrix \
1277 F (const M1& m1, const M2& m2) \ 1277 F (const M1& m1, const M2& m2) \
1278 { \ 1278 { \
1279 SparseBoolMatrix r; \ 1279 SparseBoolMatrix r; \
1280 \ 1280 \
1281 octave_idx_type m1_nr = m1.rows (); \ 1281 octave_idx_type m1_nr = m1.rows (); \
1282 octave_idx_type m1_nc = m1.cols (); \ 1282 octave_idx_type m1_nc = m1.cols (); \
1283 \ 1283 \
1284 octave_idx_type m2_nr = m2.rows (); \ 1284 octave_idx_type m2_nr = m2.rows (); \
1285 octave_idx_type m2_nc = m2.cols (); \ 1285 octave_idx_type m2_nc = m2.cols (); \
1286 \ 1286 \
1287 if (m1_nr == 1 && m1_nc == 1) \ 1287 if (m1_nr == 1 && m1_nc == 1) \
1288 r = SparseBoolMatrix (F (m1.elem (0,0), m2)); \ 1288 r = SparseBoolMatrix (F (m1.elem (0,0), m2)); \
1289 else if (m1_nr == m2_nr && m1_nc == m2_nc) \ 1289 else if (m1_nr == m2_nr && m1_nc == m2_nc) \
1290 { \ 1290 { \
1291 if (m1_nr != 0 || m1_nc != 0) \ 1291 if (m1_nr != 0 || m1_nc != 0) \
1292 { \ 1292 { \
1293 /* Count num of nonzero elements */ \ 1293 /* Count num of nonzero elements */ \
1294 octave_idx_type nel = 0; \ 1294 octave_idx_type nel = 0; \
1295 for (octave_idx_type j = 0; j < m1_nc; j++) \ 1295 for (octave_idx_type j = 0; j < m1_nc; j++) \
1296 for (octave_idx_type i = 0; i < m1_nr; i++) \ 1296 for (octave_idx_type i = 0; i < m1_nr; i++) \
1297 if (C1 (m1.elem (i, j)) OP C2 (m2.elem (i, j))) \ 1297 if (C1 (m1.elem (i, j)) OP C2 (m2.elem (i, j))) \
1298 nel++; \ 1298 nel++; \
1299 \ 1299 \
1300 r = SparseBoolMatrix (m1_nr, m1_nc, nel); \ 1300 r = SparseBoolMatrix (m1_nr, m1_nc, nel); \
1301 \ 1301 \
1302 octave_idx_type ii = 0; \ 1302 octave_idx_type ii = 0; \
1303 r.cidx (0) = 0; \ 1303 r.cidx (0) = 0; \
1304 for (octave_idx_type j = 0; j < m1_nc; j++) \ 1304 for (octave_idx_type j = 0; j < m1_nc; j++) \
1305 { \ 1305 { \
1306 for (octave_idx_type i = 0; i < m1_nr; i++) \ 1306 for (octave_idx_type i = 0; i < m1_nr; i++) \
1307 { \ 1307 { \
1308 bool el = C1 (m1.elem (i, j)) OP C2 (m2.elem (i, j)); \ 1308 bool el = C1 (m1.elem (i, j)) OP C2 (m2.elem (i, j)); \
1309 if (el) \ 1309 if (el) \
1310 { \ 1310 { \
1311 r.data (ii) = el; \ 1311 r.data (ii) = el; \
1312 r.ridx (ii++) = i; \ 1312 r.ridx (ii++) = i; \
1313 } \ 1313 } \
1314 } \ 1314 } \
1315 r.cidx (j+1) = ii; \ 1315 r.cidx (j+1) = ii; \
1316 } \ 1316 } \
1317 } \ 1317 } \
1318 } \ 1318 } \
1319 else \ 1319 else \
1320 { \ 1320 { \
1321 if ((m1_nr != 0 || m1_nc != 0) && (m2_nr != 0 || m2_nc != 0)) \ 1321 if ((m1_nr != 0 || m1_nc != 0) && (m2_nr != 0 || m2_nc != 0)) \
1322 err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \ 1322 err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \
1323 } \ 1323 } \
1324 return r; \ 1324 return r; \
1325 } 1325 }
1326 1326
1327 #define SPARSE_SMM_CMP_OPS(M1, Z1, C1, M2, Z2, C2) \ 1327 #define SPARSE_SMM_CMP_OPS(M1, Z1, C1, M2, Z2, C2) \
1328 SPARSE_SMM_CMP_OP (mx_el_lt, <, M1, , M2, ) \ 1328 SPARSE_SMM_CMP_OP (mx_el_lt, <, M1, , M2, ) \
1329 SPARSE_SMM_CMP_OP (mx_el_le, <=, M1, , M2, ) \ 1329 SPARSE_SMM_CMP_OP (mx_el_le, <=, M1, , M2, ) \
1330 SPARSE_SMM_CMP_OP (mx_el_ge, >=, M1, , M2, ) \ 1330 SPARSE_SMM_CMP_OP (mx_el_ge, >=, M1, , M2, ) \
1331 SPARSE_SMM_CMP_OP (mx_el_gt, >, M1, , M2, ) \ 1331 SPARSE_SMM_CMP_OP (mx_el_gt, >, M1, , M2, ) \
1332 SPARSE_SMM_CMP_OP (mx_el_eq, ==, M1, , M2, ) \ 1332 SPARSE_SMM_CMP_OP (mx_el_eq, ==, M1, , M2, ) \
1333 SPARSE_SMM_CMP_OP (mx_el_ne, !=, M1, , M2, ) 1333 SPARSE_SMM_CMP_OP (mx_el_ne, !=, M1, , M2, )
1334 1334
1335 #define SPARSE_SMM_EQNE_OPS(M1, Z1, C1, M2, Z2, C2) \ 1335 #define SPARSE_SMM_EQNE_OPS(M1, Z1, C1, M2, Z2, C2) \
1336 SPARSE_SMM_CMP_OP (mx_el_eq, ==, M1, , M2, ) \ 1336 SPARSE_SMM_CMP_OP (mx_el_eq, ==, M1, , M2, ) \
1337 SPARSE_SMM_CMP_OP (mx_el_ne, !=, M1, , M2, ) 1337 SPARSE_SMM_CMP_OP (mx_el_ne, !=, M1, , M2, )
1338 1338
1339 #define SPARSE_SMM_BOOL_OP(F, OP, M1, M2, LHS_ZERO, RHS_ZERO) \ 1339 #define SPARSE_SMM_BOOL_OP(F, OP, M1, M2, LHS_ZERO, RHS_ZERO) \
1340 SparseBoolMatrix \ 1340 SparseBoolMatrix \
1341 F (const M1& m1, const M2& m2) \ 1341 F (const M1& m1, const M2& m2) \
1342 { \ 1342 { \
1343 SparseBoolMatrix r; \ 1343 SparseBoolMatrix r; \
1344 \ 1344 \
1345 octave_idx_type m1_nr = m1.rows (); \ 1345 octave_idx_type m1_nr = m1.rows (); \
1346 octave_idx_type m1_nc = m1.cols (); \ 1346 octave_idx_type m1_nc = m1.cols (); \
1347 \ 1347 \
1348 octave_idx_type m2_nr = m2.rows (); \ 1348 octave_idx_type m2_nr = m2.rows (); \
1349 octave_idx_type m2_nc = m2.cols (); \ 1349 octave_idx_type m2_nc = m2.cols (); \
1350 \ 1350 \
1351 if (m1_nr == 1 && m1_nc == 1) \ 1351 if (m1_nr == 1 && m1_nc == 1) \
1352 r = SparseBoolMatrix (F (m1.elem (0,0), m2)); \ 1352 r = SparseBoolMatrix (F (m1.elem (0,0), m2)); \
1353 else if (m1_nr == m2_nr && m1_nc == m2_nc) \ 1353 else if (m1_nr == m2_nr && m1_nc == m2_nc) \
1354 { \ 1354 { \
1355 if (m1_nr != 0 || m1_nc != 0) \ 1355 if (m1_nr != 0 || m1_nc != 0) \
1356 { \ 1356 { \
1357 /* Count num of nonzero elements */ \ 1357 /* Count num of nonzero elements */ \
1358 octave_idx_type nel = 0; \ 1358 octave_idx_type nel = 0; \
1359 for (octave_idx_type j = 0; j < m1_nc; j++) \ 1359 for (octave_idx_type j = 0; j < m1_nc; j++) \
1360 for (octave_idx_type i = 0; i < m1_nr; i++) \ 1360 for (octave_idx_type i = 0; i < m1_nr; i++) \
1361 if ((m1.elem (i, j) != LHS_ZERO) \ 1361 if ((m1.elem (i, j) != LHS_ZERO) \
1362 OP (m2.elem (i, j) != RHS_ZERO)) \ 1362 OP (m2.elem (i, j) != RHS_ZERO)) \
1363 nel++; \ 1363 nel++; \
1364 \ 1364 \
1365 r = SparseBoolMatrix (m1_nr, m1_nc, nel); \ 1365 r = SparseBoolMatrix (m1_nr, m1_nc, nel); \
1366 \ 1366 \
1367 octave_idx_type ii = 0; \ 1367 octave_idx_type ii = 0; \
1368 r.cidx (0) = 0; \ 1368 r.cidx (0) = 0; \
1369 for (octave_idx_type j = 0; j < m1_nc; j++) \ 1369 for (octave_idx_type j = 0; j < m1_nc; j++) \
1370 { \ 1370 { \
1371 for (octave_idx_type i = 0; i < m1_nr; i++) \ 1371 for (octave_idx_type i = 0; i < m1_nr; i++) \
1372 { \ 1372 { \
1373 bool el = (m1.elem (i, j) != LHS_ZERO) \ 1373 bool el = (m1.elem (i, j) != LHS_ZERO) \
1374 OP (m2.elem (i, j) != RHS_ZERO); \ 1374 OP (m2.elem (i, j) != RHS_ZERO); \
1375 if (el) \ 1375 if (el) \
1376 { \ 1376 { \
1377 r.data (ii) = el; \ 1377 r.data (ii) = el; \
1378 r.ridx (ii++) = i; \ 1378 r.ridx (ii++) = i; \
1379 } \ 1379 } \
1380 } \ 1380 } \
1381 r.cidx (j+1) = ii; \ 1381 r.cidx (j+1) = ii; \
1382 } \ 1382 } \
1383 } \ 1383 } \
1384 } \ 1384 } \
1385 else \ 1385 else \
1386 { \ 1386 { \
1387 if ((m1_nr != 0 || m1_nc != 0) && (m2_nr != 0 || m2_nc != 0)) \ 1387 if ((m1_nr != 0 || m1_nc != 0) && (m2_nr != 0 || m2_nc != 0)) \
1388 err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \ 1388 err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \
1389 } \ 1389 } \
1390 return r; \ 1390 return r; \
1391 } 1391 }
1392 1392
1393 #define SPARSE_SMM_BOOL_OPS2(M1, M2, LHS_ZERO, RHS_ZERO) \ 1393 #define SPARSE_SMM_BOOL_OPS2(M1, M2, LHS_ZERO, RHS_ZERO) \
1394 SPARSE_SMM_BOOL_OP (mx_el_and, &&, M1, M2, LHS_ZERO, RHS_ZERO) \ 1394 SPARSE_SMM_BOOL_OP (mx_el_and, &&, M1, M2, LHS_ZERO, RHS_ZERO) \
1395 SPARSE_SMM_BOOL_OP (mx_el_or, ||, M1, M2, LHS_ZERO, RHS_ZERO) 1395 SPARSE_SMM_BOOL_OP (mx_el_or, ||, M1, M2, LHS_ZERO, RHS_ZERO)
1396 1396
1397 #define SPARSE_SMM_BOOL_OPS(M1, M2, ZERO) \ 1397 #define SPARSE_SMM_BOOL_OPS(M1, M2, ZERO) \
1398 SPARSE_SMM_BOOL_OPS2(M1, M2, ZERO, ZERO) 1398 SPARSE_SMM_BOOL_OPS2(M1, M2, ZERO, ZERO)
1399 1399
1400 // Avoid some code duplication. Maybe we should use templates. 1400 // Avoid some code duplication. Maybe we should use templates.
1401 1401
1402 #define SPARSE_CUMSUM(RET_TYPE, ELT_TYPE, FCN) \ 1402 #define SPARSE_CUMSUM(RET_TYPE, ELT_TYPE, FCN) \
1403 \ 1403 \
1404 octave_idx_type nr = rows (); \ 1404 octave_idx_type nr = rows (); \
1405 octave_idx_type nc = cols (); \ 1405 octave_idx_type nc = cols (); \
1406 \ 1406 \
1407 RET_TYPE retval; \ 1407 RET_TYPE retval; \
1408 \ 1408 \
1409 if (nr > 0 && nc > 0) \ 1409 if (nr > 0 && nc > 0) \
1410 { \ 1410 { \
1411 if ((nr == 1 && dim == -1) || dim == 1) \ 1411 if ((nr == 1 && dim == -1) || dim == 1) \
1412 /* Ugly!! Is there a better way? */ \ 1412 /* Ugly!! Is there a better way? */ \
1413 retval = transpose (). FCN (0) .transpose (); \ 1413 retval = transpose (). FCN (0) .transpose (); \
1414 else \ 1414 else \
1415 { \ 1415 { \
1416 octave_idx_type nel = 0; \ 1416 octave_idx_type nel = 0; \
1417 for (octave_idx_type i = 0; i < nc; i++) \ 1417 for (octave_idx_type i = 0; i < nc; i++) \
1418 { \ 1418 { \
1419 ELT_TYPE t = ELT_TYPE (); \ 1419 ELT_TYPE t = ELT_TYPE (); \
1420 for (octave_idx_type j = cidx (i); j < cidx (i+1); j++) \ 1420 for (octave_idx_type j = cidx (i); j < cidx (i+1); j++) \
1421 { \ 1421 { \
1422 t += data (j); \ 1422 t += data (j); \
1423 if (t != ELT_TYPE ()) \ 1423 if (t != ELT_TYPE ()) \
1424 { \ 1424 { \
1425 if (j == cidx (i+1) - 1) \ 1425 if (j == cidx (i+1) - 1) \
1426 nel += nr - ridx (j); \ 1426 nel += nr - ridx (j); \
1427 else \ 1427 else \
1428 nel += ridx (j+1) - ridx (j); \ 1428 nel += ridx (j+1) - ridx (j); \
1429 } \ 1429 } \
1430 } \ 1430 } \
1431 } \ 1431 } \
1432 retval = RET_TYPE (nr, nc, nel); \ 1432 retval = RET_TYPE (nr, nc, nel); \
1433 retval.cidx (0) = 0; \ 1433 retval.cidx (0) = 0; \
1434 octave_idx_type ii = 0; \ 1434 octave_idx_type ii = 0; \
1435 for (octave_idx_type i = 0; i < nc; i++) \ 1435 for (octave_idx_type i = 0; i < nc; i++) \
1436 { \ 1436 { \
1437 ELT_TYPE t = ELT_TYPE (); \ 1437 ELT_TYPE t = ELT_TYPE (); \
1438 for (octave_idx_type j = cidx (i); j < cidx (i+1); j++) \ 1438 for (octave_idx_type j = cidx (i); j < cidx (i+1); j++) \
1439 { \ 1439 { \
1440 t += data (j); \ 1440 t += data (j); \
1441 if (t != ELT_TYPE ()) \ 1441 if (t != ELT_TYPE ()) \
1442 { \ 1442 { \
1443 if (j == cidx (i+1) - 1) \ 1443 if (j == cidx (i+1) - 1) \
1444 { \ 1444 { \
1445 for (octave_idx_type k = ridx (j); k < nr; k++) \ 1445 for (octave_idx_type k = ridx (j); k < nr; k++) \
1446 { \ 1446 { \
1447 retval.data (ii) = t; \ 1447 retval.data (ii) = t; \
1448 retval.ridx (ii++) = k; \ 1448 retval.ridx (ii++) = k; \
1449 } \ 1449 } \
1450 } \ 1450 } \
1451 else \ 1451 else \
1452 { \ 1452 { \
1453 for (octave_idx_type k = ridx (j); k < ridx (j+1); k++) \ 1453 for (octave_idx_type k = ridx (j); k < ridx (j+1); k++) \
1454 { \ 1454 { \
1455 retval.data (ii) = t; \ 1455 retval.data (ii) = t; \
1456 retval.ridx (ii++) = k; \ 1456 retval.ridx (ii++) = k; \
1457 } \ 1457 } \
1458 } \ 1458 } \
1459 } \ 1459 } \
1460 } \ 1460 } \
1461 retval.cidx (i+1) = ii; \ 1461 retval.cidx (i+1) = ii; \
1462 } \ 1462 } \
1463 } \ 1463 } \
1464 } \ 1464 } \
1465 else \ 1465 else \
1466 retval = RET_TYPE (nr,nc); \ 1466 retval = RET_TYPE (nr,nc); \
1467 \ 1467 \
1468 return retval 1468 return retval
1469 1469
1470 #define SPARSE_CUMPROD(RET_TYPE, ELT_TYPE, FCN) \ 1470 #define SPARSE_CUMPROD(RET_TYPE, ELT_TYPE, FCN) \
1471 \ 1471 \
1472 octave_idx_type nr = rows (); \ 1472 octave_idx_type nr = rows (); \
1473 octave_idx_type nc = cols (); \ 1473 octave_idx_type nc = cols (); \
1474 \ 1474 \
1475 RET_TYPE retval; \ 1475 RET_TYPE retval; \
1476 \ 1476 \
1477 if (nr > 0 && nc > 0) \ 1477 if (nr > 0 && nc > 0) \
1478 { \ 1478 { \
1479 if ((nr == 1 && dim == -1) || dim == 1) \ 1479 if ((nr == 1 && dim == -1) || dim == 1) \
1480 /* Ugly!! Is there a better way? */ \ 1480 /* Ugly!! Is there a better way? */ \
1481 retval = transpose (). FCN (0) .transpose (); \ 1481 retval = transpose (). FCN (0) .transpose (); \
1482 else \ 1482 else \
1483 { \ 1483 { \
1484 octave_idx_type nel = 0; \ 1484 octave_idx_type nel = 0; \
1485 for (octave_idx_type i = 0; i < nc; i++) \ 1485 for (octave_idx_type i = 0; i < nc; i++) \
1486 { \ 1486 { \
1487 octave_idx_type jj = 0; \ 1487 octave_idx_type jj = 0; \
1488 for (octave_idx_type j = cidx (i); j < cidx (i+1); j++) \ 1488 for (octave_idx_type j = cidx (i); j < cidx (i+1); j++) \
1489 { \ 1489 { \
1490 if (jj == ridx (j)) \ 1490 if (jj == ridx (j)) \
1491 { \ 1491 { \
1492 nel++; \ 1492 nel++; \
1493 jj++; \ 1493 jj++; \
1494 } \ 1494 } \
1495 else \ 1495 else \
1496 break; \ 1496 break; \
1497 } \ 1497 } \
1498 } \ 1498 } \
1499 retval = RET_TYPE (nr, nc, nel); \ 1499 retval = RET_TYPE (nr, nc, nel); \
1500 retval.cidx (0) = 0; \ 1500 retval.cidx (0) = 0; \
1501 octave_idx_type ii = 0; \ 1501 octave_idx_type ii = 0; \
1502 for (octave_idx_type i = 0; i < nc; i++) \ 1502 for (octave_idx_type i = 0; i < nc; i++) \
1503 { \ 1503 { \
1504 ELT_TYPE t = ELT_TYPE (1.); \ 1504 ELT_TYPE t = ELT_TYPE (1.); \
1505 octave_idx_type jj = 0; \ 1505 octave_idx_type jj = 0; \
1506 for (octave_idx_type j = cidx (i); j < cidx (i+1); j++) \ 1506 for (octave_idx_type j = cidx (i); j < cidx (i+1); j++) \
1507 { \ 1507 { \
1508 if (jj == ridx (j)) \ 1508 if (jj == ridx (j)) \
1509 { \ 1509 { \
1510 t *= data (j); \ 1510 t *= data (j); \
1511 retval.data (ii) = t; \ 1511 retval.data (ii) = t; \
1512 retval.ridx (ii++) = jj++; \ 1512 retval.ridx (ii++) = jj++; \
1513 } \ 1513 } \
1514 else \ 1514 else \
1515 break; \ 1515 break; \
1516 } \ 1516 } \
1517 retval.cidx (i+1) = ii; \ 1517 retval.cidx (i+1) = ii; \
1518 } \ 1518 } \
1519 } \ 1519 } \
1520 } \ 1520 } \
1521 else \ 1521 else \
1522 retval = RET_TYPE (nr,nc); \ 1522 retval = RET_TYPE (nr,nc); \
1523 \ 1523 \
1524 return retval 1524 return retval
1525 1525
1526 #define SPARSE_BASE_REDUCTION_OP(RET_TYPE, EL_TYPE, ROW_EXPR, COL_EXPR, \ 1526 #define SPARSE_BASE_REDUCTION_OP(RET_TYPE, EL_TYPE, ROW_EXPR, COL_EXPR, \
1527 INIT_VAL, MT_RESULT) \ 1527 INIT_VAL, MT_RESULT) \
1528 \ 1528 \
1529 octave_idx_type nr = rows (); \ 1529 octave_idx_type nr = rows (); \
1530 octave_idx_type nc = cols (); \ 1530 octave_idx_type nc = cols (); \
1531 \ 1531 \
1532 RET_TYPE retval; \ 1532 RET_TYPE retval; \
1533 \ 1533 \
1534 if (nr > 0 && nc > 0) \ 1534 if (nr > 0 && nc > 0) \
1535 { \ 1535 { \
1536 if ((nr == 1 && dim == -1) || dim == 1) \ 1536 if ((nr == 1 && dim == -1) || dim == 1) \
1537 { \ 1537 { \
1538 /* Define j here to allow fancy definition for prod method */ \ 1538 /* Define j here to allow fancy definition for prod method */ \
1539 octave_idx_type j = 0; \ 1539 octave_idx_type j = 0; \
1540 OCTAVE_LOCAL_BUFFER (EL_TYPE, tmp, nr); \ 1540 OCTAVE_LOCAL_BUFFER (EL_TYPE, tmp, nr); \
1541 \ 1541 \
1542 for (octave_idx_type i = 0; i < nr; i++) \ 1542 for (octave_idx_type i = 0; i < nr; i++) \
1543 tmp[i] = INIT_VAL; \ 1543 tmp[i] = INIT_VAL; \
1544 for (j = 0; j < nc; j++) \ 1544 for (j = 0; j < nc; j++) \
1545 { \ 1545 { \
1546 for (octave_idx_type i = cidx (j); i < cidx (j + 1); i++) \ 1546 for (octave_idx_type i = cidx (j); i < cidx (j + 1); i++) \
1547 { \ 1547 { \
1548 ROW_EXPR; \ 1548 ROW_EXPR; \
1549 } \ 1549 } \
1550 } \ 1550 } \
1551 octave_idx_type nel = 0; \ 1551 octave_idx_type nel = 0; \
1552 for (octave_idx_type i = 0; i < nr; i++) \ 1552 for (octave_idx_type i = 0; i < nr; i++) \
1553 if (tmp[i] != EL_TYPE ()) \ 1553 if (tmp[i] != EL_TYPE ()) \
1554 nel++; \ 1554 nel++; \
1555 retval = RET_TYPE (nr, static_cast<octave_idx_type> (1), nel); \ 1555 retval = RET_TYPE (nr, static_cast<octave_idx_type> (1), nel); \
1556 retval.cidx (0) = 0; \ 1556 retval.cidx (0) = 0; \
1557 retval.cidx (1) = nel; \ 1557 retval.cidx (1) = nel; \
1558 nel = 0; \ 1558 nel = 0; \
1559 for (octave_idx_type i = 0; i < nr; i++) \ 1559 for (octave_idx_type i = 0; i < nr; i++) \
1560 if (tmp[i] != EL_TYPE ()) \ 1560 if (tmp[i] != EL_TYPE ()) \
1561 { \ 1561 { \
1562 retval.data (nel) = tmp[i]; \ 1562 retval.data (nel) = tmp[i]; \
1563 retval.ridx (nel++) = i; \ 1563 retval.ridx (nel++) = i; \
1564 } \ 1564 } \
1565 } \ 1565 } \
1566 else \ 1566 else \
1567 { \ 1567 { \
1568 OCTAVE_LOCAL_BUFFER (EL_TYPE, tmp, nc); \ 1568 OCTAVE_LOCAL_BUFFER (EL_TYPE, tmp, nc); \
1569 \ 1569 \
1570 for (octave_idx_type j = 0; j < nc; j++) \ 1570 for (octave_idx_type j = 0; j < nc; j++) \
1571 { \ 1571 { \
1572 tmp[j] = INIT_VAL; \ 1572 tmp[j] = INIT_VAL; \
1573 for (octave_idx_type i = cidx (j); i < cidx (j + 1); i++) \ 1573 for (octave_idx_type i = cidx (j); i < cidx (j + 1); i++) \
1574 { \ 1574 { \
1575 COL_EXPR; \ 1575 COL_EXPR; \
1576 } \ 1576 } \
1577 } \ 1577 } \
1578 octave_idx_type nel = 0; \ 1578 octave_idx_type nel = 0; \
1579 for (octave_idx_type i = 0; i < nc; i++) \ 1579 for (octave_idx_type i = 0; i < nc; i++) \
1580 if (tmp[i] != EL_TYPE ()) \ 1580 if (tmp[i] != EL_TYPE ()) \
1581 nel++; \ 1581 nel++; \
1582 retval = RET_TYPE (static_cast<octave_idx_type> (1), nc, nel); \ 1582 retval = RET_TYPE (static_cast<octave_idx_type> (1), nc, nel); \
1583 retval.cidx (0) = 0; \ 1583 retval.cidx (0) = 0; \
1584 nel = 0; \ 1584 nel = 0; \
1585 for (octave_idx_type i = 0; i < nc; i++) \ 1585 for (octave_idx_type i = 0; i < nc; i++) \
1586 if (tmp[i] != EL_TYPE ()) \ 1586 if (tmp[i] != EL_TYPE ()) \
1587 { \ 1587 { \
1588 retval.data (nel) = tmp[i]; \ 1588 retval.data (nel) = tmp[i]; \
1589 retval.ridx (nel++) = 0; \ 1589 retval.ridx (nel++) = 0; \
1590 retval.cidx (i+1) = retval.cidx (i) + 1; \ 1590 retval.cidx (i+1) = retval.cidx (i) + 1; \
1591 } \ 1591 } \
1592 else \ 1592 else \
1593 retval.cidx (i+1) = retval.cidx (i); \ 1593 retval.cidx (i+1) = retval.cidx (i); \
1594 } \ 1594 } \
1595 } \ 1595 } \
1596 else if (nc == 0 && (nr == 0 || (nr == 1 && dim == -1))) \ 1596 else if (nc == 0 && (nr == 0 || (nr == 1 && dim == -1))) \
1597 { \ 1597 { \
1598 if (MT_RESULT) \ 1598 if (MT_RESULT) \
1599 { \ 1599 { \
1600 retval = RET_TYPE (static_cast<octave_idx_type> (1), \ 1600 retval = RET_TYPE (static_cast<octave_idx_type> (1), \
1601 static_cast<octave_idx_type> (1), \ 1601 static_cast<octave_idx_type> (1), \
1602 static_cast<octave_idx_type> (1)); \ 1602 static_cast<octave_idx_type> (1)); \
1603 retval.cidx (0) = 0; \ 1603 retval.cidx (0) = 0; \
1604 retval.cidx (1) = 1; \ 1604 retval.cidx (1) = 1; \
1605 retval.ridx (0) = 0; \ 1605 retval.ridx (0) = 0; \
1606 retval.data (0) = MT_RESULT; \ 1606 retval.data (0) = MT_RESULT; \
1607 } \ 1607 } \
1608 else \ 1608 else \
1609 retval = RET_TYPE (static_cast<octave_idx_type> (1), \ 1609 retval = RET_TYPE (static_cast<octave_idx_type> (1), \
1610 static_cast<octave_idx_type> (1), \ 1610 static_cast<octave_idx_type> (1), \
1611 static_cast<octave_idx_type> (0)); \ 1611 static_cast<octave_idx_type> (0)); \
1612 } \ 1612 } \
1613 else if (nr == 0 && (dim == 0 || dim == -1)) \ 1613 else if (nr == 0 && (dim == 0 || dim == -1)) \
1614 { \ 1614 { \
1615 if (MT_RESULT) \ 1615 if (MT_RESULT) \
1616 { \ 1616 { \
1617 retval = RET_TYPE (static_cast<octave_idx_type> (1), nc, nc); \ 1617 retval = RET_TYPE (static_cast<octave_idx_type> (1), nc, nc); \
1618 retval.cidx (0) = 0; \ 1618 retval.cidx (0) = 0; \
1619 for (octave_idx_type i = 0; i < nc ; i++) \ 1619 for (octave_idx_type i = 0; i < nc ; i++) \
1620 { \ 1620 { \
1621 retval.ridx (i) = 0; \ 1621 retval.ridx (i) = 0; \
1622 retval.cidx (i+1) = i+1; \ 1622 retval.cidx (i+1) = i+1; \
1623 retval.data (i) = MT_RESULT; \ 1623 retval.data (i) = MT_RESULT; \
1624 } \ 1624 } \
1625 } \ 1625 } \
1626 else \ 1626 else \
1627 retval = RET_TYPE (static_cast<octave_idx_type> (1), nc, \ 1627 retval = RET_TYPE (static_cast<octave_idx_type> (1), nc, \
1628 static_cast<octave_idx_type> (0)); \ 1628 static_cast<octave_idx_type> (0)); \
1629 } \ 1629 } \
1630 else if (nc == 0 && dim == 1) \ 1630 else if (nc == 0 && dim == 1) \
1631 { \ 1631 { \
1632 if (MT_RESULT) \ 1632 if (MT_RESULT) \
1633 { \ 1633 { \
1634 retval = RET_TYPE (nr, static_cast<octave_idx_type> (1), nr); \ 1634 retval = RET_TYPE (nr, static_cast<octave_idx_type> (1), nr); \
1635 retval.cidx (0) = 0; \ 1635 retval.cidx (0) = 0; \
1636 retval.cidx (1) = nr; \ 1636 retval.cidx (1) = nr; \
1637 for (octave_idx_type i = 0; i < nr; i++) \ 1637 for (octave_idx_type i = 0; i < nr; i++) \
1638 { \ 1638 { \
1639 retval.ridx (i) = i; \ 1639 retval.ridx (i) = i; \
1640 retval.data (i) = MT_RESULT; \ 1640 retval.data (i) = MT_RESULT; \
1641 } \ 1641 } \
1642 } \ 1642 } \
1643 else \ 1643 else \
1644 retval = RET_TYPE (nr, static_cast<octave_idx_type> (1), \ 1644 retval = RET_TYPE (nr, static_cast<octave_idx_type> (1), \
1645 static_cast<octave_idx_type> (0)); \ 1645 static_cast<octave_idx_type> (0)); \
1646 } \ 1646 } \
1647 else \ 1647 else \
1648 retval.resize (nr > 0, nc > 0); \ 1648 retval.resize (nr > 0, nc > 0); \
1649 \ 1649 \
1650 return retval 1650 return retval
1651 1651
1652 #define SPARSE_REDUCTION_OP_ROW_EXPR(OP) \ 1652 #define SPARSE_REDUCTION_OP_ROW_EXPR(OP) \
1653 tmp[ridx (i)] OP data (i) 1653 tmp[ridx (i)] OP data (i)
1654 1654
1655 #define SPARSE_REDUCTION_OP_COL_EXPR(OP) \ 1655 #define SPARSE_REDUCTION_OP_COL_EXPR(OP) \
1656 tmp[j] OP data (i) 1656 tmp[j] OP data (i)
1657 1657
1658 #define SPARSE_REDUCTION_OP(RET_TYPE, EL_TYPE, OP, INIT_VAL, MT_RESULT) \ 1658 #define SPARSE_REDUCTION_OP(RET_TYPE, EL_TYPE, OP, INIT_VAL, MT_RESULT) \
1659 SPARSE_BASE_REDUCTION_OP (RET_TYPE, EL_TYPE, \ 1659 SPARSE_BASE_REDUCTION_OP (RET_TYPE, EL_TYPE, \
1660 SPARSE_REDUCTION_OP_ROW_EXPR (OP), \ 1660 SPARSE_REDUCTION_OP_ROW_EXPR (OP), \
1661 SPARSE_REDUCTION_OP_COL_EXPR (OP), \ 1661 SPARSE_REDUCTION_OP_COL_EXPR (OP), \
1662 INIT_VAL, MT_RESULT) 1662 INIT_VAL, MT_RESULT)
1663 1663
1664 // Don't break from this loop if the test succeeds because 1664 // Don't break from this loop if the test succeeds because
1665 // we are looping over the rows and not the columns in the inner loop. 1665 // we are looping over the rows and not the columns in the inner loop.
1666 #define SPARSE_ANY_ALL_OP_ROW_CODE(TEST_OP, TEST_TRUE_VAL) \ 1666 #define SPARSE_ANY_ALL_OP_ROW_CODE(TEST_OP, TEST_TRUE_VAL) \
1667 if (data (i) TEST_OP 0.0) \ 1667 if (data (i) TEST_OP 0.0) \
1668 tmp[ridx (i)] = TEST_TRUE_VAL; 1668 tmp[ridx (i)] = TEST_TRUE_VAL;
1669 1669
1670 #define SPARSE_ANY_ALL_OP_COL_CODE(TEST_OP, TEST_TRUE_VAL) \ 1670 #define SPARSE_ANY_ALL_OP_COL_CODE(TEST_OP, TEST_TRUE_VAL) \
1671 if (data (i) TEST_OP 0.0) \ 1671 if (data (i) TEST_OP 0.0) \
1672 { \ 1672 { \
1673 tmp[j] = TEST_TRUE_VAL; \ 1673 tmp[j] = TEST_TRUE_VAL; \
1674 break; \ 1674 break; \
1675 } 1675 }
1676 1676
1677 #define SPARSE_ANY_ALL_OP(DIM, INIT_VAL, MT_RESULT, TEST_OP, TEST_TRUE_VAL) \ 1677 #define SPARSE_ANY_ALL_OP(DIM, INIT_VAL, MT_RESULT, TEST_OP, TEST_TRUE_VAL) \
1678 SPARSE_BASE_REDUCTION_OP (SparseBoolMatrix, char, \ 1678 SPARSE_BASE_REDUCTION_OP (SparseBoolMatrix, char, \
1679 SPARSE_ANY_ALL_OP_ROW_CODE (TEST_OP, TEST_TRUE_VAL), \ 1679 SPARSE_ANY_ALL_OP_ROW_CODE (TEST_OP, TEST_TRUE_VAL), \
1680 SPARSE_ANY_ALL_OP_COL_CODE (TEST_OP, TEST_TRUE_VAL), \ 1680 SPARSE_ANY_ALL_OP_COL_CODE (TEST_OP, TEST_TRUE_VAL), \
1681 INIT_VAL, MT_RESULT) 1681 INIT_VAL, MT_RESULT)
1682 1682
1683 #define SPARSE_ALL_OP(DIM) \ 1683 #define SPARSE_ALL_OP(DIM) \
1684 if ((rows () == 1 && dim == -1) || dim == 1) \ 1684 if ((rows () == 1 && dim == -1) || dim == 1) \
1685 return transpose (). all (0). transpose (); \ 1685 return transpose (). all (0). transpose (); \
1686 else \ 1686 else \
1687 { \ 1687 { \
1688 SPARSE_ANY_ALL_OP (DIM, (cidx (j+1) - cidx (j) < nr ? false : true), \ 1688 SPARSE_ANY_ALL_OP (DIM, (cidx (j+1) - cidx (j) < nr ? false : true), \
1689 true, ==, false); \ 1689 true, ==, false); \
1690 } 1690 }
1691 1691
1692 #define SPARSE_ANY_OP(DIM) SPARSE_ANY_ALL_OP (DIM, false, false, !=, true) 1692 #define SPARSE_ANY_OP(DIM) SPARSE_ANY_ALL_OP (DIM, false, false, !=, true)
1693 1693
1694 #define SPARSE_SPARSE_MUL(RET_TYPE, RET_EL_TYPE, EL_TYPE) \ 1694 #define SPARSE_SPARSE_MUL(RET_TYPE, RET_EL_TYPE, EL_TYPE) \
1695 octave_idx_type nr = m.rows (); \ 1695 octave_idx_type nr = m.rows (); \
1696 octave_idx_type nc = m.cols (); \ 1696 octave_idx_type nc = m.cols (); \
1697 \ 1697 \
1698 octave_idx_type a_nr = a.rows (); \ 1698 octave_idx_type a_nr = a.rows (); \
1699 octave_idx_type a_nc = a.cols (); \ 1699 octave_idx_type a_nc = a.cols (); \
1700 \ 1700 \
1701 if (nr == 1 && nc == 1) \ 1701 if (nr == 1 && nc == 1) \
1702 { \ 1702 { \
1703 RET_EL_TYPE s = m.elem (0,0); \ 1703 RET_EL_TYPE s = m.elem (0,0); \
1704 octave_idx_type nz = a.nnz (); \ 1704 octave_idx_type nz = a.nnz (); \
1705 RET_TYPE r (a_nr, a_nc, nz); \ 1705 RET_TYPE r (a_nr, a_nc, nz); \
1706 \ 1706 \
1707 for (octave_idx_type i = 0; i < nz; i++) \ 1707 for (octave_idx_type i = 0; i < nz; i++) \
1708 { \ 1708 { \
1709 octave_quit (); \ 1709 octave_quit (); \
1710 r.data (i) = s * a.data (i); \ 1710 r.data (i) = s * a.data (i); \
1711 r.ridx (i) = a.ridx (i); \ 1711 r.ridx (i) = a.ridx (i); \
1712 } \ 1712 } \
1713 for (octave_idx_type i = 0; i < a_nc + 1; i++) \ 1713 for (octave_idx_type i = 0; i < a_nc + 1; i++) \
1714 { \ 1714 { \
1715 octave_quit (); \ 1715 octave_quit (); \
1716 r.cidx (i) = a.cidx (i); \ 1716 r.cidx (i) = a.cidx (i); \
1717 } \ 1717 } \
1718 \ 1718 \
1719 r.maybe_compress (true); \ 1719 r.maybe_compress (true); \
1720 return r; \ 1720 return r; \
1721 } \ 1721 } \
1722 else if (a_nr == 1 && a_nc == 1) \ 1722 else if (a_nr == 1 && a_nc == 1) \
1723 { \ 1723 { \
1724 RET_EL_TYPE s = a.elem (0,0); \ 1724 RET_EL_TYPE s = a.elem (0,0); \
1725 octave_idx_type nz = m.nnz (); \ 1725 octave_idx_type nz = m.nnz (); \
1726 RET_TYPE r (nr, nc, nz); \ 1726 RET_TYPE r (nr, nc, nz); \
1727 \ 1727 \
1728 for (octave_idx_type i = 0; i < nz; i++) \ 1728 for (octave_idx_type i = 0; i < nz; i++) \
1729 { \ 1729 { \
1730 octave_quit (); \ 1730 octave_quit (); \
1731 r.data (i) = m.data (i) * s; \ 1731 r.data (i) = m.data (i) * s; \
1732 r.ridx (i) = m.ridx (i); \ 1732 r.ridx (i) = m.ridx (i); \
1733 } \ 1733 } \
1734 for (octave_idx_type i = 0; i < nc + 1; i++) \ 1734 for (octave_idx_type i = 0; i < nc + 1; i++) \
1735 { \ 1735 { \
1736 octave_quit (); \ 1736 octave_quit (); \
1737 r.cidx (i) = m.cidx (i); \ 1737 r.cidx (i) = m.cidx (i); \
1738 } \ 1738 } \
1739 \ 1739 \
1740 r.maybe_compress (true); \ 1740 r.maybe_compress (true); \
1741 return r; \ 1741 return r; \
1742 } \ 1742 } \
1743 else if (nc != a_nr) \ 1743 else if (nc != a_nr) \
1744 err_nonconformant ("operator *", nr, nc, a_nr, a_nc); \ 1744 err_nonconformant ("operator *", nr, nc, a_nr, a_nc); \
1745 else \ 1745 else \
1746 { \ 1746 { \
1747 OCTAVE_LOCAL_BUFFER (octave_idx_type, w, nr); \ 1747 OCTAVE_LOCAL_BUFFER (octave_idx_type, w, nr); \
1748 RET_TYPE retval (nr, a_nc, static_cast<octave_idx_type> (0)); \ 1748 RET_TYPE retval (nr, a_nc, static_cast<octave_idx_type> (0)); \
1749 for (octave_idx_type i = 0; i < nr; i++) \ 1749 for (octave_idx_type i = 0; i < nr; i++) \
1750 w[i] = 0; \ 1750 w[i] = 0; \
1751 retval.xcidx (0) = 0; \ 1751 retval.xcidx (0) = 0; \
1752 \ 1752 \
1753 octave_idx_type nel = 0; \ 1753 octave_idx_type nel = 0; \
1754 \ 1754 \
1755 for (octave_idx_type i = 0; i < a_nc; i++) \ 1755 for (octave_idx_type i = 0; i < a_nc; i++) \
1756 { \ 1756 { \
1757 for (octave_idx_type j = a.cidx (i); j < a.cidx (i+1); j++) \ 1757 for (octave_idx_type j = a.cidx (i); j < a.cidx (i+1); j++) \
1758 { \ 1758 { \
1759 octave_idx_type col = a.ridx (j); \ 1759 octave_idx_type col = a.ridx (j); \
1760 for (octave_idx_type k = m.cidx (col) ; k < m.cidx (col+1); k++) \ 1760 for (octave_idx_type k = m.cidx (col) ; k < m.cidx (col+1); k++) \
1761 { \ 1761 { \
1762 if (w[m.ridx (k)] < i + 1) \ 1762 if (w[m.ridx (k)] < i + 1) \
1763 { \ 1763 { \
1764 w[m.ridx (k)] = i + 1; \ 1764 w[m.ridx (k)] = i + 1; \
1765 nel++; \ 1765 nel++; \
1766 } \ 1766 } \
1767 octave_quit (); \ 1767 octave_quit (); \
1768 } \ 1768 } \
1769 } \ 1769 } \
1770 retval.xcidx (i+1) = nel; \ 1770 retval.xcidx (i+1) = nel; \
1771 } \ 1771 } \
1772 \ 1772 \
1773 if (nel == 0) \ 1773 if (nel == 0) \
1774 return RET_TYPE (nr, a_nc); \ 1774 return RET_TYPE (nr, a_nc); \
1775 else \ 1775 else \
1776 { \ 1776 { \
1777 for (octave_idx_type i = 0; i < nr; i++) \ 1777 for (octave_idx_type i = 0; i < nr; i++) \
1778 w[i] = 0; \ 1778 w[i] = 0; \
1779 \ 1779 \
1780 OCTAVE_LOCAL_BUFFER (RET_EL_TYPE, Xcol, nr); \ 1780 OCTAVE_LOCAL_BUFFER (RET_EL_TYPE, Xcol, nr); \
1781 \ 1781 \
1782 retval.change_capacity (nel); \ 1782 retval.change_capacity (nel); \
1783 /* The optimal break-point as estimated from simulations */ \ 1783 /* The optimal break-point as estimated from simulations */ \
1784 /* Note that Mergesort is O(nz log(nz)) while searching all */ \ 1784 /* Note that Mergesort is O(nz log(nz)) while searching all */ \
1785 /* values is O(nr), where nz here is nonzero per row of */ \ 1785 /* values is O(nr), where nz here is nonzero per row of */ \
1786 /* length nr. The test itself was then derived from the */ \ 1786 /* length nr. The test itself was then derived from the */ \
1787 /* simulation with random square matrices and the observation */ \ 1787 /* simulation with random square matrices and the observation */ \
1788 /* of the number of nonzero elements in the output matrix */ \ 1788 /* of the number of nonzero elements in the output matrix */ \
1789 /* it was found that the breakpoints were */ \ 1789 /* it was found that the breakpoints were */ \
1790 /* nr: 500 1000 2000 5000 10000 */ \ 1790 /* nr: 500 1000 2000 5000 10000 */ \
1791 /* nz: 6 25 97 585 2202 */ \ 1791 /* nz: 6 25 97 585 2202 */ \
1792 /* The below is a simplication of the 'polyfit'-ed parameters */ \ 1792 /* The below is a simplication of the 'polyfit'-ed parameters */ \
1793 /* to these breakpoints */ \ 1793 /* to these breakpoints */ \
1794 octave_idx_type n_per_col = (a_nc > 43000 ? 43000 : \ 1794 octave_idx_type n_per_col = (a_nc > 43000 ? 43000 : \
1795 (a_nc * a_nc) / 43000); \ 1795 (a_nc * a_nc) / 43000); \
1796 octave_idx_type ii = 0; \ 1796 octave_idx_type ii = 0; \
1797 octave_idx_type *ri = retval.xridx (); \ 1797 octave_idx_type *ri = retval.xridx (); \
1798 octave_sort<octave_idx_type> sort; \ 1798 octave_sort<octave_idx_type> sort; \
1799 \ 1799 \
1800 for (octave_idx_type i = 0; i < a_nc ; i++) \ 1800 for (octave_idx_type i = 0; i < a_nc ; i++) \
1801 { \ 1801 { \
1802 if (retval.xcidx (i+1) - retval.xcidx (i) > n_per_col) \ 1802 if (retval.xcidx (i+1) - retval.xcidx (i) > n_per_col) \
1803 { \ 1803 { \
1804 for (octave_idx_type j = a.cidx (i); j < a.cidx (i+1); j++) \ 1804 for (octave_idx_type j = a.cidx (i); j < a.cidx (i+1); j++) \
1805 { \ 1805 { \
1806 octave_idx_type col = a.ridx (j); \ 1806 octave_idx_type col = a.ridx (j); \
1807 EL_TYPE tmpval = a.data (j); \ 1807 EL_TYPE tmpval = a.data (j); \
1808 for (octave_idx_type k = m.cidx (col) ; \ 1808 for (octave_idx_type k = m.cidx (col) ; \
1809 k < m.cidx (col+1); k++) \ 1809 k < m.cidx (col+1); k++) \
1810 { \ 1810 { \
1811 octave_quit (); \ 1811 octave_quit (); \
1812 octave_idx_type row = m.ridx (k); \ 1812 octave_idx_type row = m.ridx (k); \
1813 if (w[row] < i + 1) \ 1813 if (w[row] < i + 1) \
1814 { \ 1814 { \
1815 w[row] = i + 1; \ 1815 w[row] = i + 1; \
1816 Xcol[row] = tmpval * m.data (k); \ 1816 Xcol[row] = tmpval * m.data (k); \
1817 } \ 1817 } \
1818 else \ 1818 else \
1819 Xcol[row] += tmpval * m.data (k); \ 1819 Xcol[row] += tmpval * m.data (k); \
1820 } \ 1820 } \
1821 } \ 1821 } \
1822 for (octave_idx_type k = 0; k < nr; k++) \ 1822 for (octave_idx_type k = 0; k < nr; k++) \
1823 if (w[k] == i + 1) \ 1823 if (w[k] == i + 1) \
1824 { \ 1824 { \
1825 retval.xdata (ii) = Xcol[k]; \ 1825 retval.xdata (ii) = Xcol[k]; \
1826 retval.xridx (ii++) = k; \ 1826 retval.xridx (ii++) = k; \
1827 } \ 1827 } \
1828 } \ 1828 } \
1829 else \ 1829 else \
1830 { \ 1830 { \
1831 for (octave_idx_type j = a.cidx (i); j < a.cidx (i+1); j++) \ 1831 for (octave_idx_type j = a.cidx (i); j < a.cidx (i+1); j++) \
1832 { \ 1832 { \
1833 octave_idx_type col = a.ridx (j); \ 1833 octave_idx_type col = a.ridx (j); \
1834 EL_TYPE tmpval = a.data (j); \ 1834 EL_TYPE tmpval = a.data (j); \
1835 for (octave_idx_type k = m.cidx (col) ; \ 1835 for (octave_idx_type k = m.cidx (col) ; \
1836 k < m.cidx (col+1); k++) \ 1836 k < m.cidx (col+1); k++) \
1837 { \ 1837 { \
1838 octave_quit (); \ 1838 octave_quit (); \
1839 octave_idx_type row = m.ridx (k); \ 1839 octave_idx_type row = m.ridx (k); \
1840 if (w[row] < i + 1) \ 1840 if (w[row] < i + 1) \
1841 { \ 1841 { \
1842 w[row] = i + 1; \ 1842 w[row] = i + 1; \
1843 retval.xridx (ii++) = row;\ 1843 retval.xridx (ii++) = row; \
1844 Xcol[row] = tmpval * m.data (k); \ 1844 Xcol[row] = tmpval * m.data (k); \
1845 } \ 1845 } \
1846 else \ 1846 else \
1847 Xcol[row] += tmpval * m.data (k); \ 1847 Xcol[row] += tmpval * m.data (k); \
1848 } \ 1848 } \
1849 } \ 1849 } \
1850 sort.sort (ri + retval.xcidx (i), ii - retval.xcidx (i)); \ 1850 sort.sort (ri + retval.xcidx (i), ii - retval.xcidx (i)); \
1851 for (octave_idx_type k = retval.xcidx (i); k < ii; k++) \ 1851 for (octave_idx_type k = retval.xcidx (i); k < ii; k++) \
1852 retval.xdata (k) = Xcol[retval.xridx (k)]; \ 1852 retval.xdata (k) = Xcol[retval.xridx (k)]; \
1853 } \ 1853 } \
1854 } \ 1854 } \
1855 retval.maybe_compress (true);\ 1855 retval.maybe_compress (true); \
1856 return retval; \ 1856 return retval; \
1857 } \ 1857 } \
1858 } 1858 }
1859 1859
1860 #define SPARSE_FULL_MUL(RET_TYPE, EL_TYPE, ZERO) \ 1860 #define SPARSE_FULL_MUL(RET_TYPE, EL_TYPE, ZERO) \
1861 octave_idx_type nr = m.rows (); \ 1861 octave_idx_type nr = m.rows (); \
1862 octave_idx_type nc = m.cols (); \ 1862 octave_idx_type nc = m.cols (); \
1863 \ 1863 \
1864 octave_idx_type a_nr = a.rows (); \ 1864 octave_idx_type a_nr = a.rows (); \
1865 octave_idx_type a_nc = a.cols (); \ 1865 octave_idx_type a_nc = a.cols (); \
1866 \ 1866 \
1867 if (nr == 1 && nc == 1) \ 1867 if (nr == 1 && nc == 1) \
1868 { \ 1868 { \
1869 RET_TYPE retval = m.elem (0,0) * a; \ 1869 RET_TYPE retval = m.elem (0,0) * a; \
1870 return retval; \ 1870 return retval; \
1871 } \ 1871 } \
1872 else if (nc != a_nr) \ 1872 else if (nc != a_nr) \
1873 err_nonconformant ("operator *", nr, nc, a_nr, a_nc); \ 1873 err_nonconformant ("operator *", nr, nc, a_nr, a_nc); \
1874 else \ 1874 else \
1875 { \ 1875 { \
1876 RET_TYPE retval (nr, a_nc, ZERO); \ 1876 RET_TYPE retval (nr, a_nc, ZERO); \
1877 \ 1877 \
1878 for (octave_idx_type i = 0; i < a_nc ; i++) \ 1878 for (octave_idx_type i = 0; i < a_nc ; i++) \
1879 { \ 1879 { \
1880 for (octave_idx_type j = 0; j < a_nr; j++) \ 1880 for (octave_idx_type j = 0; j < a_nr; j++) \
1881 { \ 1881 { \
1882 octave_quit (); \ 1882 octave_quit (); \
1883 \ 1883 \
1884 EL_TYPE tmpval = a.elem (j,i); \ 1884 EL_TYPE tmpval = a.elem (j,i); \
1885 for (octave_idx_type k = m.cidx (j) ; k < m.cidx (j+1); k++) \ 1885 for (octave_idx_type k = m.cidx (j) ; k < m.cidx (j+1); k++) \
1886 retval.elem (m.ridx (k),i) += tmpval * m.data (k); \ 1886 retval.elem (m.ridx (k),i) += tmpval * m.data (k); \
1887 } \ 1887 } \
1888 } \ 1888 } \
1889 return retval; \ 1889 return retval; \
1890 } 1890 }
1891 1891
1892 #define SPARSE_FULL_TRANS_MUL(RET_TYPE, EL_TYPE, ZERO, CONJ_OP) \ 1892 #define SPARSE_FULL_TRANS_MUL(RET_TYPE, EL_TYPE, ZERO, CONJ_OP) \
1893 octave_idx_type nr = m.rows (); \ 1893 octave_idx_type nr = m.rows (); \
1894 octave_idx_type nc = m.cols (); \ 1894 octave_idx_type nc = m.cols (); \
1895 \ 1895 \
1896 octave_idx_type a_nr = a.rows (); \ 1896 octave_idx_type a_nr = a.rows (); \
1897 octave_idx_type a_nc = a.cols (); \ 1897 octave_idx_type a_nc = a.cols (); \
1898 \ 1898 \
1899 if (nr == 1 && nc == 1) \ 1899 if (nr == 1 && nc == 1) \
1900 { \ 1900 { \
1901 RET_TYPE retval = CONJ_OP (m.elem (0,0)) * a; \ 1901 RET_TYPE retval = CONJ_OP (m.elem (0,0)) * a; \
1902 return retval; \ 1902 return retval; \
1903 } \ 1903 } \
1904 else if (nr != a_nr) \ 1904 else if (nr != a_nr) \
1905 err_nonconformant ("operator *", nc, nr, a_nr, a_nc); \ 1905 err_nonconformant ("operator *", nc, nr, a_nr, a_nc); \
1906 else \ 1906 else \
1907 { \ 1907 { \
1908 RET_TYPE retval (nc, a_nc); \ 1908 RET_TYPE retval (nc, a_nc); \
1909 \ 1909 \
1910 for (octave_idx_type i = 0; i < a_nc ; i++) \ 1910 for (octave_idx_type i = 0; i < a_nc ; i++) \
1911 { \ 1911 { \
1912 for (octave_idx_type j = 0; j < nc; j++) \ 1912 for (octave_idx_type j = 0; j < nc; j++) \
1913 { \ 1913 { \
1914 octave_quit (); \ 1914 octave_quit (); \
1915 \ 1915 \
1916 EL_TYPE acc = ZERO; \ 1916 EL_TYPE acc = ZERO; \
1917 for (octave_idx_type k = m.cidx (j) ; k < m.cidx (j+1); k++) \ 1917 for (octave_idx_type k = m.cidx (j) ; k < m.cidx (j+1); k++) \
1918 acc += a.elem (m.ridx (k),i) * CONJ_OP (m.data (k)); \ 1918 acc += a.elem (m.ridx (k),i) * CONJ_OP (m.data (k)); \
1919 retval.xelem (j,i) = acc; \ 1919 retval.xelem (j,i) = acc; \
1920 } \ 1920 } \
1921 } \ 1921 } \
1922 return retval; \ 1922 return retval; \
1923 } 1923 }
1924 1924
1925 #define FULL_SPARSE_MUL(RET_TYPE, EL_TYPE, ZERO) \ 1925 #define FULL_SPARSE_MUL(RET_TYPE, EL_TYPE, ZERO) \
1926 octave_idx_type nr = m.rows (); \ 1926 octave_idx_type nr = m.rows (); \
1927 octave_idx_type nc = m.cols (); \ 1927 octave_idx_type nc = m.cols (); \
1928 \ 1928 \
1929 octave_idx_type a_nr = a.rows (); \ 1929 octave_idx_type a_nr = a.rows (); \
1930 octave_idx_type a_nc = a.cols (); \ 1930 octave_idx_type a_nc = a.cols (); \
1931 \ 1931 \
1932 if (a_nr == 1 && a_nc == 1) \ 1932 if (a_nr == 1 && a_nc == 1) \
1933 { \ 1933 { \
1934 RET_TYPE retval = m * a.elem (0,0); \ 1934 RET_TYPE retval = m * a.elem (0,0); \
1935 return retval; \ 1935 return retval; \
1936 } \ 1936 } \
1937 else if (nc != a_nr) \ 1937 else if (nc != a_nr) \
1938 err_nonconformant ("operator *", nr, nc, a_nr, a_nc); \ 1938 err_nonconformant ("operator *", nr, nc, a_nr, a_nc); \
1939 else \ 1939 else \
1940 { \ 1940 { \
1941 RET_TYPE retval (nr, a_nc, ZERO); \ 1941 RET_TYPE retval (nr, a_nc, ZERO); \
1942 \ 1942 \
1943 for (octave_idx_type i = 0; i < a_nc ; i++) \ 1943 for (octave_idx_type i = 0; i < a_nc ; i++) \
1944 { \ 1944 { \
1945 octave_quit (); \ 1945 octave_quit (); \
1946 for (octave_idx_type j = a.cidx (i); j < a.cidx (i+1); j++) \ 1946 for (octave_idx_type j = a.cidx (i); j < a.cidx (i+1); j++) \
1947 { \ 1947 { \
1948 octave_idx_type col = a.ridx (j); \ 1948 octave_idx_type col = a.ridx (j); \
1949 EL_TYPE tmpval = a.data (j); \ 1949 EL_TYPE tmpval = a.data (j); \
1950 \ 1950 \
1951 for (octave_idx_type k = 0 ; k < nr; k++) \ 1951 for (octave_idx_type k = 0 ; k < nr; k++) \
1952 retval.xelem (k,i) += tmpval * m.elem (k,col); \ 1952 retval.xelem (k,i) += tmpval * m.elem (k,col); \
1953 } \ 1953 } \
1954 } \ 1954 } \
1955 return retval; \ 1955 return retval; \
1956 } 1956 }
1957 1957
1958 #define FULL_SPARSE_MUL_TRANS(RET_TYPE, EL_TYPE, ZERO, CONJ_OP) \ 1958 #define FULL_SPARSE_MUL_TRANS(RET_TYPE, EL_TYPE, ZERO, CONJ_OP) \
1959 octave_idx_type nr = m.rows (); \ 1959 octave_idx_type nr = m.rows (); \
1960 octave_idx_type nc = m.cols (); \ 1960 octave_idx_type nc = m.cols (); \
1961 \ 1961 \
1962 octave_idx_type a_nr = a.rows (); \ 1962 octave_idx_type a_nr = a.rows (); \
1963 octave_idx_type a_nc = a.cols (); \ 1963 octave_idx_type a_nc = a.cols (); \
1964 \ 1964 \
1965 if (a_nr == 1 && a_nc == 1) \ 1965 if (a_nr == 1 && a_nc == 1) \
1966 { \ 1966 { \
1967 RET_TYPE retval = m * CONJ_OP (a.elem (0,0)); \ 1967 RET_TYPE retval = m * CONJ_OP (a.elem (0,0)); \
1968 return retval; \ 1968 return retval; \
1969 } \ 1969 } \
1970 else if (nc != a_nc) \ 1970 else if (nc != a_nc) \
1971 err_nonconformant ("operator *", nr, nc, a_nc, a_nr); \ 1971 err_nonconformant ("operator *", nr, nc, a_nc, a_nr); \
1972 else \ 1972 else \
1973 { \ 1973 { \
1974 RET_TYPE retval (nr, a_nr, ZERO); \ 1974 RET_TYPE retval (nr, a_nr, ZERO); \
1975 \ 1975 \
1976 for (octave_idx_type i = 0; i < a_nc ; i++) \ 1976 for (octave_idx_type i = 0; i < a_nc ; i++) \
1977 { \ 1977 { \
1978 octave_quit (); \ 1978 octave_quit (); \
1979 for (octave_idx_type j = a.cidx (i); j < a.cidx (i+1); j++) \ 1979 for (octave_idx_type j = a.cidx (i); j < a.cidx (i+1); j++) \
1980 { \ 1980 { \
1981 octave_idx_type col = a.ridx (j); \ 1981 octave_idx_type col = a.ridx (j); \
1982 EL_TYPE tmpval = CONJ_OP (a.data (j)); \ 1982 EL_TYPE tmpval = CONJ_OP (a.data (j)); \
1983 for (octave_idx_type k = 0 ; k < nr; k++) \ 1983 for (octave_idx_type k = 0 ; k < nr; k++) \
1984 retval.xelem (k,col) += tmpval * m.elem (k,i); \ 1984 retval.xelem (k,col) += tmpval * m.elem (k,i); \
1985 } \ 1985 } \
1986 } \ 1986 } \
1987 return retval; \ 1987 return retval; \
1988 } 1988 }
1989 1989
1990 #endif 1990 #endif