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