|
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 { \ |
|
|
1907 RET_TYPE retval (a_nr, a_nc, ZERO); \ |
|
|
1908 for (octave_idx_type i = 0; i < a_nc ; i++) \ |
|
|
1909 { \ |
|
|
1910 for (octave_idx_type j = 0; j < a_nr; j++) \ |
|
|
1911 { \ |
|
|
1912 OCTAVE_QUIT; \ |
|
|
1913 retval.elem (j,i) += a.elem(j,i) * m.elem(0,0); \ |
|
|
1914 } \ |
|
|
1915 } \ |
|
|
1916 return retval; \ |
|
|
1917 } \ |
|
|
1918 else if (nc != a_nr) \ |
|
5429
|
1919 { \ |
|
|
1920 gripe_nonconformant ("operator *", nr, nc, a_nr, a_nc); \ |
|
|
1921 return RET_TYPE (); \ |
|
|
1922 } \ |
|
|
1923 else \ |
|
|
1924 { \ |
|
5681
|
1925 RET_TYPE retval (nr, a_nc, ZERO); \ |
|
5429
|
1926 \ |
|
|
1927 for (octave_idx_type i = 0; i < a_nc ; i++) \ |
|
|
1928 { \ |
|
|
1929 for (octave_idx_type j = 0; j < a_nr; j++) \ |
|
|
1930 { \ |
|
|
1931 OCTAVE_QUIT; \ |
|
|
1932 \ |
|
|
1933 EL_TYPE tmpval = a.elem(j,i); \ |
|
|
1934 for (octave_idx_type k = m.cidx(j) ; k < m.cidx(j+1); k++) \ |
|
|
1935 retval.elem (m.ridx(k),i) += tmpval * m.data(k); \ |
|
|
1936 } \ |
|
|
1937 } \ |
|
|
1938 return retval; \ |
|
|
1939 } |
|
|
1940 |
|
5681
|
1941 #define FULL_SPARSE_MUL( RET_TYPE, EL_TYPE, ZERO ) \ |
|
5429
|
1942 octave_idx_type nr = m.rows (); \ |
|
|
1943 octave_idx_type nc = m.cols (); \ |
|
|
1944 \ |
|
|
1945 octave_idx_type a_nr = a.rows (); \ |
|
|
1946 octave_idx_type a_nc = a.cols (); \ |
|
|
1947 \ |
|
6221
|
1948 if (a_nr == 1 && a_nc == 1) \ |
|
|
1949 { \ |
|
|
1950 RET_TYPE retval (nr, nc, ZERO); \ |
|
|
1951 for (octave_idx_type i = 0; i < nc ; i++) \ |
|
|
1952 { \ |
|
|
1953 for (octave_idx_type j = 0; j < nr; j++) \ |
|
|
1954 { \ |
|
|
1955 OCTAVE_QUIT; \ |
|
|
1956 retval.elem (j,i) += a.elem(0,0) * m.elem(j,i); \ |
|
|
1957 } \ |
|
|
1958 } \ |
|
|
1959 return retval; \ |
|
|
1960 } \ |
|
|
1961 else if (nc != a_nr) \ |
|
5429
|
1962 { \ |
|
|
1963 gripe_nonconformant ("operator *", nr, nc, a_nr, a_nc); \ |
|
|
1964 return RET_TYPE (); \ |
|
|
1965 } \ |
|
|
1966 else \ |
|
|
1967 { \ |
|
5681
|
1968 RET_TYPE retval (nr, a_nc, ZERO); \ |
|
5429
|
1969 \ |
|
|
1970 for (octave_idx_type i = 0; i < a_nc ; i++) \ |
|
|
1971 { \ |
|
|
1972 for (octave_idx_type j = a.cidx(i); j < a.cidx(i+1); j++) \ |
|
|
1973 { \ |
|
|
1974 octave_idx_type col = a.ridx(j); \ |
|
|
1975 EL_TYPE tmpval = a.data(j); \ |
|
|
1976 OCTAVE_QUIT; \ |
|
|
1977 \ |
|
|
1978 for (octave_idx_type k = 0 ; k < nr; k++) \ |
|
|
1979 retval.elem (k,i) += tmpval * m.elem(k,col); \ |
|
|
1980 } \ |
|
|
1981 } \ |
|
|
1982 return retval; \ |
|
|
1983 } |
|
|
1984 |
|
5164
|
1985 #endif |
|
|
1986 |
|
|
1987 /* |
|
|
1988 ;;; Local Variables: *** |
|
|
1989 ;;; mode: C++ *** |
|
|
1990 ;;; End: *** |
|
|
1991 */ |
