|
5164
|
1 /* |
|
|
2 |
|
|
3 Copyright (C) 2004 David Bateman |
|
|
4 Copyright (C) 1998-2004 Andy Adler |
|
|
5 |
|
|
6 Octave is free software; you can redistribute it and/or modify it |
|
|
7 under the terms of the GNU General Public License as published by the |
|
|
8 Free Software Foundation; either version 2, or (at your option) any |
|
|
9 later version. |
|
|
10 |
|
|
11 Octave is distributed in the hope that it will be useful, but WITHOUT |
|
|
12 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or |
|
|
13 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License |
|
|
14 for more details. |
|
|
15 |
|
|
16 You should have received a copy of the GNU General Public License |
|
|
17 along with this program; see the file COPYING. If not, write to the Free |
|
|
18 Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. |
|
|
19 |
|
|
20 */ |
|
|
21 |
|
|
22 #if !defined (octave_sparse_op_defs_h) |
|
|
23 #define octave_sparse_op_defs_h 1 |
|
|
24 |
|
|
25 #include "Array-util.h" |
|
|
26 |
|
|
27 #define SPARSE_BIN_OP_DECL(R, OP, X, Y) \ |
|
|
28 extern R OP (const X&, const Y&) |
|
|
29 |
|
|
30 #define SPARSE_CMP_OP_DECL(OP, X, Y) \ |
|
|
31 extern SparseBoolMatrix OP (const X&, const Y&) |
|
|
32 |
|
|
33 #define SPARSE_BOOL_OP_DECL(OP, X, Y) \ |
|
|
34 extern SparseBoolMatrix OP (const X&, const Y&) |
|
|
35 |
|
|
36 // matrix by scalar operations. |
|
|
37 |
|
|
38 #define SPARSE_SMS_BIN_OP_DECLS(R1, R2, M, S) \ |
|
|
39 SPARSE_BIN_OP_DECL (R1, operator +, M, S); \ |
|
|
40 SPARSE_BIN_OP_DECL (R1, operator -, M, S); \ |
|
|
41 SPARSE_BIN_OP_DECL (R2, operator *, M, S); \ |
|
|
42 SPARSE_BIN_OP_DECL (R2, operator /, M, S); |
|
|
43 |
|
|
44 #define SPARSE_SMS_BIN_OP_1(R, F, OP, M, S) \ |
|
|
45 R \ |
|
|
46 F (const M& m, const S& s) \ |
|
|
47 { \ |
|
|
48 int nr = m.rows (); \ |
|
|
49 int nc = m.cols (); \ |
|
|
50 \ |
|
|
51 R r (nr, nc, (0.0 OP s)); \ |
|
|
52 \ |
|
|
53 for (int j = 0; j < nc; j++) \ |
|
|
54 for (int i = m.cidx (j); i < m.cidx (j+1); i++) \ |
|
|
55 r.elem (m.ridx (i), j) = m.data (i) OP s; \ |
|
|
56 return r; \ |
|
|
57 } |
|
|
58 |
|
|
59 #define SPARSE_SMS_BIN_OP_2(R, F, OP, M, S) \ |
|
|
60 R \ |
|
|
61 F (const M& m, const S& s) \ |
|
|
62 { \ |
|
|
63 int nr = m.rows (); \ |
|
|
64 int nc = m.cols (); \ |
|
|
65 int nz = m.nnz (); \ |
|
|
66 \ |
|
|
67 R r (nr, nc, nz); \ |
|
|
68 \ |
|
|
69 for (int i = 0; i < nz; i++) \ |
|
|
70 { \ |
|
|
71 r.data(i) = m.data(i) OP s; \ |
|
|
72 r.ridx(i) = m.ridx(i); \ |
|
|
73 } \ |
|
|
74 for (int i = 0; i < nc + 1; i++) \ |
|
|
75 r.cidx(i) = m.cidx(i); \ |
|
|
76 \ |
|
|
77 r.maybe_compress (true); \ |
|
|
78 return r; \ |
|
|
79 } |
|
|
80 |
|
|
81 #define SPARSE_SMS_BIN_OPS(R1, R2, M, S) \ |
|
|
82 SPARSE_SMS_BIN_OP_1 (R1, operator +, +, M, S) \ |
|
|
83 SPARSE_SMS_BIN_OP_1 (R1, operator -, -, M, S) \ |
|
|
84 SPARSE_SMS_BIN_OP_2 (R2, operator *, *, M, S) \ |
|
|
85 SPARSE_SMS_BIN_OP_2 (R2, operator /, /, M, S) |
|
|
86 |
|
|
87 #define SPARSE_SMS_CMP_OP_DECLS(M, S) \ |
|
|
88 SPARSE_CMP_OP_DECL (mx_el_lt, M, S); \ |
|
|
89 SPARSE_CMP_OP_DECL (mx_el_le, M, S); \ |
|
|
90 SPARSE_CMP_OP_DECL (mx_el_ge, M, S); \ |
|
|
91 SPARSE_CMP_OP_DECL (mx_el_gt, M, S); \ |
|
|
92 SPARSE_CMP_OP_DECL (mx_el_eq, M, S); \ |
|
|
93 SPARSE_CMP_OP_DECL (mx_el_ne, M, S); |
|
|
94 |
|
|
95 #define SPARSE_SMS_EQNE_OP_DECLS(M, S) \ |
|
|
96 SPARSE_CMP_OP_DECL (mx_el_eq, M, S); \ |
|
|
97 SPARSE_CMP_OP_DECL (mx_el_ne, M, S); |
|
|
98 |
|
|
99 #define SPARSE_SMS_CMP_OP(F, OP, M, MZ, MC, S, SZ, SC) \ |
|
|
100 SparseBoolMatrix \ |
|
|
101 F (const M& m, const S& s) \ |
|
|
102 { \ |
|
|
103 /* Count num of non-zero elements */ \ |
|
|
104 int nel = 0; \ |
|
|
105 int nz = m.nnz (); \ |
|
|
106 if (MC (MZ) OP SC (s)) \ |
|
|
107 nel += m.numel() - nz; \ |
|
|
108 for (int i = 0; i < nz; i++) \ |
|
|
109 if (MC (m.data (i)) OP SC (s)) \ |
|
|
110 nel++; \ |
|
|
111 \ |
|
|
112 int nr = m.rows (); \ |
|
|
113 int nc = m.cols (); \ |
|
|
114 SparseBoolMatrix r (nr, nc, nel); \ |
|
|
115 \ |
|
|
116 if (nr > 0 && nc > 0) \ |
|
|
117 { \ |
|
|
118 if (MC (MZ) OP SC (s)) \ |
|
|
119 { \ |
|
|
120 int ii = 0; \ |
|
|
121 r.cidx (0) = 0; \ |
|
|
122 for (int j = 0; j < nc; j++) \ |
|
|
123 { \ |
|
|
124 for (int i = 0; i < nr; i++) \ |
|
|
125 { \ |
|
|
126 bool el = MC (m.elem(i, j)) OP SC (s); \ |
|
|
127 if (el) \ |
|
|
128 { \ |
|
|
129 r.data(ii) = el; \ |
|
|
130 r.ridx(ii++) = i; \ |
|
|
131 } \ |
|
|
132 } \ |
|
|
133 r.cidx(j+1) = ii; \ |
|
|
134 } \ |
|
|
135 } \ |
|
|
136 else \ |
|
|
137 { \ |
|
|
138 int ii = 0; \ |
|
|
139 r.cidx (0) = 0; \ |
|
|
140 for (int j = 0; j < nc; j++) \ |
|
|
141 { \ |
|
|
142 for (int i = m.cidx(j); i < m.cidx(j+1); i++) \ |
|
|
143 { \ |
|
|
144 bool el = MC (m.data(i)) OP SC (s); \ |
|
|
145 if (el) \ |
|
|
146 { \ |
|
|
147 r.data(ii) = el; \ |
|
|
148 r.ridx(ii++) = m.ridx(i); \ |
|
|
149 } \ |
|
|
150 } \ |
|
|
151 r.cidx(j+1) = ii; \ |
|
|
152 } \ |
|
|
153 } \ |
|
|
154 } \ |
|
|
155 return r; \ |
|
|
156 } |
|
|
157 |
|
|
158 #define SPARSE_SMS_CMP_OPS(M, MZ, CM, S, SZ, CS) \ |
|
|
159 SPARSE_SMS_CMP_OP (mx_el_lt, <, M, MZ, CM, S, SZ, CS) \ |
|
|
160 SPARSE_SMS_CMP_OP (mx_el_le, <=, M, MZ, CM, S, SZ, CS) \ |
|
|
161 SPARSE_SMS_CMP_OP (mx_el_ge, >=, M, MZ, CM, S, SZ, CS) \ |
|
|
162 SPARSE_SMS_CMP_OP (mx_el_gt, >, M, MZ, CM, S, SZ, CS) \ |
|
|
163 SPARSE_SMS_CMP_OP (mx_el_eq, ==, M, MZ, , S, SZ, ) \ |
|
|
164 SPARSE_SMS_CMP_OP (mx_el_ne, !=, M, MZ, , S, SZ, ) |
|
|
165 |
|
|
166 #define SPARSE_SMS_EQNE_OPS(M, MZ, CM, S, SZ, CS) \ |
|
|
167 SPARSE_SMS_CMP_OP (mx_el_eq, ==, M, MZ, , S, SZ, ) \ |
|
|
168 SPARSE_SMS_CMP_OP (mx_el_ne, !=, M, MZ, , S, SZ, ) |
|
|
169 |
|
|
170 #define SPARSE_SMS_BOOL_OP_DECLS(M, S) \ |
|
|
171 SPARSE_BOOL_OP_DECL (mx_el_and, M, S); \ |
|
|
172 SPARSE_BOOL_OP_DECL (mx_el_or, M, S); |
|
|
173 |
|
|
174 #define SPARSE_SMS_BOOL_OP(F, OP, M, S, LHS_ZERO, RHS_ZERO) \ |
|
|
175 SparseBoolMatrix \ |
|
|
176 F (const M& m, const S& s) \ |
|
|
177 { \ |
|
|
178 /* Count num of non-zero elements */ \ |
|
|
179 int nel = 0; \ |
|
|
180 int nz = m.nnz (); \ |
|
|
181 if (LHS_ZERO OP (s != RHS_ZERO)) \ |
|
|
182 nel += m.numel() - nz; \ |
|
|
183 for (int i = 0; i < nz; i++) \ |
|
|
184 if ((m.data(i) != LHS_ZERO) OP (s != RHS_ZERO))\ |
|
|
185 nel++; \ |
|
|
186 \ |
|
|
187 int nr = m.rows (); \ |
|
|
188 int nc = m.cols (); \ |
|
|
189 SparseBoolMatrix r (nr, nc, nel); \ |
|
|
190 \ |
|
|
191 if (nr > 0 && nc > 0) \ |
|
|
192 { \ |
|
|
193 if (LHS_ZERO OP (s != RHS_ZERO)) \ |
|
|
194 { \ |
|
|
195 int ii = 0; \ |
|
|
196 r.cidx (0) = 0; \ |
|
|
197 for (int j = 0; j < nc; j++) \ |
|
|
198 { \ |
|
|
199 for (int i = 0; i < nr; i++) \ |
|
|
200 { \ |
|
|
201 bool el = (m.elem(i, j) != LHS_ZERO) OP (s != RHS_ZERO); \ |
|
|
202 if (el) \ |
|
|
203 { \ |
|
|
204 r.data(ii) = el; \ |
|
|
205 r.ridx(ii++) = i; \ |
|
|
206 } \ |
|
|
207 } \ |
|
|
208 r.cidx(j+1) = ii; \ |
|
|
209 } \ |
|
|
210 } \ |
|
|
211 else \ |
|
|
212 { \ |
|
|
213 int ii = 0; \ |
|
|
214 r.cidx (0) = 0; \ |
|
|
215 for (int j = 0; j < nc; j++) \ |
|
|
216 { \ |
|
|
217 for (int i = m.cidx(j); i < m.cidx(j+1); i++) \ |
|
|
218 { \ |
|
|
219 bool el = (m.data(i) != LHS_ZERO) OP (s != RHS_ZERO); \ |
|
|
220 if (el) \ |
|
|
221 { \ |
|
|
222 r.data(ii) = el; \ |
|
|
223 r.ridx(ii++) = m.ridx(i); \ |
|
|
224 } \ |
|
|
225 } \ |
|
|
226 r.cidx(j+1) = ii; \ |
|
|
227 } \ |
|
|
228 } \ |
|
|
229 } \ |
|
|
230 return r; \ |
|
|
231 } |
|
|
232 |
|
|
233 #define SPARSE_SMS_BOOL_OPS2(M, S, LHS_ZERO, RHS_ZERO) \ |
|
|
234 SPARSE_SMS_BOOL_OP (mx_el_and, &&, M, S, LHS_ZERO, RHS_ZERO) \ |
|
|
235 SPARSE_SMS_BOOL_OP (mx_el_or, ||, M, S, LHS_ZERO, RHS_ZERO) |
|
|
236 |
|
|
237 #define SPARSE_SMS_BOOL_OPS(M, S, ZERO) \ |
|
|
238 SPARSE_SMS_BOOL_OPS2(M, S, ZERO, ZERO) |
|
|
239 |
|
|
240 #define SPARSE_SMS_OP_DECLS(R1, R2, M, S) \ |
|
|
241 SPARSE_SMS_BIN_OP_DECLS (R1, R2, M, S) \ |
|
|
242 SPARSE_SMS_CMP_OP_DECLS (M, S) \ |
|
|
243 SPARSE_SMS_BOOL_OP_DECLS (M, S) |
|
|
244 |
|
|
245 // scalar by matrix operations. |
|
|
246 |
|
|
247 #define SPARSE_SSM_BIN_OP_DECLS(R1, R2, S, M) \ |
|
|
248 SPARSE_BIN_OP_DECL (R1, operator +, S, M); \ |
|
|
249 SPARSE_BIN_OP_DECL (R1, operator -, S, M); \ |
|
|
250 SPARSE_BIN_OP_DECL (R2, operator *, S, M); \ |
|
|
251 SPARSE_BIN_OP_DECL (R2, operator /, S, M); |
|
|
252 |
|
|
253 #define SPARSE_SSM_BIN_OP_1(R, F, OP, S, M) \ |
|
|
254 R \ |
|
|
255 F (const S& s, const M& m) \ |
|
|
256 { \ |
|
|
257 int nr = m.rows (); \ |
|
|
258 int nc = m.cols (); \ |
|
|
259 \ |
|
|
260 R r (nr, nc, (s OP 0.0)); \ |
|
|
261 \ |
|
|
262 for (int j = 0; j < nc; j++) \ |
|
|
263 for (int i = m.cidx (j); i < m.cidx (j+1); i++) \ |
|
|
264 r.elem (m.ridx (i), j) = s OP m.data (i); \ |
|
|
265 \ |
|
|
266 return r; \ |
|
|
267 } |
|
|
268 |
|
|
269 #define SPARSE_SSM_BIN_OP_2(R, F, OP, S, M) \ |
|
|
270 R \ |
|
|
271 F (const S& s, const M& m) \ |
|
|
272 { \ |
|
|
273 int nr = m.rows (); \ |
|
|
274 int nc = m.cols (); \ |
|
|
275 int nz = m.nnz (); \ |
|
|
276 \ |
|
|
277 R r (nr, nc, nz); \ |
|
|
278 \ |
|
|
279 for (int i = 0; i < nz; i++) \ |
|
|
280 { \ |
|
|
281 r.data(i) = s OP m.data(i); \ |
|
|
282 r.ridx(i) = m.ridx(i); \ |
|
|
283 } \ |
|
|
284 for (int i = 0; i < nc + 1; i++) \ |
|
|
285 r.cidx(i) = m.cidx(i); \ |
|
|
286 \ |
|
|
287 r.maybe_compress(true); \ |
|
|
288 return r; \ |
|
|
289 } |
|
|
290 |
|
|
291 #define SPARSE_SSM_BIN_OPS(R1, R2, S, M) \ |
|
|
292 SPARSE_SSM_BIN_OP_1 (R1, operator +, +, S, M) \ |
|
|
293 SPARSE_SSM_BIN_OP_1 (R1, operator -, -, S, M) \ |
|
|
294 SPARSE_SSM_BIN_OP_2 (R2, operator *, *, S, M) \ |
|
|
295 SPARSE_SSM_BIN_OP_2 (R2, operator /, /, S, M) |
|
|
296 |
|
|
297 #define SPARSE_SSM_CMP_OP_DECLS(S, M) \ |
|
|
298 SPARSE_CMP_OP_DECL (mx_el_lt, S, M); \ |
|
|
299 SPARSE_CMP_OP_DECL (mx_el_le, S, M); \ |
|
|
300 SPARSE_CMP_OP_DECL (mx_el_ge, S, M); \ |
|
|
301 SPARSE_CMP_OP_DECL (mx_el_gt, S, M); \ |
|
|
302 SPARSE_CMP_OP_DECL (mx_el_eq, S, M); \ |
|
|
303 SPARSE_CMP_OP_DECL (mx_el_ne, S, M); |
|
|
304 |
|
|
305 #define SPARSE_SSM_EQNE_OP_DECLS(S, M) \ |
|
|
306 SPARSE_CMP_OP_DECL (mx_el_eq, S, M); \ |
|
|
307 SPARSE_CMP_OP_DECL (mx_el_ne, S, M); |
|
|
308 |
|
|
309 #define SPARSE_SSM_CMP_OP(F, OP, S, SZ, SC, M, MZ, MC) \ |
|
|
310 SparseBoolMatrix \ |
|
|
311 F (const S& s, const M& m) \ |
|
|
312 { \ |
|
|
313 /* Count num of non-zero elements */ \ |
|
|
314 int nel = 0; \ |
|
|
315 int nz = m.nnz (); \ |
|
|
316 if (SC (s) OP MC (MZ)) \ |
|
|
317 nel += m.numel() - nz; \ |
|
|
318 for (int i = 0; i < nz; i++) \ |
|
|
319 if (SC (s) OP MC (m.data (i))) \ |
|
|
320 nel++; \ |
|
|
321 \ |
|
|
322 int nr = m.rows (); \ |
|
|
323 int nc = m.cols (); \ |
|
|
324 SparseBoolMatrix r (nr, nc, nel); \ |
|
|
325 \ |
|
|
326 if (nr > 0 && nc > 0) \ |
|
|
327 { \ |
|
|
328 if (SC (s) OP MC (MZ))\ |
|
|
329 { \ |
|
|
330 int ii = 0; \ |
|
|
331 r.cidx (0) = 0; \ |
|
|
332 for (int j = 0; j < nc; j++) \ |
|
|
333 { \ |
|
|
334 for (int i = 0; i < nr; i++) \ |
|
|
335 { \ |
|
|
336 bool el = SC (s) OP MC (m.elem(i, j)); \ |
|
|
337 if (el) \ |
|
|
338 { \ |
|
|
339 r.data(ii) = el; \ |
|
|
340 r.ridx(ii++) = i; \ |
|
|
341 } \ |
|
|
342 } \ |
|
|
343 r.cidx(j+1) = ii; \ |
|
|
344 } \ |
|
|
345 } \ |
|
|
346 else \ |
|
|
347 { \ |
|
|
348 int ii = 0; \ |
|
|
349 r.cidx (0) = 0; \ |
|
|
350 for (int j = 0; j < nc; j++) \ |
|
|
351 { \ |
|
|
352 for (int i = m.cidx(j); i < m.cidx(j+1); i++) \ |
|
|
353 { \ |
|
|
354 bool el = SC (s) OP MC (m.data(i)); \ |
|
|
355 if (el) \ |
|
|
356 { \ |
|
|
357 r.data(ii) = el; \ |
|
|
358 r.ridx(ii++) = m.ridx(i); \ |
|
|
359 } \ |
|
|
360 } \ |
|
|
361 r.cidx(j+1) = ii; \ |
|
|
362 } \ |
|
|
363 } \ |
|
|
364 } \ |
|
|
365 return r; \ |
|
|
366 } |
|
|
367 |
|
|
368 #define SPARSE_SSM_CMP_OPS(S, SZ, SC, M, MZ, MC) \ |
|
|
369 SPARSE_SSM_CMP_OP (mx_el_lt, <, S, SZ, SC, M, MZ, MC) \ |
|
|
370 SPARSE_SSM_CMP_OP (mx_el_le, <=, S, SZ, SC, M, MZ, MC) \ |
|
|
371 SPARSE_SSM_CMP_OP (mx_el_ge, >=, S, SZ, SC, M, MZ, MC) \ |
|
|
372 SPARSE_SSM_CMP_OP (mx_el_gt, >, S, SZ, SC, M, MZ, MC) \ |
|
|
373 SPARSE_SSM_CMP_OP (mx_el_eq, ==, S, SZ, , M, MZ, ) \ |
|
|
374 SPARSE_SSM_CMP_OP (mx_el_ne, !=, S, SZ, , M, MZ, ) |
|
|
375 |
|
|
376 #define SPARSE_SSM_EQNE_OPS(S, SZ, SC, M, MZ, MC) \ |
|
|
377 SPARSE_SSM_CMP_OP (mx_el_eq, ==, S, SZ, , M, MZ, ) \ |
|
|
378 SPARSE_SSM_CMP_OP (mx_el_ne, !=, S, SZ, , M, MZ, ) |
|
|
379 |
|
|
380 #define SPARSE_SSM_BOOL_OP_DECLS(S, M) \ |
|
|
381 SPARSE_BOOL_OP_DECL (mx_el_and, S, M); \ |
|
|
382 SPARSE_BOOL_OP_DECL (mx_el_or, S, M); \ |
|
|
383 |
|
|
384 #define SPARSE_SSM_BOOL_OP(F, OP, S, M, LHS_ZERO, RHS_ZERO) \ |
|
|
385 SparseBoolMatrix \ |
|
|
386 F (const S& s, const M& m) \ |
|
|
387 { \ |
|
|
388 /* Count num of non-zero elements */ \ |
|
|
389 int nel = 0; \ |
|
|
390 int nz = m.nnz (); \ |
|
|
391 if ((s != LHS_ZERO) OP RHS_ZERO) \ |
|
|
392 nel += m.numel() - nz; \ |
|
|
393 for (int i = 0; i < nz; i++) \ |
|
|
394 if ((s != LHS_ZERO) OP m.data(i) != RHS_ZERO) \ |
|
|
395 nel++; \ |
|
|
396 \ |
|
|
397 int nr = m.rows (); \ |
|
|
398 int nc = m.cols (); \ |
|
|
399 SparseBoolMatrix r (nr, nc, nel); \ |
|
|
400 \ |
|
|
401 if (nr > 0 && nc > 0) \ |
|
|
402 { \ |
|
|
403 if ((s != LHS_ZERO) OP RHS_ZERO) \ |
|
|
404 { \ |
|
|
405 int ii = 0; \ |
|
|
406 r.cidx (0) = 0; \ |
|
|
407 for (int j = 0; j < nc; j++) \ |
|
|
408 { \ |
|
|
409 for (int i = 0; i < nr; i++) \ |
|
|
410 { \ |
|
|
411 bool el = (s != LHS_ZERO) OP (m.elem(i, j) != RHS_ZERO); \ |
|
|
412 if (el) \ |
|
|
413 { \ |
|
|
414 r.data(ii) = el; \ |
|
|
415 r.ridx(ii++) = i; \ |
|
|
416 } \ |
|
|
417 } \ |
|
|
418 r.cidx(j+1) = ii; \ |
|
|
419 } \ |
|
|
420 } \ |
|
|
421 else \ |
|
|
422 { \ |
|
|
423 int ii = 0; \ |
|
|
424 r.cidx (0) = 0; \ |
|
|
425 for (int j = 0; j < nc; j++) \ |
|
|
426 { \ |
|
|
427 for (int i = m.cidx(j); i < m.cidx(j+1); i++) \ |
|
|
428 { \ |
|
|
429 bool el = (s != LHS_ZERO) OP (m.data(i) != RHS_ZERO); \ |
|
|
430 if (el) \ |
|
|
431 { \ |
|
|
432 r.data(ii) = el; \ |
|
|
433 r.ridx(ii++) = m.ridx(i); \ |
|
|
434 } \ |
|
|
435 } \ |
|
|
436 r.cidx(j+1) = ii; \ |
|
|
437 } \ |
|
|
438 } \ |
|
|
439 } \ |
|
|
440 return r; \ |
|
|
441 } |
|
|
442 |
|
|
443 #define SPARSE_SSM_BOOL_OPS2(S, M, LHS_ZERO, RHS_ZERO) \ |
|
|
444 SPARSE_SSM_BOOL_OP (mx_el_and, &&, S, M, LHS_ZERO, RHS_ZERO) \ |
|
|
445 SPARSE_SSM_BOOL_OP (mx_el_or, ||, S, M, LHS_ZERO, RHS_ZERO) |
|
|
446 |
|
|
447 #define SPARSE_SSM_BOOL_OPS(S, M, ZERO) \ |
|
|
448 SPARSE_SSM_BOOL_OPS2(S, M, ZERO, ZERO) |
|
|
449 |
|
|
450 #define SPARSE_SSM_OP_DECLS(R1, R2, S, M) \ |
|
|
451 SPARSE_SSM_BIN_OP_DECLS (R1, R2, S, M) \ |
|
|
452 SPARSE_SSM_CMP_OP_DECLS (S, M) \ |
|
|
453 SPARSE_SSM_BOOL_OP_DECLS (S, M) \ |
|
|
454 |
|
|
455 // matrix by matrix operations. |
|
|
456 |
|
|
457 #define SPARSE_SMSM_BIN_OP_DECLS(R1, R2, M1, M2) \ |
|
|
458 SPARSE_BIN_OP_DECL (R1, operator +, M1, M2); \ |
|
|
459 SPARSE_BIN_OP_DECL (R1, operator -, M1, M2); \ |
|
|
460 SPARSE_BIN_OP_DECL (R2, product, M1, M2); \ |
|
|
461 SPARSE_BIN_OP_DECL (R2, quotient, M1, M2); |
|
|
462 |
|
|
463 #define SPARSE_SMSM_BIN_OP_1(R, F, OP, M1, M2) \ |
|
|
464 R \ |
|
|
465 F (const M1& m1, const M2& m2) \ |
|
|
466 { \ |
|
|
467 R r; \ |
|
|
468 \ |
|
|
469 int m1_nr = m1.rows (); \ |
|
|
470 int m1_nc = m1.cols (); \ |
|
|
471 \ |
|
|
472 int m2_nr = m2.rows (); \ |
|
|
473 int m2_nc = m2.cols (); \ |
|
|
474 \ |
|
|
475 if (m1_nr != m2_nr || m1_nc != m2_nc) \ |
|
|
476 gripe_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \ |
|
|
477 else \ |
|
|
478 { \ |
|
|
479 r = R (m1_nr, m1_nc, (m1.nnz () + m2.nnz ())); \ |
|
|
480 \ |
|
|
481 int jx = 0; \ |
|
|
482 r.cidx (0) = 0; \ |
|
|
483 for (int i = 0 ; i < m1_nc ; i++) \ |
|
|
484 { \ |
|
|
485 int ja = m1.cidx(i); \ |
|
|
486 int ja_max = m1.cidx(i+1); \ |
|
|
487 bool ja_lt_max= ja < ja_max; \ |
|
|
488 \ |
|
|
489 int jb = m2.cidx(i); \ |
|
|
490 int jb_max = m2.cidx(i+1); \ |
|
|
491 bool jb_lt_max = jb < jb_max; \ |
|
|
492 \ |
|
|
493 while (ja_lt_max || jb_lt_max ) \ |
|
|
494 { \ |
|
|
495 OCTAVE_QUIT; \ |
|
|
496 if ((! jb_lt_max) || \ |
|
|
497 (ja_lt_max && (m1.ridx(ja) < m2.ridx(jb)))) \ |
|
|
498 { \ |
|
|
499 r.ridx(jx) = m1.ridx(ja); \ |
|
|
500 r.data(jx) = m1.data(ja) OP 0.; \ |
|
|
501 jx++; \ |
|
|
502 ja++; \ |
|
|
503 ja_lt_max= ja < ja_max; \ |
|
|
504 } \ |
|
|
505 else if (( !ja_lt_max ) || \ |
|
|
506 (jb_lt_max && (m2.ridx(jb) < m1.ridx(ja)) ) ) \ |
|
|
507 { \ |
|
|
508 r.ridx(jx) = m2.ridx(jb); \ |
|
|
509 r.data(jx) = 0. OP m2.data(jb); \ |
|
|
510 jx++; \ |
|
|
511 jb++; \ |
|
|
512 jb_lt_max= jb < jb_max; \ |
|
|
513 } \ |
|
|
514 else \ |
|
|
515 { \ |
|
|
516 if ((m1.data(ja) OP m2.data(jb)) != 0.) \ |
|
|
517 { \ |
|
|
518 r.data(jx) = m1.data(ja) OP m2.data(jb); \ |
|
|
519 r.ridx(jx) = m1.ridx(ja); \ |
|
|
520 jx++; \ |
|
|
521 } \ |
|
|
522 ja++; \ |
|
|
523 ja_lt_max= ja < ja_max; \ |
|
|
524 jb++; \ |
|
|
525 jb_lt_max= jb < jb_max; \ |
|
|
526 } \ |
|
|
527 } \ |
|
|
528 r.cidx(i+1) = jx; \ |
|
|
529 } \ |
|
|
530 \ |
|
|
531 r.maybe_compress (); \ |
|
|
532 } \ |
|
|
533 \ |
|
|
534 return r; \ |
|
|
535 } |
|
|
536 |
|
|
537 #define SPARSE_SMSM_BIN_OP_2(R, F, OP, M1, M2) \ |
|
|
538 R \ |
|
|
539 F (const M1& m1, const M2& m2) \ |
|
|
540 { \ |
|
|
541 R r; \ |
|
|
542 \ |
|
|
543 int m1_nr = m1.rows (); \ |
|
|
544 int m1_nc = m1.cols (); \ |
|
|
545 \ |
|
|
546 int m2_nr = m2.rows (); \ |
|
|
547 int m2_nc = m2.cols (); \ |
|
|
548 \ |
|
|
549 if (m1_nr != m2_nr || m1_nc != m2_nc) \ |
|
|
550 gripe_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \ |
|
|
551 else \ |
|
|
552 { \ |
|
|
553 r = R (m1_nr, m1_nc, (m1.nnz () > m2.nnz () ? m1.nnz () : m2.nnz ())); \ |
|
|
554 \ |
|
|
555 int jx = 0; \ |
|
|
556 r.cidx (0) = 0; \ |
|
|
557 for (int i = 0 ; i < m1_nc ; i++) \ |
|
|
558 { \ |
|
|
559 int ja = m1.cidx(i); \ |
|
|
560 int ja_max = m1.cidx(i+1); \ |
|
|
561 bool ja_lt_max= ja < ja_max; \ |
|
|
562 \ |
|
|
563 int jb = m2.cidx(i); \ |
|
|
564 int jb_max = m2.cidx(i+1); \ |
|
|
565 bool jb_lt_max = jb < jb_max; \ |
|
|
566 \ |
|
|
567 while (ja_lt_max || jb_lt_max ) \ |
|
|
568 { \ |
|
|
569 OCTAVE_QUIT; \ |
|
|
570 if ((! jb_lt_max) || \ |
|
|
571 (ja_lt_max && (m1.ridx(ja) < m2.ridx(jb)))) \ |
|
|
572 { \ |
|
|
573 ja++; ja_lt_max= ja < ja_max; \ |
|
|
574 } \ |
|
|
575 else if (( !ja_lt_max ) || \ |
|
|
576 (jb_lt_max && (m2.ridx(jb) < m1.ridx(ja)) ) ) \ |
|
|
577 { \ |
|
|
578 jb++; jb_lt_max= jb < jb_max; \ |
|
|
579 } \ |
|
|
580 else \ |
|
|
581 { \ |
|
|
582 if ((m1.data(ja) OP m2.data(jb)) != 0.) \ |
|
|
583 { \ |
|
|
584 r.data(jx) = m1.data(ja) OP m2.data(jb); \ |
|
|
585 r.ridx(jx) = m1.ridx(ja); \ |
|
|
586 jx++; \ |
|
|
587 } \ |
|
|
588 ja++; ja_lt_max= ja < ja_max; \ |
|
|
589 jb++; jb_lt_max= jb < jb_max; \ |
|
|
590 } \ |
|
|
591 } \ |
|
|
592 r.cidx(i+1) = jx; \ |
|
|
593 } \ |
|
|
594 \ |
|
|
595 r.maybe_compress (); \ |
|
|
596 } \ |
|
|
597 \ |
|
|
598 return r; \ |
|
|
599 } |
|
|
600 |
|
|
601 #define SPARSE_SMSM_BIN_OP_3(R, F, OP, M1, M2) \ |
|
|
602 R \ |
|
|
603 F (const M1& m1, const M2& m2) \ |
|
|
604 { \ |
|
|
605 R r; \ |
|
|
606 \ |
|
|
607 int m1_nr = m1.rows (); \ |
|
|
608 int m1_nc = m1.cols (); \ |
|
|
609 \ |
|
|
610 int m2_nr = m2.rows (); \ |
|
|
611 int m2_nc = m2.cols (); \ |
|
|
612 \ |
|
|
613 if (m1_nr != m2_nr || m1_nc != m2_nc) \ |
|
|
614 gripe_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \ |
|
|
615 else \ |
|
|
616 { \ |
|
|
617 \ |
|
|
618 /* XXX FIXME XXX Kludge... Always double/Complex, so Complex () */ \ |
|
|
619 r = R (m1_nr, m1_nc, (Complex () OP Complex ())); \ |
|
|
620 \ |
|
|
621 for (int i = 0 ; i < m1_nc ; i++) \ |
|
|
622 { \ |
|
|
623 int ja = m1.cidx(i); \ |
|
|
624 int ja_max = m1.cidx(i+1); \ |
|
|
625 bool ja_lt_max= ja < ja_max; \ |
|
|
626 \ |
|
|
627 int jb = m2.cidx(i); \ |
|
|
628 int jb_max = m2.cidx(i+1); \ |
|
|
629 bool jb_lt_max = jb < jb_max; \ |
|
|
630 \ |
|
|
631 while (ja_lt_max || jb_lt_max ) \ |
|
|
632 { \ |
|
|
633 OCTAVE_QUIT; \ |
|
|
634 if ((! jb_lt_max) || \ |
|
|
635 (ja_lt_max && (m1.ridx(ja) < m2.ridx(jb)))) \ |
|
|
636 { \ |
|
|
637 /* keep those kludges coming */ \ |
|
|
638 r.elem(m1.ridx(ja),i) = m1.data(ja) OP Complex (); \ |
|
|
639 ja++; \ |
|
|
640 ja_lt_max= ja < ja_max; \ |
|
|
641 } \ |
|
|
642 else if (( !ja_lt_max ) || \ |
|
|
643 (jb_lt_max && (m2.ridx(jb) < m1.ridx(ja)) ) ) \ |
|
|
644 { \ |
|
|
645 /* keep those kludges coming */ \ |
|
|
646 r.elem(m2.ridx(jb),i) = Complex () OP m2.data(jb); \ |
|
|
647 jb++; \ |
|
|
648 jb_lt_max= jb < jb_max; \ |
|
|
649 } \ |
|
|
650 else \ |
|
|
651 { \ |
|
|
652 r.elem(m1.ridx(ja),i) = m1.data(ja) OP m2.data(jb); \ |
|
|
653 ja++; \ |
|
|
654 ja_lt_max= ja < ja_max; \ |
|
|
655 jb++; \ |
|
|
656 jb_lt_max= jb < jb_max; \ |
|
|
657 } \ |
|
|
658 } \ |
|
|
659 } \ |
|
|
660 r.maybe_compress (true); \ |
|
|
661 } \ |
|
|
662 \ |
|
|
663 return r; \ |
|
|
664 } |
|
|
665 |
|
|
666 // Note that SM ./ SM needs to take into account the NaN and Inf values |
|
|
667 // implied by the division by zero. |
|
|
668 // XXX FIXME XXX Are the NaNs double(NaN) or Complex(NaN,Nan) in the complex |
|
|
669 // case? |
|
|
670 #define SPARSE_SMSM_BIN_OPS(R1, R2, M1, M2) \ |
|
|
671 SPARSE_SMSM_BIN_OP_1 (R1, operator +, +, M1, M2) \ |
|
|
672 SPARSE_SMSM_BIN_OP_1 (R1, operator -, -, M1, M2) \ |
|
|
673 SPARSE_SMSM_BIN_OP_2 (R2, product, *, M1, M2) \ |
|
|
674 SPARSE_SMSM_BIN_OP_3 (R2, quotient, /, M1, M2) |
|
|
675 |
|
|
676 #define SPARSE_SMSM_CMP_OP_DECLS(M1, M2) \ |
|
|
677 SPARSE_CMP_OP_DECL (mx_el_lt, M1, M2); \ |
|
|
678 SPARSE_CMP_OP_DECL (mx_el_le, M1, M2); \ |
|
|
679 SPARSE_CMP_OP_DECL (mx_el_ge, M1, M2); \ |
|
|
680 SPARSE_CMP_OP_DECL (mx_el_gt, M1, M2); \ |
|
|
681 SPARSE_CMP_OP_DECL (mx_el_eq, M1, M2); \ |
|
|
682 SPARSE_CMP_OP_DECL (mx_el_ne, M1, M2); |
|
|
683 |
|
|
684 #define SPARSE_SMSM_EQNE_OP_DECLS(M1, M2) \ |
|
|
685 SPARSE_CMP_OP_DECL (mx_el_eq, M1, M2); \ |
|
|
686 SPARSE_CMP_OP_DECL (mx_el_ne, M1, M2); |
|
|
687 |
|
|
688 #define SPARSE_SMSM_CMP_OP(F, OP, M1, C1, M2, C2) \ |
|
|
689 SparseBoolMatrix \ |
|
|
690 F (const M1& m1, const M2& m2) \ |
|
|
691 { \ |
|
|
692 SparseBoolMatrix r; \ |
|
|
693 \ |
|
|
694 int m1_nr = m1.rows (); \ |
|
|
695 int m1_nc = m1.cols (); \ |
|
|
696 \ |
|
|
697 int m2_nr = m2.rows (); \ |
|
|
698 int m2_nc = m2.cols (); \ |
|
|
699 \ |
|
|
700 if (m1_nr == m2_nr && m1_nc == m2_nc) \ |
|
|
701 { \ |
|
|
702 if (m1_nr != 0 || m1_nc != 0) \ |
|
|
703 { \ |
|
|
704 /* Count num of non-zero elements */ \ |
|
|
705 int nel = 0; \ |
|
|
706 for (int j = 0; j < m1_nc; j++) \ |
|
|
707 for (int i = 0; i < m1_nr; i++) \ |
|
|
708 if (C1 (m1.elem(i, j)) OP C2 (m2.elem(i, j))) \ |
|
|
709 nel++; \ |
|
|
710 \ |
|
|
711 r = SparseBoolMatrix (m1_nr, m1_nc, nel); \ |
|
|
712 \ |
|
|
713 int ii = 0; \ |
|
|
714 r.cidx (0) = 0; \ |
|
|
715 for (int j = 0; j < m1_nc; j++) \ |
|
|
716 { \ |
|
|
717 for (int i = 0; i < m1_nr; i++) \ |
|
|
718 { \ |
|
|
719 bool el = C1 (m1.elem(i, j)) OP C2 (m2.elem(i, j)); \ |
|
|
720 if (el) \ |
|
|
721 { \ |
|
|
722 r.data(ii) = el; \ |
|
|
723 r.ridx(ii++) = i; \ |
|
|
724 } \ |
|
|
725 } \ |
|
|
726 r.cidx(j+1) = ii; \ |
|
|
727 } \ |
|
|
728 } \ |
|
|
729 } \ |
|
|
730 else \ |
|
|
731 { \ |
|
|
732 if ((m1_nr != 0 || m1_nc != 0) && (m2_nr != 0 || m2_nc != 0)) \ |
|
|
733 gripe_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \ |
|
|
734 } \ |
|
|
735 return r; \ |
|
|
736 } |
|
|
737 |
|
|
738 #define SPARSE_SMSM_CMP_OPS(M1, Z1, C1, M2, Z2, C2) \ |
|
|
739 SPARSE_SMSM_CMP_OP (mx_el_lt, <, M1, C1, M2, C2) \ |
|
|
740 SPARSE_SMSM_CMP_OP (mx_el_le, <=, M1, C1, M2, C2) \ |
|
|
741 SPARSE_SMSM_CMP_OP (mx_el_ge, >=, M1, C1, M2, C2) \ |
|
|
742 SPARSE_SMSM_CMP_OP (mx_el_gt, >, M1, C1, M2, C2) \ |
|
|
743 SPARSE_SMSM_CMP_OP (mx_el_eq, ==, M1, , M2, ) \ |
|
|
744 SPARSE_SMSM_CMP_OP (mx_el_ne, !=, M1, , M2, ) |
|
|
745 |
|
|
746 #define SPARSE_SMSM_EQNE_OPS(M1, Z1, C1, M2, Z2, C2) \ |
|
|
747 SPARSE_SMSM_CMP_OP (mx_el_eq, ==, M1, , M2, ) \ |
|
|
748 SPARSE_SMSM_CMP_OP (mx_el_ne, !=, M1, , M2, ) |
|
|
749 |
|
|
750 #define SPARSE_SMSM_BOOL_OP_DECLS(M1, M2) \ |
|
|
751 SPARSE_BOOL_OP_DECL (mx_el_and, M1, M2); \ |
|
|
752 SPARSE_BOOL_OP_DECL (mx_el_or, M1, M2); |
|
|
753 |
|
|
754 #define SPARSE_SMSM_BOOL_OP(F, OP, M1, M2, LHS_ZERO, RHS_ZERO) \ |
|
|
755 SparseBoolMatrix \ |
|
|
756 F (const M1& m1, const M2& m2) \ |
|
|
757 { \ |
|
|
758 SparseBoolMatrix r; \ |
|
|
759 \ |
|
|
760 int m1_nr = m1.rows (); \ |
|
|
761 int m1_nc = m1.cols (); \ |
|
|
762 \ |
|
|
763 int m2_nr = m2.rows (); \ |
|
|
764 int m2_nc = m2.cols (); \ |
|
|
765 \ |
|
|
766 if (m1_nr == m2_nr && m1_nc == m2_nc) \ |
|
|
767 { \ |
|
|
768 if (m1_nr != 0 || m1_nc != 0) \ |
|
|
769 { \ |
|
|
770 /* Count num of non-zero elements */ \ |
|
|
771 int nel = 0; \ |
|
|
772 for (int j = 0; j < m1_nc; j++) \ |
|
|
773 for (int i = 0; i < m1_nr; i++) \ |
|
|
774 if ((m1.elem(i, j) != LHS_ZERO) \ |
|
|
775 OP (m2.elem(i, j) != RHS_ZERO)) \ |
|
|
776 nel++; \ |
|
|
777 \ |
|
|
778 r = SparseBoolMatrix (m1_nr, m1_nc, nel); \ |
|
|
779 \ |
|
|
780 int ii = 0; \ |
|
|
781 r.cidx (0) = 0; \ |
|
|
782 for (int j = 0; j < m1_nc; j++) \ |
|
|
783 { \ |
|
|
784 for (int i = 0; i < m1_nr; i++) \ |
|
|
785 { \ |
|
|
786 bool el = (m1.elem(i, j) != LHS_ZERO) \ |
|
|
787 OP (m2.elem(i, j) != RHS_ZERO); \ |
|
|
788 if (el) \ |
|
|
789 { \ |
|
|
790 r.data(ii) = el; \ |
|
|
791 r.ridx(ii++) = i; \ |
|
|
792 } \ |
|
|
793 } \ |
|
|
794 r.cidx(j+1) = ii; \ |
|
|
795 } \ |
|
|
796 } \ |
|
|
797 } \ |
|
|
798 else \ |
|
|
799 { \ |
|
|
800 if ((m1_nr != 0 || m1_nc != 0) && (m2_nr != 0 || m2_nc != 0)) \ |
|
|
801 gripe_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \ |
|
|
802 } \ |
|
|
803 return r; \ |
|
|
804 } |
|
|
805 |
|
|
806 #define SPARSE_SMSM_BOOL_OPS2(M1, M2, LHS_ZERO, RHS_ZERO) \ |
|
|
807 SPARSE_SMSM_BOOL_OP (mx_el_and, &&, M1, M2, LHS_ZERO, RHS_ZERO) \ |
|
|
808 SPARSE_SMSM_BOOL_OP (mx_el_or, ||, M1, M2, LHS_ZERO, RHS_ZERO) \ |
|
|
809 |
|
|
810 #define SPARSE_SMSM_BOOL_OPS(M1, M2, ZERO) \ |
|
|
811 SPARSE_SMSM_BOOL_OPS2(M1, M2, ZERO, ZERO) |
|
|
812 |
|
|
813 #define SPARSE_SMSM_OP_DECLS(R1, R2, M1, M2) \ |
|
|
814 SPARSE_SMSM_BIN_OP_DECLS (R1, R2, M1, M2) \ |
|
|
815 SPARSE_SMSM_CMP_OP_DECLS (M1, M2) \ |
|
|
816 SPARSE_SMSM_BOOL_OP_DECLS (M1, M2) |
|
|
817 |
|
|
818 // matrix by matrix operations. |
|
|
819 |
|
|
820 #define SPARSE_MSM_BIN_OP_DECLS(R1, R2, M1, M2) \ |
|
|
821 SPARSE_BIN_OP_DECL (R1, operator +, M1, M2); \ |
|
|
822 SPARSE_BIN_OP_DECL (R1, operator -, M1, M2); \ |
|
|
823 SPARSE_BIN_OP_DECL (R2, product, M1, M2); \ |
|
|
824 SPARSE_BIN_OP_DECL (R2, quotient, M1, M2); |
|
|
825 |
|
|
826 #define SPARSE_MSM_BIN_OP_1(R, F, OP, M1, M2) \ |
|
|
827 R \ |
|
|
828 F (const M1& m1, const M2& m2) \ |
|
|
829 { \ |
|
|
830 R r; \ |
|
|
831 \ |
|
|
832 int m1_nr = m1.rows (); \ |
|
|
833 int m1_nc = m1.cols (); \ |
|
|
834 \ |
|
|
835 int m2_nr = m2.rows (); \ |
|
|
836 int m2_nc = m2.cols (); \ |
|
|
837 \ |
|
|
838 if (m1_nr != m2_nr || m1_nc != m2_nc) \ |
|
|
839 gripe_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \ |
|
|
840 else \ |
|
|
841 { \ |
|
|
842 r = R (m1_nr, m1_nc); \ |
|
|
843 \ |
|
|
844 for (int j = 0; j < m1_nc; j++) \ |
|
|
845 for (int i = 0; i < m1_nr; i++) \ |
|
|
846 r.elem (i, j) = m1.elem (i, j) OP m2.elem (i, j); \ |
|
|
847 } \ |
|
|
848 return r; \ |
|
|
849 } |
|
|
850 |
|
|
851 #define SPARSE_MSM_BIN_OP_2(R, F, OP, M1, M2, ZERO) \ |
|
|
852 R \ |
|
|
853 F (const M1& m1, const M2& m2) \ |
|
|
854 { \ |
|
|
855 R r; \ |
|
|
856 \ |
|
|
857 int m1_nr = m1.rows (); \ |
|
|
858 int m1_nc = m1.cols (); \ |
|
|
859 \ |
|
|
860 int m2_nr = m2.rows (); \ |
|
|
861 int m2_nc = m2.cols (); \ |
|
|
862 \ |
|
|
863 if (m1_nr != m2_nr || m1_nc != m2_nc) \ |
|
|
864 gripe_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \ |
|
|
865 else \ |
|
|
866 { \ |
|
|
867 /* Count num of non-zero elements */ \ |
|
|
868 int nel = 0; \ |
|
|
869 for (int j = 0; j < m1_nc; j++) \ |
|
|
870 for (int i = 0; i < m1_nr; i++) \ |
|
|
871 if ((m1.elem(i, j) OP m2.elem(i, j)) != ZERO) \ |
|
|
872 nel++; \ |
|
|
873 \ |
|
|
874 r = R (m1_nr, m1_nc, nel); \ |
|
|
875 \ |
|
|
876 int ii = 0; \ |
|
|
877 r.cidx (0) = 0; \ |
|
|
878 for (int j = 0 ; j < m1_nc ; j++) \ |
|
|
879 { \ |
|
|
880 for (int i = 0 ; i < m1_nr ; i++) \ |
|
|
881 { \ |
|
|
882 if ((m1.elem(i, j) OP m2.elem(i, j)) != ZERO) \ |
|
|
883 { \ |
|
|
884 r.data (ii) = m1.elem(i, j) OP m2.elem(i,j); \ |
|
|
885 r.ridx (ii++) = i; \ |
|
|
886 } \ |
|
|
887 } \ |
|
|
888 r.cidx(j+1) = ii; \ |
|
|
889 } \ |
|
|
890 } \ |
|
|
891 \ |
|
|
892 return r; \ |
|
|
893 } |
|
|
894 |
|
|
895 // XXX FIXME XXX Pass a specific ZERO value |
|
|
896 #define SPARSE_MSM_BIN_OPS(R1, R2, M1, M2) \ |
|
|
897 SPARSE_MSM_BIN_OP_1 (R1, operator +, +, M1, M2) \ |
|
|
898 SPARSE_MSM_BIN_OP_1 (R1, operator -, -, M1, M2) \ |
|
|
899 SPARSE_MSM_BIN_OP_2 (R2, product, *, M1, M2, 0.0) \ |
|
|
900 SPARSE_MSM_BIN_OP_2 (R2, quotient, /, M1, M2, 0.0) |
|
|
901 |
|
|
902 #define SPARSE_MSM_CMP_OP_DECLS(M1, M2) \ |
|
|
903 SPARSE_CMP_OP_DECL (mx_el_lt, M1, M2); \ |
|
|
904 SPARSE_CMP_OP_DECL (mx_el_le, M1, M2); \ |
|
|
905 SPARSE_CMP_OP_DECL (mx_el_ge, M1, M2); \ |
|
|
906 SPARSE_CMP_OP_DECL (mx_el_gt, M1, M2); \ |
|
|
907 SPARSE_CMP_OP_DECL (mx_el_eq, M1, M2); \ |
|
|
908 SPARSE_CMP_OP_DECL (mx_el_ne, M1, M2); |
|
|
909 |
|
|
910 #define SPARSE_MSM_EQNE_OP_DECLS(M1, M2) \ |
|
|
911 SPARSE_CMP_OP_DECL (mx_el_eq, M1, M2); \ |
|
|
912 SPARSE_CMP_OP_DECL (mx_el_ne, M1, M2); |
|
|
913 |
|
|
914 #define SPARSE_MSM_CMP_OP(F, OP, M1, C1, M2, C2) \ |
|
|
915 SparseBoolMatrix \ |
|
|
916 F (const M1& m1, const M2& m2) \ |
|
|
917 { \ |
|
|
918 SparseBoolMatrix r; \ |
|
|
919 \ |
|
|
920 int m1_nr = m1.rows (); \ |
|
|
921 int m1_nc = m1.cols (); \ |
|
|
922 \ |
|
|
923 int m2_nr = m2.rows (); \ |
|
|
924 int m2_nc = m2.cols (); \ |
|
|
925 \ |
|
|
926 if (m1_nr == m2_nr && m1_nc == m2_nc) \ |
|
|
927 { \ |
|
|
928 if (m1_nr != 0 || m1_nc != 0) \ |
|
|
929 { \ |
|
|
930 /* Count num of non-zero elements */ \ |
|
|
931 int nel = 0; \ |
|
|
932 for (int j = 0; j < m1_nc; j++) \ |
|
|
933 for (int i = 0; i < m1_nr; i++) \ |
|
|
934 if (C1 (m1.elem(i, j)) OP C2 (m2.elem(i, j))) \ |
|
|
935 nel++; \ |
|
|
936 \ |
|
|
937 r = SparseBoolMatrix (m1_nr, m1_nc, nel); \ |
|
|
938 \ |
|
|
939 int ii = 0; \ |
|
|
940 r.cidx (0) = 0; \ |
|
|
941 for (int j = 0; j < m1_nc; j++) \ |
|
|
942 { \ |
|
|
943 for (int i = 0; i < m1_nr; i++) \ |
|
|
944 { \ |
|
|
945 bool el = C1 (m1.elem(i, j)) OP C2 (m2.elem(i, j)); \ |
|
|
946 if (el) \ |
|
|
947 { \ |
|
|
948 r.data(ii) = el; \ |
|
|
949 r.ridx(ii++) = i; \ |
|
|
950 } \ |
|
|
951 } \ |
|
|
952 r.cidx(j+1) = ii; \ |
|
|
953 } \ |
|
|
954 } \ |
|
|
955 } \ |
|
|
956 else \ |
|
|
957 { \ |
|
|
958 if ((m1_nr != 0 || m1_nc != 0) && (m2_nr != 0 || m2_nc != 0)) \ |
|
|
959 gripe_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \ |
|
|
960 } \ |
|
|
961 return r; \ |
|
|
962 } |
|
|
963 |
|
|
964 #define SPARSE_MSM_CMP_OPS(M1, Z1, C1, M2, Z2, C2) \ |
|
|
965 SPARSE_MSM_CMP_OP (mx_el_lt, <, M1, C1, M2, C2) \ |
|
|
966 SPARSE_MSM_CMP_OP (mx_el_le, <=, M1, C1, M2, C2) \ |
|
|
967 SPARSE_MSM_CMP_OP (mx_el_ge, >=, M1, C1, M2, C2) \ |
|
|
968 SPARSE_MSM_CMP_OP (mx_el_gt, >, M1, C1, M2, C2) \ |
|
|
969 SPARSE_MSM_CMP_OP (mx_el_eq, ==, M1, , M2, ) \ |
|
|
970 SPARSE_MSM_CMP_OP (mx_el_ne, !=, M1, , M2, ) |
|
|
971 |
|
|
972 #define SPARSE_MSM_EQNE_OPS(M1, Z1, C1, M2, Z2, C2) \ |
|
|
973 SPARSE_MSM_CMP_OP (mx_el_eq, ==, M1, , M2, ) \ |
|
|
974 SPARSE_MSM_CMP_OP (mx_el_ne, !=, M1, , M2, ) |
|
|
975 |
|
|
976 #define SPARSE_MSM_BOOL_OP_DECLS(M1, M2) \ |
|
|
977 SPARSE_BOOL_OP_DECL (mx_el_and, M1, M2); \ |
|
|
978 SPARSE_BOOL_OP_DECL (mx_el_or, M1, M2); |
|
|
979 |
|
|
980 #define SPARSE_MSM_BOOL_OP(F, OP, M1, M2, LHS_ZERO, RHS_ZERO) \ |
|
|
981 SparseBoolMatrix \ |
|
|
982 F (const M1& m1, const M2& m2) \ |
|
|
983 { \ |
|
|
984 SparseBoolMatrix r; \ |
|
|
985 \ |
|
|
986 int m1_nr = m1.rows (); \ |
|
|
987 int m1_nc = m1.cols (); \ |
|
|
988 \ |
|
|
989 int m2_nr = m2.rows (); \ |
|
|
990 int m2_nc = m2.cols (); \ |
|
|
991 \ |
|
|
992 if (m1_nr == m2_nr && m1_nc == m2_nc) \ |
|
|
993 { \ |
|
|
994 if (m1_nr != 0 || m1_nc != 0) \ |
|
|
995 { \ |
|
|
996 /* Count num of non-zero elements */ \ |
|
|
997 int nel = 0; \ |
|
|
998 for (int j = 0; j < m1_nc; j++) \ |
|
|
999 for (int i = 0; i < m1_nr; i++) \ |
|
|
1000 if ((m1.elem(i, j) != LHS_ZERO) \ |
|
|
1001 OP (m2.elem(i, j) != RHS_ZERO)) \ |
|
|
1002 nel++; \ |
|
|
1003 \ |
|
|
1004 r = SparseBoolMatrix (m1_nr, m1_nc, nel); \ |
|
|
1005 \ |
|
|
1006 int ii = 0; \ |
|
|
1007 r.cidx (0) = 0; \ |
|
|
1008 for (int j = 0; j < m1_nc; j++) \ |
|
|
1009 { \ |
|
|
1010 for (int i = 0; i < m1_nr; i++) \ |
|
|
1011 { \ |
|
|
1012 bool el = (m1.elem(i, j) != LHS_ZERO) \ |
|
|
1013 OP (m2.elem(i, j) != RHS_ZERO); \ |
|
|
1014 if (el) \ |
|
|
1015 { \ |
|
|
1016 r.data(ii) = el; \ |
|
|
1017 r.ridx(ii++) = i; \ |
|
|
1018 } \ |
|
|
1019 } \ |
|
|
1020 r.cidx(j+1) = ii; \ |
|
|
1021 } \ |
|
|
1022 } \ |
|
|
1023 } \ |
|
|
1024 else \ |
|
|
1025 { \ |
|
|
1026 if ((m1_nr != 0 || m1_nc != 0) && (m2_nr != 0 || m2_nc != 0)) \ |
|
|
1027 gripe_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \ |
|
|
1028 } \ |
|
|
1029 return r; \ |
|
|
1030 } |
|
|
1031 |
|
|
1032 #define SPARSE_MSM_BOOL_OPS2(M1, M2, LHS_ZERO, RHS_ZERO) \ |
|
|
1033 SPARSE_MSM_BOOL_OP (mx_el_and, &&, M1, M2, LHS_ZERO, RHS_ZERO) \ |
|
|
1034 SPARSE_MSM_BOOL_OP (mx_el_or, ||, M1, M2, LHS_ZERO, RHS_ZERO) \ |
|
|
1035 |
|
|
1036 #define SPARSE_MSM_BOOL_OPS(M1, M2, ZERO) \ |
|
|
1037 SPARSE_MSM_BOOL_OPS2(M1, M2, ZERO, ZERO) |
|
|
1038 |
|
|
1039 #define SPARSE_MSM_OP_DECLS(R1, R2, M1, M2) \ |
|
|
1040 SPARSE_MSM_BIN_OP_DECLS (R1, R2, M1, M2) \ |
|
|
1041 SPARSE_MSM_CMP_OP_DECLS (M1, M2) \ |
|
|
1042 SPARSE_MSM_BOOL_OP_DECLS (M1, M2) |
|
|
1043 |
|
|
1044 // matrix by matrix operations. |
|
|
1045 |
|
|
1046 #define SPARSE_SMM_BIN_OP_DECLS(R1, R2, M1, M2) \ |
|
|
1047 SPARSE_BIN_OP_DECL (R1, operator +, M1, M2); \ |
|
|
1048 SPARSE_BIN_OP_DECL (R1, operator -, M1, M2); \ |
|
|
1049 SPARSE_BIN_OP_DECL (R2, product, M1, M2); \ |
|
|
1050 SPARSE_BIN_OP_DECL (R2, quotient, M1, M2); |
|
|
1051 |
|
|
1052 #define SPARSE_SMM_BIN_OP_1(R, F, OP, M1, M2) \ |
|
|
1053 R \ |
|
|
1054 F (const M1& m1, const M2& m2) \ |
|
|
1055 { \ |
|
|
1056 R r; \ |
|
|
1057 \ |
|
|
1058 int m1_nr = m1.rows (); \ |
|
|
1059 int m1_nc = m1.cols (); \ |
|
|
1060 \ |
|
|
1061 int m2_nr = m2.rows (); \ |
|
|
1062 int m2_nc = m2.cols (); \ |
|
|
1063 \ |
|
|
1064 if (m1_nr != m2_nr || m1_nc != m2_nc) \ |
|
|
1065 gripe_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \ |
|
|
1066 else \ |
|
|
1067 { \ |
|
|
1068 r = R (m1_nr, m1_nc); \ |
|
|
1069 \ |
|
|
1070 for (int j = 0; j < m1_nc; j++) \ |
|
|
1071 for (int i = 0; i < m1_nr; i++) \ |
|
|
1072 r.elem (i, j) = m1.elem (i, j) OP m2.elem (i, j); \ |
|
|
1073 } \ |
|
|
1074 return r; \ |
|
|
1075 } |
|
|
1076 |
|
|
1077 #define SPARSE_SMM_BIN_OP_2(R, F, OP, M1, M2, ZERO) \ |
|
|
1078 R \ |
|
|
1079 F (const M1& m1, const M2& m2) \ |
|
|
1080 { \ |
|
|
1081 R r; \ |
|
|
1082 \ |
|
|
1083 int m1_nr = m1.rows (); \ |
|
|
1084 int m1_nc = m1.cols (); \ |
|
|
1085 \ |
|
|
1086 int m2_nr = m2.rows (); \ |
|
|
1087 int m2_nc = m2.cols (); \ |
|
|
1088 \ |
|
|
1089 if (m1_nr != m2_nr || m1_nc != m2_nc) \ |
|
|
1090 gripe_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \ |
|
|
1091 else \ |
|
|
1092 { \ |
|
|
1093 /* Count num of non-zero elements */ \ |
|
|
1094 int nel = 0; \ |
|
|
1095 for (int j = 0; j < m1_nc; j++) \ |
|
|
1096 for (int i = 0; i < m1_nr; i++) \ |
|
|
1097 if ((m1.elem(i, j) OP m2.elem(i, j)) != ZERO) \ |
|
|
1098 nel++; \ |
|
|
1099 \ |
|
|
1100 r = R (m1_nr, m1_nc, nel); \ |
|
|
1101 \ |
|
|
1102 int ii = 0; \ |
|
|
1103 r.cidx (0) = 0; \ |
|
|
1104 for (int j = 0 ; j < m1_nc ; j++) \ |
|
|
1105 { \ |
|
|
1106 for (int i = 0 ; i < m1_nr ; i++) \ |
|
|
1107 { \ |
|
|
1108 if ((m1.elem(i, j) OP m2.elem(i, j)) != ZERO) \ |
|
|
1109 { \ |
|
|
1110 r.data (ii) = m1.elem(i, j) OP m2.elem(i,j); \ |
|
|
1111 r.ridx (ii++) = i; \ |
|
|
1112 } \ |
|
|
1113 } \ |
|
|
1114 r.cidx(j+1) = ii; \ |
|
|
1115 } \ |
|
|
1116 } \ |
|
|
1117 \ |
|
|
1118 return r; \ |
|
|
1119 } |
|
|
1120 |
|
|
1121 // XXX FIXME XXX Pass a specific ZERO value |
|
|
1122 #define SPARSE_SMM_BIN_OPS(R1, R2, M1, M2) \ |
|
|
1123 SPARSE_SMM_BIN_OP_1 (R1, operator +, +, M1, M2) \ |
|
|
1124 SPARSE_SMM_BIN_OP_1 (R1, operator -, -, M1, M2) \ |
|
|
1125 SPARSE_SMM_BIN_OP_2 (R2, product, *, M1, M2, 0.0) \ |
|
|
1126 SPARSE_SMM_BIN_OP_2 (R2, quotient, /, M1, M2, 0.0) |
|
|
1127 |
|
|
1128 #define SPARSE_SMM_CMP_OP_DECLS(M1, M2) \ |
|
|
1129 SPARSE_CMP_OP_DECL (mx_el_lt, M1, M2); \ |
|
|
1130 SPARSE_CMP_OP_DECL (mx_el_le, M1, M2); \ |
|
|
1131 SPARSE_CMP_OP_DECL (mx_el_ge, M1, M2); \ |
|
|
1132 SPARSE_CMP_OP_DECL (mx_el_gt, M1, M2); \ |
|
|
1133 SPARSE_CMP_OP_DECL (mx_el_eq, M1, M2); \ |
|
|
1134 SPARSE_CMP_OP_DECL (mx_el_ne, M1, M2); |
|
|
1135 |
|
|
1136 #define SPARSE_SMM_EQNE_OP_DECLS(M1, M2) \ |
|
|
1137 SPARSE_CMP_OP_DECL (mx_el_eq, M1, M2); \ |
|
|
1138 SPARSE_CMP_OP_DECL (mx_el_ne, M1, M2); |
|
|
1139 |
|
|
1140 #define SPARSE_SMM_CMP_OP(F, OP, M1, C1, M2, C2) \ |
|
|
1141 SparseBoolMatrix \ |
|
|
1142 F (const M1& m1, const M2& m2) \ |
|
|
1143 { \ |
|
|
1144 SparseBoolMatrix r; \ |
|
|
1145 \ |
|
|
1146 int m1_nr = m1.rows (); \ |
|
|
1147 int m1_nc = m1.cols (); \ |
|
|
1148 \ |
|
|
1149 int m2_nr = m2.rows (); \ |
|
|
1150 int m2_nc = m2.cols (); \ |
|
|
1151 \ |
|
|
1152 if (m1_nr == m2_nr && m1_nc == m2_nc) \ |
|
|
1153 { \ |
|
|
1154 if (m1_nr != 0 || m1_nc != 0) \ |
|
|
1155 { \ |
|
|
1156 /* Count num of non-zero elements */ \ |
|
|
1157 int nel = 0; \ |
|
|
1158 for (int j = 0; j < m1_nc; j++) \ |
|
|
1159 for (int i = 0; i < m1_nr; i++) \ |
|
|
1160 if (C1 (m1.elem(i, j)) OP C2 (m2.elem(i, j))) \ |
|
|
1161 nel++; \ |
|
|
1162 \ |
|
|
1163 r = SparseBoolMatrix (m1_nr, m1_nc, nel); \ |
|
|
1164 \ |
|
|
1165 int ii = 0; \ |
|
|
1166 r.cidx (0) = 0; \ |
|
|
1167 for (int j = 0; j < m1_nc; j++) \ |
|
|
1168 { \ |
|
|
1169 for (int i = 0; i < m1_nr; i++) \ |
|
|
1170 { \ |
|
|
1171 bool el = C1 (m1.elem(i, j)) OP C2 (m2.elem(i, j)); \ |
|
|
1172 if (el) \ |
|
|
1173 { \ |
|
|
1174 r.data(ii) = el; \ |
|
|
1175 r.ridx(ii++) = i; \ |
|
|
1176 } \ |
|
|
1177 } \ |
|
|
1178 r.cidx(j+1) = ii; \ |
|
|
1179 } \ |
|
|
1180 } \ |
|
|
1181 } \ |
|
|
1182 else \ |
|
|
1183 { \ |
|
|
1184 if ((m1_nr != 0 || m1_nc != 0) && (m2_nr != 0 || m2_nc != 0)) \ |
|
|
1185 gripe_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \ |
|
|
1186 } \ |
|
|
1187 return r; \ |
|
|
1188 } |
|
|
1189 |
|
|
1190 #define SPARSE_SMM_CMP_OPS(M1, Z1, C1, M2, Z2, C2) \ |
|
|
1191 SPARSE_SMM_CMP_OP (mx_el_lt, <, M1, C1, M2, C2) \ |
|
|
1192 SPARSE_SMM_CMP_OP (mx_el_le, <=, M1, C1, M2, C2) \ |
|
|
1193 SPARSE_SMM_CMP_OP (mx_el_ge, >=, M1, C1, M2, C2) \ |
|
|
1194 SPARSE_SMM_CMP_OP (mx_el_gt, >, M1, C1, M2, C2) \ |
|
|
1195 SPARSE_SMM_CMP_OP (mx_el_eq, ==, M1, , M2, ) \ |
|
|
1196 SPARSE_SMM_CMP_OP (mx_el_ne, !=, M1, , M2, ) |
|
|
1197 |
|
|
1198 #define SPARSE_SMM_EQNE_OPS(M1, Z1, C1, M2, Z2, C2) \ |
|
|
1199 SPARSE_SMM_CMP_OP (mx_el_eq, ==, M1, , M2, ) \ |
|
|
1200 SPARSE_SMM_CMP_OP (mx_el_ne, !=, M1, , M2, ) |
|
|
1201 |
|
|
1202 #define SPARSE_SMM_BOOL_OP_DECLS(M1, M2) \ |
|
|
1203 SPARSE_BOOL_OP_DECL (mx_el_and, M1, M2); \ |
|
|
1204 SPARSE_BOOL_OP_DECL (mx_el_or, M1, M2); |
|
|
1205 |
|
|
1206 #define SPARSE_SMM_BOOL_OP(F, OP, M1, M2, LHS_ZERO, RHS_ZERO) \ |
|
|
1207 SparseBoolMatrix \ |
|
|
1208 F (const M1& m1, const M2& m2) \ |
|
|
1209 { \ |
|
|
1210 SparseBoolMatrix r; \ |
|
|
1211 \ |
|
|
1212 int m1_nr = m1.rows (); \ |
|
|
1213 int m1_nc = m1.cols (); \ |
|
|
1214 \ |
|
|
1215 int m2_nr = m2.rows (); \ |
|
|
1216 int m2_nc = m2.cols (); \ |
|
|
1217 \ |
|
|
1218 if (m1_nr == m2_nr && m1_nc == m2_nc) \ |
|
|
1219 { \ |
|
|
1220 if (m1_nr != 0 || m1_nc != 0) \ |
|
|
1221 { \ |
|
|
1222 /* Count num of non-zero elements */ \ |
|
|
1223 int nel = 0; \ |
|
|
1224 for (int j = 0; j < m1_nc; j++) \ |
|
|
1225 for (int i = 0; i < m1_nr; i++) \ |
|
|
1226 if ((m1.elem(i, j) != LHS_ZERO) \ |
|
|
1227 OP (m2.elem(i, j) != RHS_ZERO)) \ |
|
|
1228 nel++; \ |
|
|
1229 \ |
|
|
1230 r = SparseBoolMatrix (m1_nr, m1_nc, nel); \ |
|
|
1231 \ |
|
|
1232 int ii = 0; \ |
|
|
1233 r.cidx (0) = 0; \ |
|
|
1234 for (int j = 0; j < m1_nc; j++) \ |
|
|
1235 { \ |
|
|
1236 for (int i = 0; i < m1_nr; i++) \ |
|
|
1237 { \ |
|
|
1238 bool el = (m1.elem(i, j) != LHS_ZERO) \ |
|
|
1239 OP (m2.elem(i, j) != RHS_ZERO); \ |
|
|
1240 if (el) \ |
|
|
1241 { \ |
|
|
1242 r.data(ii) = el; \ |
|
|
1243 r.ridx(ii++) = i; \ |
|
|
1244 } \ |
|
|
1245 } \ |
|
|
1246 r.cidx(j+1) = ii; \ |
|
|
1247 } \ |
|
|
1248 } \ |
|
|
1249 } \ |
|
|
1250 else \ |
|
|
1251 { \ |
|
|
1252 if ((m1_nr != 0 || m1_nc != 0) && (m2_nr != 0 || m2_nc != 0)) \ |
|
|
1253 gripe_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \ |
|
|
1254 } \ |
|
|
1255 return r; \ |
|
|
1256 } |
|
|
1257 |
|
|
1258 #define SPARSE_SMM_BOOL_OPS2(M1, M2, LHS_ZERO, RHS_ZERO) \ |
|
|
1259 SPARSE_SMM_BOOL_OP (mx_el_and, &&, M1, M2, LHS_ZERO, RHS_ZERO) \ |
|
|
1260 SPARSE_SMM_BOOL_OP (mx_el_or, ||, M1, M2, LHS_ZERO, RHS_ZERO) \ |
|
|
1261 |
|
|
1262 #define SPARSE_SMM_BOOL_OPS(M1, M2, ZERO) \ |
|
|
1263 SPARSE_SMM_BOOL_OPS2(M1, M2, ZERO, ZERO) |
|
|
1264 |
|
|
1265 #define SPARSE_SMM_OP_DECLS(R1, R2, M1, M2) \ |
|
|
1266 SPARSE_SMM_BIN_OP_DECLS (R1, R2, M1, M2) \ |
|
|
1267 SPARSE_SMM_CMP_OP_DECLS (M1, M2) \ |
|
|
1268 SPARSE_SMM_BOOL_OP_DECLS (M1, M2) |
|
|
1269 |
|
|
1270 // Avoid some code duplication. Maybe we should use templates. |
|
|
1271 |
|
|
1272 #define SPARSE_CUMSUM(RET_TYPE, ELT_TYPE, FCN) \ |
|
|
1273 \ |
|
|
1274 int nr = rows (); \ |
|
|
1275 int nc = cols (); \ |
|
|
1276 \ |
|
|
1277 RET_TYPE retval; \ |
|
|
1278 \ |
|
|
1279 if (nr > 0 && nc > 0) \ |
|
|
1280 { \ |
|
|
1281 if ((nr == 1 && dim == -1) || dim == 1) \ |
|
|
1282 /* Ugly!! Is there a better way? */ \ |
|
|
1283 retval = transpose (). FCN (0) .transpose (); \ |
|
|
1284 else \ |
|
|
1285 { \ |
|
|
1286 int nel = 0; \ |
|
|
1287 for (int i = 0; i < nc; i++) \ |
|
|
1288 { \ |
|
|
1289 ELT_TYPE t = ELT_TYPE (); \ |
|
|
1290 for (int j = cidx (i); j < cidx (i+1); j++) \ |
|
|
1291 { \ |
|
|
1292 t += data(j); \ |
|
|
1293 if (t != ELT_TYPE ()) \ |
|
|
1294 if (j == cidx(i+1) - 1) \ |
|
|
1295 nel += nr - ridx(j); \ |
|
|
1296 else \ |
|
|
1297 nel += ridx(j+1) - ridx(j); \ |
|
|
1298 } \ |
|
|
1299 } \ |
|
|
1300 retval = RET_TYPE (nr, nc, nel); \ |
|
|
1301 retval.cidx(0) = 0; \ |
|
|
1302 int ii = 0; \ |
|
|
1303 for (int i = 0; i < nc; i++) \ |
|
|
1304 { \ |
|
|
1305 ELT_TYPE t = ELT_TYPE (); \ |
|
|
1306 for (int j = cidx (i); j < cidx (i+1); j++) \ |
|
|
1307 { \ |
|
|
1308 t += data(j); \ |
|
|
1309 if (t != ELT_TYPE ()) \ |
|
|
1310 { \ |
|
|
1311 if (j == cidx(i+1) - 1) \ |
|
|
1312 { \ |
|
|
1313 for (int k = ridx(j); k < nr; k++) \ |
|
|
1314 { \ |
|
|
1315 retval.data (ii) = t; \ |
|
|
1316 retval.ridx (ii++) = k; \ |
|
|
1317 } \ |
|
|
1318 } \ |
|
|
1319 else \ |
|
|
1320 { \ |
|
|
1321 for (int k = ridx(j); k < ridx(j+1); k++) \ |
|
|
1322 { \ |
|
|
1323 retval.data (ii) = t; \ |
|
|
1324 retval.ridx (ii++) = k; \ |
|
|
1325 } \ |
|
|
1326 } \ |
|
|
1327 } \ |
|
|
1328 } \ |
|
|
1329 retval.cidx(i+1) = ii; \ |
|
|
1330 } \ |
|
|
1331 } \ |
|
|
1332 } \ |
|
|
1333 else \ |
|
|
1334 retval = RET_TYPE (nr,nc); \ |
|
|
1335 \ |
|
|
1336 return retval |
|
|
1337 |
|
|
1338 |
|
|
1339 #define SPARSE_CUMPROD(RET_TYPE, ELT_TYPE, FCN) \ |
|
|
1340 \ |
|
|
1341 int nr = rows (); \ |
|
|
1342 int nc = cols (); \ |
|
|
1343 \ |
|
|
1344 RET_TYPE retval; \ |
|
|
1345 \ |
|
|
1346 if (nr > 0 && nc > 0) \ |
|
|
1347 { \ |
|
|
1348 if ((nr == 1 && dim == -1) || dim == 1) \ |
|
|
1349 /* Ugly!! Is there a better way? */ \ |
|
|
1350 retval = transpose (). FCN (0) .transpose (); \ |
|
|
1351 else \ |
|
|
1352 { \ |
|
|
1353 int nel = 0; \ |
|
|
1354 for (int i = 0; i < nc; i++) \ |
|
|
1355 { \ |
|
|
1356 int jj = 0; \ |
|
|
1357 for (int j = cidx (i); j < cidx (i+1); j++) \ |
|
|
1358 { \ |
|
|
1359 if (jj == ridx(j)) \ |
|
|
1360 { \ |
|
|
1361 nel++; \ |
|
|
1362 jj++; \ |
|
|
1363 } \ |
|
|
1364 else \ |
|
|
1365 break; \ |
|
|
1366 } \ |
|
|
1367 } \ |
|
|
1368 retval = RET_TYPE (nr, nc, nel); \ |
|
|
1369 retval.cidx(0) = 0; \ |
|
|
1370 int ii = 0; \ |
|
|
1371 for (int i = 0; i < nc; i++) \ |
|
|
1372 { \ |
|
|
1373 ELT_TYPE t = ELT_TYPE (1.); \ |
|
|
1374 int jj = 0; \ |
|
|
1375 for (int j = cidx (i); j < cidx (i+1); j++) \ |
|
|
1376 { \ |
|
|
1377 if (jj == ridx(j)) \ |
|
|
1378 { \ |
|
|
1379 t *= data(j); \ |
|
|
1380 retval.data(ii) = t; \ |
|
|
1381 retval.ridx(ii++) = jj++; \ |
|
|
1382 } \ |
|
|
1383 else \ |
|
|
1384 break; \ |
|
|
1385 } \ |
|
|
1386 retval.cidx(i+1) = ii; \ |
|
|
1387 } \ |
|
|
1388 } \ |
|
|
1389 } \ |
|
|
1390 else \ |
|
|
1391 retval = RET_TYPE (nr,nc); \ |
|
|
1392 \ |
|
|
1393 return retval |
|
|
1394 |
|
|
1395 #define SPARSE_BASE_REDUCTION_OP(RET_TYPE, EL_TYPE, ROW_EXPR, COL_EXPR, \ |
|
|
1396 INIT_VAL, MT_RESULT) \ |
|
|
1397 \ |
|
|
1398 int nr = rows (); \ |
|
|
1399 int nc = cols (); \ |
|
|
1400 \ |
|
|
1401 RET_TYPE retval; \ |
|
|
1402 \ |
|
|
1403 if (nr > 0 && nc > 0) \ |
|
|
1404 { \ |
|
|
1405 if ((nr == 1 && dim == -1) || dim == 1) \ |
|
|
1406 { \ |
|
|
1407 OCTAVE_LOCAL_BUFFER (EL_TYPE, tmp, nr); \ |
|
|
1408 \ |
|
|
1409 for (int i = 0; i < nr; i++) \ |
|
|
1410 { \ |
|
|
1411 tmp[i] = INIT_VAL; \ |
|
|
1412 for (int j = 0; j < nc; j++) \ |
|
|
1413 { \ |
|
|
1414 ROW_EXPR; \ |
|
|
1415 } \ |
|
|
1416 } \ |
|
|
1417 int nel = 0; \ |
|
|
1418 for (int i = 0; i < nr; i++) \ |
|
|
1419 if (tmp[i] != EL_TYPE ()) \ |
|
|
1420 nel++ ; \ |
|
|
1421 retval = RET_TYPE (nr, 1, nel); \ |
|
|
1422 retval.cidx(0) = 0; \ |
|
|
1423 retval.cidx(1) = nel; \ |
|
|
1424 nel = 0; \ |
|
|
1425 for (int i = 0; i < nr; i++) \ |
|
|
1426 if (tmp[i] != EL_TYPE ()) \ |
|
|
1427 { \ |
|
|
1428 retval.data(nel) = tmp[i]; \ |
|
|
1429 retval.ridx(nel++) = i; \ |
|
|
1430 } \ |
|
|
1431 } \ |
|
|
1432 else \ |
|
|
1433 { \ |
|
|
1434 OCTAVE_LOCAL_BUFFER (EL_TYPE, tmp, nc); \ |
|
|
1435 \ |
|
|
1436 for (int j = 0; j < nc; j++) \ |
|
|
1437 { \ |
|
|
1438 tmp[j] = INIT_VAL; \ |
|
|
1439 for (int i = 0; i < nr; i++) \ |
|
|
1440 { \ |
|
|
1441 COL_EXPR; \ |
|
|
1442 } \ |
|
|
1443 } \ |
|
|
1444 int nel = 0; \ |
|
|
1445 for (int i = 0; i < nc; i++) \ |
|
|
1446 if (tmp[i] != EL_TYPE ()) \ |
|
|
1447 nel++ ; \ |
|
|
1448 retval = RET_TYPE (1, nc, nel); \ |
|
|
1449 retval.cidx(0) = 0; \ |
|
|
1450 nel = 0; \ |
|
|
1451 for (int i = 0; i < nc; i++) \ |
|
|
1452 if (tmp[i] != EL_TYPE ()) \ |
|
|
1453 { \ |
|
|
1454 retval.data(nel) = tmp[i]; \ |
|
|
1455 retval.ridx(nel++) = 0; \ |
|
|
1456 retval.cidx(i+1) = retval.cidx(i) + 1; \ |
|
|
1457 } \ |
|
|
1458 else \ |
|
|
1459 retval.cidx(i+1) = retval.cidx(i); \ |
|
|
1460 } \ |
|
|
1461 } \ |
|
|
1462 else if (nc == 0 && (nr == 0 || (nr == 1 && dim == -1))) \ |
|
|
1463 { \ |
|
|
1464 retval = RET_TYPE (1, 1, 1); \ |
|
|
1465 retval.cidx(0) = 0; \ |
|
|
1466 retval.cidx(1) = 1; \ |
|
|
1467 retval.ridx(0) = 0; \ |
|
|
1468 retval.data(0) = MT_RESULT; \ |
|
|
1469 } \ |
|
|
1470 else if (nr == 0 && (dim == 0 || dim == -1)) \ |
|
|
1471 { \ |
|
|
1472 retval = RET_TYPE (1, nc, nc); \ |
|
|
1473 retval.cidx (0) = 0; \ |
|
|
1474 for (int i = 0; i < nc ; i++) \ |
|
|
1475 { \ |
|
|
1476 retval.ridx (i) = 0; \ |
|
|
1477 retval.cidx (i+1) = i; \ |
|
|
1478 retval.data (i) = MT_RESULT; \ |
|
|
1479 } \ |
|
|
1480 } \ |
|
|
1481 else if (nc == 0 && dim == 1) \ |
|
|
1482 { \ |
|
|
1483 retval = RET_TYPE (nr, 1, nr); \ |
|
|
1484 retval.cidx(0) = 0; \ |
|
|
1485 retval.cidx(1) = nr; \ |
|
|
1486 for (int i = 0; i < nr; i++) \ |
|
|
1487 { \ |
|
|
1488 retval.ridx(i) = i; \ |
|
|
1489 retval.data(i) = MT_RESULT; \ |
|
|
1490 } \ |
|
|
1491 } \ |
|
|
1492 else \ |
|
|
1493 retval.resize (nr > 0, nc > 0); \ |
|
|
1494 \ |
|
|
1495 return retval |
|
|
1496 |
|
|
1497 #define SPARSE_REDUCTION_OP_ROW_EXPR(OP) \ |
|
|
1498 tmp[i] OP elem (i, j) |
|
|
1499 |
|
|
1500 #define SPARSE_REDUCTION_OP_COL_EXPR(OP) \ |
|
|
1501 tmp[j] OP elem (i, j) |
|
|
1502 |
|
|
1503 #define SPARSE_REDUCTION_OP(RET_TYPE, EL_TYPE, OP, INIT_VAL, MT_RESULT) \ |
|
|
1504 SPARSE_BASE_REDUCTION_OP (RET_TYPE, EL_TYPE, \ |
|
|
1505 SPARSE_REDUCTION_OP_ROW_EXPR (OP), \ |
|
|
1506 SPARSE_REDUCTION_OP_COL_EXPR (OP), \ |
|
|
1507 INIT_VAL, MT_RESULT) |
|
|
1508 |
|
|
1509 #define SPARSE_ANY_ALL_OP_ROW_CODE(TEST_OP, TEST_TRUE_VAL) \ |
|
|
1510 if (elem (i, j) TEST_OP 0.0) \ |
|
|
1511 { \ |
|
|
1512 tmp[i] = TEST_TRUE_VAL; \ |
|
|
1513 break; \ |
|
|
1514 } |
|
|
1515 |
|
|
1516 #define SPARSE_ANY_ALL_OP_COL_CODE(TEST_OP, TEST_TRUE_VAL) \ |
|
|
1517 if (elem (i, j) TEST_OP 0.0) \ |
|
|
1518 { \ |
|
|
1519 tmp[j] = TEST_TRUE_VAL; \ |
|
|
1520 break; \ |
|
|
1521 } |
|
|
1522 |
|
|
1523 #define SPARSE_ANY_ALL_OP(DIM, INIT_VAL, TEST_OP, TEST_TRUE_VAL) \ |
|
|
1524 SPARSE_BASE_REDUCTION_OP (SparseBoolMatrix, char, \ |
|
|
1525 SPARSE_ANY_ALL_OP_ROW_CODE (TEST_OP, TEST_TRUE_VAL), \ |
|
|
1526 SPARSE_ANY_ALL_OP_COL_CODE (TEST_OP, TEST_TRUE_VAL), \ |
|
|
1527 INIT_VAL, INIT_VAL) |
|
|
1528 |
|
|
1529 #define SPARSE_ALL_OP(DIM) SPARSE_ANY_ALL_OP (DIM, true, ==, false) |
|
|
1530 |
|
|
1531 #define SPARSE_ANY_OP(DIM) SPARSE_ANY_ALL_OP (DIM, false, !=, true) |
|
|
1532 |
|
|
1533 #define SPARSE_SPARSE_MUL( RET_TYPE, EL_TYPE ) \ |
|
|
1534 int nr = m.rows (); \ |
|
|
1535 int nc = m.cols (); \ |
|
|
1536 \ |
|
|
1537 int a_nr = a.rows (); \ |
|
|
1538 int a_nc = a.cols (); \ |
|
|
1539 \ |
|
|
1540 if (nc != a_nr) \ |
|
|
1541 { \ |
|
|
1542 gripe_nonconformant ("operator *", nr, nc, a_nr, a_nc); \ |
|
|
1543 return RET_TYPE (); \ |
|
|
1544 } \ |
|
|
1545 else \ |
|
|
1546 { \ |
|
|
1547 OCTAVE_LOCAL_BUFFER (EL_TYPE, Xcol, nr); \ |
|
|
1548 \ |
|
|
1549 int nel = 0; \ |
|
|
1550 \ |
|
|
1551 for (int i = 0; i < a_nc; i++) \ |
|
|
1552 { \ |
|
|
1553 OCTAVE_QUIT; \ |
|
|
1554 for (int k = 0; k < nr; k++) \ |
|
|
1555 Xcol[k]= 0.; \ |
|
|
1556 for (int j = a.cidx(i); j < a.cidx(i+1); j++) \ |
|
|
1557 { \ |
|
|
1558 int col = a.ridx(j); \ |
|
|
1559 for (int k = m.cidx(col) ; k < m.cidx(col+1); k++) \ |
|
|
1560 if (Xcol[m.ridx(k)] == 0.) \ |
|
|
1561 { \ |
|
|
1562 Xcol[m.ridx(k)] = 1.; \ |
|
|
1563 nel++; \ |
|
|
1564 } \ |
|
|
1565 } \ |
|
|
1566 } \ |
|
|
1567 \ |
|
|
1568 if (nel == 0) \ |
|
|
1569 return RET_TYPE (nr, a_nc); \ |
|
|
1570 else \ |
|
|
1571 { \ |
|
|
1572 RET_TYPE retval (nr, a_nc, nel); \ |
|
|
1573 \ |
|
|
1574 int ii = 0; \ |
|
|
1575 \ |
|
|
1576 retval.cidx(0) = 0; \ |
|
|
1577 for (int i = 0; i < a_nc ; i++) \ |
|
|
1578 { \ |
|
|
1579 OCTAVE_QUIT; \ |
|
|
1580 for (int k = 0; k < nr; k++) \ |
|
|
1581 Xcol[k]= 0.; \ |
|
|
1582 for (int j = a.cidx(i); j < a.cidx(i+1); j++) \ |
|
|
1583 { \ |
|
|
1584 int col = a.ridx(j); \ |
|
|
1585 EL_TYPE tmpval = a.data(j); \ |
|
|
1586 for (int k = m.cidx(col) ; k < m.cidx(col+1); k++) \ |
|
|
1587 Xcol[m.ridx(k)] += tmpval * m.data(k); \ |
|
|
1588 } \ |
|
|
1589 for (int k = 0; k < nr; k++) \ |
|
|
1590 { \ |
|
|
1591 if (Xcol[k] !=0. ) \ |
|
|
1592 { \ |
|
|
1593 retval.ridx (ii) = k; \ |
|
|
1594 retval.data (ii++) = Xcol[k]; \ |
|
|
1595 } \ |
|
|
1596 } \ |
|
|
1597 retval.cidx(i+1) = ii; \ |
|
|
1598 } \ |
|
|
1599 return retval; \ |
|
|
1600 } \ |
|
|
1601 } |
|
|
1602 |
|
|
1603 #endif |
|
|
1604 |
|
|
1605 /* |
|
|
1606 ;;; Local Variables: *** |
|
|
1607 ;;; mode: C++ *** |
|
|
1608 ;;; End: *** |
|
|
1609 */ |
