GNU Octave 11.1.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 
Loading...
Searching...
No Matches
Sparse-op-defs.h
Go to the documentation of this file.
1////////////////////////////////////////////////////////////////////////
2//
3// Copyright (C) 1998-2026 The Octave Project Developers
4//
5// See the file COPYRIGHT.md in the top-level directory of this
6// distribution or <https://octave.org/copyright/>.
7//
8// This file is part of Octave.
9//
10// Octave is free software: you can redistribute it and/or modify it
11// under the terms of the GNU General Public License as published by
12// the Free Software Foundation, either version 3 of the License, or
13// (at your option) any later version.
14//
15// Octave is distributed in the hope that it will be useful, but
16// WITHOUT ANY WARRANTY; without even the implied warranty of
17// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18// GNU General Public License for more details.
19//
20// You should have received a copy of the GNU General Public License
21// along with Octave; see the file COPYING. If not, see
22// <https://www.gnu.org/licenses/>.
23//
24////////////////////////////////////////////////////////////////////////
25
26#if ! defined (octave_Sparse_op_defs_h)
27#define octave_Sparse_op_defs_h 1
28
29#include "octave-config.h"
30
31#include "Array-util.h"
32#include "lo-array-errwarn.h"
33#include "mx-inlines.cc"
34#include "oct-locbuf.h"
35
36// sparse matrix by scalar operations.
37
38#define SPARSE_SMS_BIN_OP_1(R, F, OP, M, S) \
39 R \
40 F (const M& m, const S& s) \
41 { \
42 octave_idx_type nr = m.rows (); \
43 octave_idx_type nc = m.cols (); \
44 \
45 R r (nr, nc, (0.0 OP s)); \
46 \
47 for (octave_idx_type j = 0; j < nc; j++) \
48 for (octave_idx_type i = m.cidx (j); i < m.cidx (j+1); i++) \
49 r.xelem (m.ridx (i), j) = m.data (i) OP s; \
50 return r; \
51 }
52
53#define SPARSE_SMS_BIN_OP_2(R, F, OP, M, S) \
54 R \
55 F (const M& m, const S& s) \
56 { \
57 octave_idx_type nr = m.rows (); \
58 octave_idx_type nc = m.cols (); \
59 octave_idx_type nz = m.nnz (); \
60 \
61 R r (nr, nc, nz); \
62 \
63 for (octave_idx_type i = 0; i < nz; i++) \
64 { \
65 r.xdata (i) = m.data (i) OP s; \
66 r.xridx (i) = m.ridx (i); \
67 } \
68 for (octave_idx_type i = 0; i < nc + 1; i++) \
69 r.xcidx (i) = m.cidx (i); \
70 \
71 r.maybe_compress (true); \
72 return r; \
73 }
74
75#define SPARSE_SMS_BIN_OPS(R1, R2, M, S) \
76 SPARSE_SMS_BIN_OP_1 (R1, operator +, +, M, S) \
77 SPARSE_SMS_BIN_OP_1 (R1, operator -, -, M, S) \
78 SPARSE_SMS_BIN_OP_2 (R2, operator *, *, M, S) \
79 SPARSE_SMS_BIN_OP_2 (R2, operator /, /, M, S)
80
81#define SPARSE_SMS_CMP_OP(F, OP, M, S) \
82 SparseBoolMatrix \
83 F (const M& m, const S& s) \
84 { \
85 octave_idx_type nr = m.rows (); \
86 octave_idx_type nc = m.cols (); \
87 SparseBoolMatrix r; \
88 \
89 M::element_type m_zero = M::element_type (); \
90 \
91 if (m_zero OP s) \
92 { \
93 r = SparseBoolMatrix (nr, nc, true); \
94 for (octave_idx_type j = 0; j < nc; j++) \
95 for (octave_idx_type i = m.cidx (j); i < m.cidx (j+1); i++) \
96 if (! (m.data (i) OP s)) \
97 r.data (m.ridx (i) + j * nr) = false; \
98 r.maybe_compress (true); \
99 } \
100 else \
101 { \
102 r = SparseBoolMatrix (nr, nc, m.nnz ()); \
103 r.cidx (0) = static_cast<octave_idx_type> (0); \
104 octave_idx_type nel = 0; \
105 for (octave_idx_type j = 0; j < nc; j++) \
106 { \
107 for (octave_idx_type i = m.cidx (j); i < m.cidx (j+1); i++) \
108 if (m.data (i) OP s) \
109 { \
110 r.ridx (nel) = m.ridx (i); \
111 r.data (nel++) = true; \
112 } \
113 r.cidx (j + 1) = nel; \
114 } \
115 r.maybe_compress (false); \
116 } \
117 return r; \
118 }
119
120#define SPARSE_SMS_CMP_OPS(M, S) \
121 SPARSE_SMS_CMP_OP (mx_el_lt, <, M, S) \
122 SPARSE_SMS_CMP_OP (mx_el_le, <=, M, S) \
123 SPARSE_SMS_CMP_OP (mx_el_ge, >=, M, S) \
124 SPARSE_SMS_CMP_OP (mx_el_gt, >, M, S) \
125 SPARSE_SMS_CMP_OP (mx_el_eq, ==, M, S) \
126 SPARSE_SMS_CMP_OP (mx_el_ne, !=, M, S)
127
128#define SPARSE_SMS_EQNE_OPS(M, S) \
129 SPARSE_SMS_CMP_OP (mx_el_eq, ==, M, S) \
130 SPARSE_SMS_CMP_OP (mx_el_ne, !=, M, S)
131
132#define SPARSE_SMS_BOOL_AND_OP(M, S) \
133 SparseBoolMatrix \
134 mx_el_and (const M& m, const S& s) \
135 { \
136 octave_idx_type nr = m.rows (); \
137 octave_idx_type nc = m.cols (); \
138 SparseBoolMatrix r; \
139 \
140 M::element_type lhs_zero = M::element_type (); \
141 S rhs_zero = S (); \
142 \
143 if (nr > 0 && nc > 0) \
144 { \
145 if (s != rhs_zero) \
146 { \
147 r = SparseBoolMatrix (nr, nc, m.nnz ()); \
148 r.cidx (0) = static_cast<octave_idx_type> (0); \
149 octave_idx_type nel = 0; \
150 for (octave_idx_type j = 0; j < nc; j++) \
151 { \
152 for (octave_idx_type i = m.cidx (j); i < m.cidx (j+1); i++) \
153 if (m.data (i) != lhs_zero) \
154 { \
155 r.ridx (nel) = m.ridx (i); \
156 r.data (nel++) = true; \
157 } \
158 r.cidx (j + 1) = nel; \
159 } \
160 r.maybe_compress (false); \
161 } \
162 else \
163 r = SparseBoolMatrix (nr, nc); \
164 } \
165 else \
166 r = SparseBoolMatrix (nr, nc); \
167 return r; \
168 }
169
170#define SPARSE_SMS_BOOL_OR_OP(M, S) \
171 boolMatrix \
172 mx_el_or (const M& m, const S& s) \
173 { \
174 octave_idx_type nr = m.rows (); \
175 octave_idx_type nc = m.cols (); \
176 boolMatrix r; \
177 \
178 M::element_type lhs_zero = M::element_type (); \
179 S rhs_zero = S (); \
180 \
181 if (nr > 0 && nc > 0) \
182 { \
183 if (s != rhs_zero) \
184 r = boolMatrix (nr, nc, true); \
185 else \
186 { \
187 r = boolMatrix (nr, nc, false); \
188 \
189 for (octave_idx_type j = 0; j < nc; j++) \
190 for (octave_idx_type i = 0; i < nr; i++) \
191 if (m.elem (i, j) != lhs_zero) \
192 r.elem (i, j) = true; \
193 } \
194 } \
195 else \
196 r = boolMatrix (nr, nc, false); \
197 return r; \
198 }
199
200#define SPARSE_SMS_BOOL_OPS(M, S) \
201 SPARSE_SMS_BOOL_AND_OP (M, S) \
202 SPARSE_SMS_BOOL_OR_OP (M, S)
203
204// scalar by sparse matrix operations.
205
206#define SPARSE_SSM_BIN_OP_1(R, F, OP, S, M) \
207 R \
208 F (const S& s, const M& m) \
209 { \
210 octave_idx_type nr = m.rows (); \
211 octave_idx_type nc = m.cols (); \
212 \
213 R r (nr, nc, (s OP 0.0)); \
214 \
215 for (octave_idx_type j = 0; j < nc; j++) \
216 for (octave_idx_type i = m.cidx (j); i < m.cidx (j+1); i++) \
217 r.xelem (m.ridx (i), j) = s OP m.data (i); \
218 \
219 return r; \
220 }
221
222#define SPARSE_SSM_BIN_OP_2(R, F, OP, S, M) \
223 R \
224 F (const S& s, const M& m) \
225 { \
226 octave_idx_type nr = m.rows (); \
227 octave_idx_type nc = m.cols (); \
228 octave_idx_type nz = m.nnz (); \
229 \
230 R r (nr, nc, nz); \
231 \
232 for (octave_idx_type i = 0; i < nz; i++) \
233 { \
234 r.xdata (i) = s OP m.data (i); \
235 r.xridx (i) = m.ridx (i); \
236 } \
237 for (octave_idx_type i = 0; i < nc + 1; i++) \
238 r.xcidx (i) = m.cidx (i); \
239 \
240 r.maybe_compress(true); \
241 return r; \
242 }
243
244#define SPARSE_SSM_BIN_OPS(R1, R2, S, M) \
245 SPARSE_SSM_BIN_OP_1 (R1, operator +, +, S, M) \
246 SPARSE_SSM_BIN_OP_1 (R1, operator -, -, S, M) \
247 SPARSE_SSM_BIN_OP_2 (R2, operator *, *, S, M) \
248 SPARSE_SSM_BIN_OP_1 (R1, operator /, /, S, M)
249
250#define SPARSE_SSM_CMP_OP(F, OP, S, M) \
251 SparseBoolMatrix \
252 F (const S& s, const M& m) \
253 { \
254 octave_idx_type nr = m.rows (); \
255 octave_idx_type nc = m.cols (); \
256 SparseBoolMatrix r; \
257 \
258 M::element_type m_zero = M::element_type (); \
259 \
260 if (s OP m_zero) \
261 { \
262 r = SparseBoolMatrix (nr, nc, true); \
263 for (octave_idx_type j = 0; j < nc; j++) \
264 for (octave_idx_type i = m.cidx (j); i < m.cidx (j+1); i++) \
265 if (! (s OP m.data (i))) \
266 r.data (m.ridx (i) + j * nr) = false; \
267 r.maybe_compress (true); \
268 } \
269 else \
270 { \
271 r = SparseBoolMatrix (nr, nc, m.nnz ()); \
272 r.cidx (0) = static_cast<octave_idx_type> (0); \
273 octave_idx_type nel = 0; \
274 for (octave_idx_type j = 0; j < nc; j++) \
275 { \
276 for (octave_idx_type i = m.cidx (j); i < m.cidx (j+1); i++) \
277 if (s OP m.data (i)) \
278 { \
279 r.ridx (nel) = m.ridx (i); \
280 r.data (nel++) = true; \
281 } \
282 r.cidx (j + 1) = nel; \
283 } \
284 r.maybe_compress (false); \
285 } \
286 return r; \
287 }
288
289#define SPARSE_SSM_CMP_OPS(S, M) \
290 SPARSE_SSM_CMP_OP (mx_el_lt, <, S, M) \
291 SPARSE_SSM_CMP_OP (mx_el_le, <=, S, M) \
292 SPARSE_SSM_CMP_OP (mx_el_ge, >=, S, M) \
293 SPARSE_SSM_CMP_OP (mx_el_gt, >, S, M) \
294 SPARSE_SSM_CMP_OP (mx_el_eq, ==, S, M) \
295 SPARSE_SSM_CMP_OP (mx_el_ne, !=, S, M)
296
297#define SPARSE_SSM_EQNE_OPS(S, M) \
298 SPARSE_SSM_CMP_OP (mx_el_eq, ==, S, M) \
299 SPARSE_SSM_CMP_OP (mx_el_ne, !=, S, M)
300
301#define SPARSE_SSM_BOOL_AND_OP(S, M) \
302 SparseBoolMatrix \
303 mx_el_and (const S& s, const M& m) \
304 { \
305 octave_idx_type nr = m.rows (); \
306 octave_idx_type nc = m.cols (); \
307 SparseBoolMatrix r; \
308 \
309 S lhs_zero = S (); \
310 M::element_type rhs_zero = M::element_type (); \
311 \
312 if (nr > 0 && nc > 0) \
313 { \
314 if (s != lhs_zero) \
315 { \
316 r = SparseBoolMatrix (nr, nc, m.nnz ()); \
317 r.cidx (0) = static_cast<octave_idx_type> (0); \
318 octave_idx_type nel = 0; \
319 for (octave_idx_type j = 0; j < nc; j++) \
320 { \
321 for (octave_idx_type i = m.cidx (j); i < m.cidx (j+1); i++) \
322 if (m.data (i) != rhs_zero) \
323 { \
324 r.ridx (nel) = m.ridx (i); \
325 r.data (nel++) = true; \
326 } \
327 r.cidx (j + 1) = nel; \
328 } \
329 r.maybe_compress (false); \
330 } \
331 else \
332 r = SparseBoolMatrix (nr, nc); \
333 } \
334 else \
335 r = SparseBoolMatrix (nr, nc); \
336 return r; \
337 }
338
339#define SPARSE_SSM_BOOL_OR_OP(S, M) \
340 boolMatrix \
341 mx_el_or (const S& s, const M& m) \
342 { \
343 octave_idx_type nr = m.rows (); \
344 octave_idx_type nc = m.cols (); \
345 boolMatrix r; \
346 \
347 S lhs_zero = S (); \
348 M::element_type rhs_zero = M::element_type (); \
349 \
350 if (nr > 0 && nc > 0) \
351 { \
352 if (s != lhs_zero) \
353 r = boolMatrix (nr, nc, true); \
354 else \
355 { \
356 r = boolMatrix (nr, nc, false); \
357 \
358 for (octave_idx_type j = 0; j < nc; j++) \
359 for (octave_idx_type i = 0; i < nr; i++) \
360 if (m.elem (i, j) != rhs_zero) \
361 r.elem (i, j) = true; \
362 } \
363 } \
364 else \
365 r = boolMatrix (nr, nc, false); \
366 return r; \
367 }
368
369#define SPARSE_SSM_BOOL_OPS(S, M) \
370 SPARSE_SSM_BOOL_AND_OP (S, M) \
371 SPARSE_SSM_BOOL_OR_OP (S, M)
372
373// sparse matrix by sparse matrix operations.
374
375#define SPARSE_SMSM_BIN_OP_1(R, F, OP, M1, M2) \
376 R \
377 F (const M1& m1, const M2& m2) \
378 { \
379 R r; \
380 \
381 octave_idx_type m1_nr = m1.rows (); \
382 octave_idx_type m1_nc = m1.cols (); \
383 \
384 octave_idx_type m2_nr = m2.rows (); \
385 octave_idx_type m2_nc = m2.cols (); \
386 \
387 if (m1_nr == 1 && m1_nc == 1) \
388 { \
389 if (m1.elem (0,0) == 0.) \
390 r = OP R (m2); \
391 else \
392 { \
393 r = R (m2_nr, m2_nc, m1.data (0) OP 0.); \
394 \
395 for (octave_idx_type j = 0 ; j < m2_nc ; j++) \
396 { \
397 octave_quit (); \
398 octave_idx_type idxj = j * m2_nr; \
399 for (octave_idx_type i = m2.cidx (j) ; i < m2.cidx (j+1) ; i++) \
400 { \
401 octave_quit (); \
402 r.data (idxj + m2.ridx (i)) = m1.data (0) OP m2.data (i); \
403 } \
404 } \
405 r.maybe_compress (); \
406 } \
407 } \
408 else if (m2_nr == 1 && m2_nc == 1) \
409 { \
410 if (m2.elem (0,0) == 0.) \
411 r = R (m1); \
412 else \
413 { \
414 r = R (m1_nr, m1_nc, 0. OP m2.data (0)); \
415 \
416 for (octave_idx_type j = 0 ; j < m1_nc ; j++) \
417 { \
418 octave_quit (); \
419 octave_idx_type idxj = j * m1_nr; \
420 for (octave_idx_type i = m1.cidx (j) ; i < m1.cidx (j+1) ; i++) \
421 { \
422 octave_quit (); \
423 r.data (idxj + m1.ridx (i)) = m1.data (i) OP m2.data (0); \
424 } \
425 } \
426 r.maybe_compress (); \
427 } \
428 } \
429 else if (m1_nr == 0 && (m2_nr == 0 || m2_nr == 1)) \
430 { \
431 if (m1_nc == 1 || m2_nc == 1 || m1_nc == m2_nc) \
432 r.resize (m1_nr, std::max (m1_nc, m2_nc)); \
433 else \
434 octave::err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \
435 } \
436 else if (m1_nc == 0 && (m2_nc == 0 || m2_nc == 1)) \
437 { \
438 if (m1_nr == 1 || m2_nr == 1 || m1_nr == m2_nr) \
439 r.resize (std::max (m1_nr, m2_nr), m1_nc); \
440 else \
441 octave::err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \
442 } \
443 else if (m2_nr == 0 && (m1_nr == 0 || m1_nr == 1)) \
444 { \
445 if (m2_nc == 1 || m1_nc == 1 || m2_nc == m1_nc) \
446 r.resize (m2_nr, std::max (m1_nc, m2_nc)); \
447 else \
448 octave::err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \
449 } \
450 else if (m2_nc == 0 && (m1_nc == 0 || m1_nc == 1)) \
451 { \
452 if (m2_nr == 1 || m1_nr == 1 || m2_nr == m1_nr) \
453 r.resize (std::max (m1_nr, m2_nr), m2_nc); \
454 else \
455 octave::err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \
456 } \
457 else if (m1_nr == m2_nr && m1_nc == m2_nc) \
458 { \
459 r = R (m1_nr, m1_nc, (m1.nnz () + m2.nnz ())); \
460 \
461 octave_idx_type jx = 0; \
462 r.cidx (0) = 0; \
463 for (octave_idx_type i = 0 ; i < m1_nc ; i++) \
464 { \
465 octave_idx_type ja = m1.cidx (i); \
466 octave_idx_type ja_max = m1.cidx (i+1); \
467 bool ja_lt_max = ja < ja_max; \
468 \
469 octave_idx_type jb = m2.cidx (i); \
470 octave_idx_type jb_max = m2.cidx (i+1); \
471 bool jb_lt_max = jb < jb_max; \
472 \
473 while (ja_lt_max || jb_lt_max) \
474 { \
475 octave_quit (); \
476 if ((! jb_lt_max) || \
477 (ja_lt_max && (m1.ridx (ja) < m2.ridx (jb)))) \
478 { \
479 r.ridx (jx) = m1.ridx (ja); \
480 r.data (jx) = m1.data (ja) OP 0.; \
481 jx++; \
482 ja++; \
483 ja_lt_max= ja < ja_max; \
484 } \
485 else if ((! ja_lt_max) || \
486 (jb_lt_max && (m2.ridx (jb) < m1.ridx (ja)))) \
487 { \
488 r.ridx (jx) = m2.ridx (jb); \
489 r.data (jx) = 0. OP m2.data (jb); \
490 jx++; \
491 jb++; \
492 jb_lt_max= jb < jb_max; \
493 } \
494 else \
495 { \
496 if ((m1.data (ja) OP m2.data (jb)) != 0.) \
497 { \
498 r.data (jx) = m1.data (ja) OP m2.data (jb); \
499 r.ridx (jx) = m1.ridx (ja); \
500 jx++; \
501 } \
502 ja++; \
503 ja_lt_max= ja < ja_max; \
504 jb++; \
505 jb_lt_max= jb < jb_max; \
506 } \
507 } \
508 r.cidx (i+1) = jx; \
509 } \
510 \
511 r.maybe_compress (); \
512 } \
513 else if (m1_nr == m2_nr && (m1_nc == 1 || m2_nc == 1)) \
514 { \
515 octave_idx_type r_nc = (m1_nc < m2_nc ? m2_nc : m1_nc); \
516 octave_idx_type rnnz = (m1_nc < m2_nc ? m1.nnz () * m2_nc + m2.nnz () : \
517 m1.nnz () + m1_nc * m2.nnz ()); \
518 r = R (m1_nr, r_nc, rnnz); \
519 \
520 octave_idx_type jx = 0; \
521 r.cidx (0) = 0; \
522 for (octave_idx_type i = 0 ; i < r_nc ; i++) \
523 { \
524 octave_idx_type ja; \
525 octave_idx_type ja_max; \
526 octave_idx_type jb; \
527 octave_idx_type jb_max; \
528 if (m1_nc == 1) \
529 { \
530 ja = m1.cidx(0); \
531 ja_max = m1.cidx(1); \
532 jb = m2.cidx (i); \
533 jb_max = m2.cidx (i+1); \
534 } \
535 else \
536 { \
537 ja = m1.cidx(i); \
538 ja_max = m1.cidx(i+1); \
539 jb = m2.cidx (0); \
540 jb_max = m2.cidx (1); \
541 } \
542 bool ja_lt_max = ja < ja_max; \
543 bool jb_lt_max = jb < jb_max; \
544 \
545 while (ja_lt_max || jb_lt_max) \
546 { \
547 octave_quit (); \
548 if ((! jb_lt_max) || \
549 (ja_lt_max && (m1.ridx (ja) < m2.ridx (jb)))) \
550 { \
551 r.ridx (jx) = m1.ridx (ja); \
552 r.data (jx) = m1.data (ja) OP 0.; \
553 jx++; \
554 ja++; \
555 ja_lt_max = ja < ja_max; \
556 } \
557 else if ((! ja_lt_max) || \
558 (m2.ridx (jb) < m1.ridx (ja))) \
559 { \
560 r.ridx (jx) = m2.ridx (jb); \
561 r.data (jx) = 0. OP m2.data (jb); \
562 jx++; \
563 jb++; \
564 jb_lt_max = jb < jb_max; \
565 } \
566 else \
567 { \
568 if (m1.data (ja) OP m2.data (jb) != 0.) \
569 { \
570 r.data (jx) = m1.data (ja) OP m2.data (jb); \
571 r.ridx (jx) = m1.ridx (ja); \
572 jx++; \
573 } \
574 ja++; \
575 ja_lt_max = ja < ja_max; \
576 jb++; \
577 jb_lt_max = jb < jb_max; \
578 } \
579 } \
580 r.cidx (i+1) = jx; \
581 } \
582 r.maybe_compress (); \
583 } \
584 else if (m1_nc == m2_nc && (m1_nr == 1 || m2_nr == 1)) \
585 r = F (m1.transpose (), m2.transpose ()).transpose (); \
586 else if (m1_nr == 1 && m2_nc == 1) \
587 { \
588 r = R (m2_nr, m1_nc, (m1.nnz () * m2_nr + m2.nnz () * m1_nc)); \
589 \
590 for (octave_idx_type j = 0; j < m1_nc; j++) \
591 { \
592 const Complex a_val = m1.elem (0, j); \
593 \
594 for (octave_idx_type i = 0; i < m2_nr; i++) \
595 { \
596 octave_quit (); \
597 \
598 const Complex b_val = m2.elem (i, 0); \
599 const Complex val = a_val OP b_val; \
600 \
601 if (val != Complex ()) \
602 r.elem (i, j) = val; \
603 } \
604 } \
605 r.maybe_compress (); \
606 } \
607 else if (m1_nc == 1 && m2_nr == 1) \
608 r = F (m1.transpose (), m2.transpose ()).transpose (); \
609 else \
610 octave::err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \
611 \
612 return r; \
613 }
614
615#define SPARSE_SMSM_BIN_OP_2(R, F, OP, M1, M2) \
616 R \
617 F (const M1& m1, const M2& m2) \
618 { \
619 R r; \
620 Complex nan = Complex (octave::numeric_limits<double>::NaN (), \
621 octave::numeric_limits<double>::NaN ()); \
622 \
623 octave_idx_type m1_nr = m1.rows (); \
624 octave_idx_type m1_nc = m1.cols (); \
625 \
626 octave_idx_type m2_nr = m2.rows (); \
627 octave_idx_type m2_nc = m2.cols (); \
628 \
629 if (m1_nr == 1 && m1_nc == 1) \
630 { \
631 if (m1.elem (0, 0) == 0.0) \
632 r = R (m2_nr, m2_nc); \
633 else if (octave::math::isnan (m1.elem (0, 0))) \
634 { \
635 r = R (m2_nr, m2_nc, m2_nr * m2_nc); \
636 for (octave_idx_type i = 0 ; i < r.numel () ; i++) \
637 r.elem(i) = nan; \
638 } \
639 else if (octave::math::isinf (m1.elem (0, 0))) \
640 { \
641 r = R (m2_nr, m2_nc, m2_nr * m2_nc); \
642 \
643 for (octave_idx_type j = 0 ; j < m2_nc ; j++) \
644 { \
645 octave_quit (); \
646 for (octave_idx_type i = 0 ; i < m2_nr ; i++) \
647 { \
648 if (m2.elem (i, j) == 0.0) \
649 r.elem (i, j) = nan; \
650 else \
651 r.elem (i, j) = m1.elem (0, 0) * m2.elem (i, j); \
652 } \
653 } \
654 r.maybe_compress (true); \
655 } \
656 else \
657 { \
658 r = R (m2); \
659 octave_idx_type m2_nnz = m2.nnz (); \
660 \
661 for (octave_idx_type i = 0 ; i < m2_nnz ; i++) \
662 { \
663 octave_quit (); \
664 r.data (i) = m1.data (0) OP r.data (i); \
665 } \
666 r.maybe_compress (); \
667 } \
668 } \
669 else if (m2_nr == 1 && m2_nc == 1) \
670 { \
671 if (m2.elem (0, 0) == 0.0) \
672 r = R (m1_nr, m1_nc); \
673 else if (octave::math::isnan (m2.elem (0, 0))) \
674 { \
675 r = R (m1_nr, m1_nc, m1_nr * m1_nc); \
676 for (octave_idx_type i = 0 ; i < r.numel () ; i++) \
677 r.elem(i) = nan; \
678 } \
679 else if (octave::math::isinf (m2.elem (0, 0))) \
680 { \
681 r = R (m1_nr, m1_nc, m1_nr * m1_nc); \
682 \
683 for (octave_idx_type j = 0 ; j < m1_nc ; j++) \
684 { \
685 octave_quit (); \
686 for (octave_idx_type i = 0 ; i < m1_nr ; i++) \
687 { \
688 if (m1.elem (i, j) == 0.0) \
689 r.elem (i, j) = nan; \
690 else \
691 r.elem (i, j) = m1.elem (i, j) * m2.elem (0, 0); \
692 } \
693 } \
694 r.maybe_compress (true); \
695 } \
696 else \
697 { \
698 r = R (m1); \
699 octave_idx_type m1_nnz = m1.nnz (); \
700 \
701 for (octave_idx_type i = 0 ; i < m1_nnz ; i++) \
702 { \
703 octave_quit (); \
704 r.data (i) = r.data (i) OP m2.data (0); \
705 } \
706 r.maybe_compress (); \
707 } \
708 } \
709 else if (m1_nr == 0 && (m2_nr == 0 || m2_nr == 1)) \
710 { \
711 if (m1_nc == 1 || m2_nc == 1 || m1_nc == m2_nc) \
712 r.resize (m1_nr, std::max (m1_nc, m2_nc)); \
713 else \
714 octave::err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \
715 } \
716 else if (m1_nc == 0 && (m2_nc == 0 || m2_nc == 1)) \
717 { \
718 if (m1_nr == 1 || m2_nr == 1 || m1_nr == m2_nr) \
719 r.resize (std::max (m1_nr, m2_nr), m1_nc); \
720 else \
721 octave::err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \
722 } \
723 else if (m2_nr == 0 && (m1_nr == 0 || m1_nr == 1)) \
724 { \
725 if (m2_nc == 1 || m1_nc == 1 || m2_nc == m1_nc) \
726 r.resize (m2_nr, std::max (m1_nc, m2_nc)); \
727 else \
728 octave::err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \
729 } \
730 else if (m2_nc == 0 && (m1_nc == 0 || m1_nc == 1)) \
731 { \
732 if (m2_nr == 1 || m1_nr == 1 || m2_nr == m1_nr) \
733 r.resize (std::max (m1_nr, m2_nr), m2_nc); \
734 else \
735 octave::err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \
736 } \
737 else if (m1_nr == m2_nr && m1_nc == m2_nc && \
738 (m1.any_element_is_inf_or_nan () || \
739 m2.any_element_is_inf_or_nan ())) \
740 { \
741 r = R (m1_nr, m1_nc, m1_nr * m1_nc); \
742 \
743 for (octave_idx_type j = 0 ; j < m1_nc ; j++) \
744 { \
745 octave_quit (); \
746 for (octave_idx_type i = 0 ; i < m1_nr ; i++) \
747 r.elem (i, j) = m1.elem (i, j) * m2.elem (i, j); \
748 } \
749 r.maybe_compress (true); \
750 } \
751 else if (m1_nr == m2_nr && m1_nc == m2_nc) \
752 { \
753 r = R (m1_nr, m1_nc, \
754 (m1.nnz () > m2.nnz () ? m1.nnz () : m2.nnz ())); \
755 \
756 octave_idx_type jx = 0; \
757 r.cidx (0) = 0; \
758 for (octave_idx_type i = 0 ; i < m1_nc ; i++) \
759 { \
760 octave_idx_type ja = m1.cidx (i); \
761 octave_idx_type ja_max = m1.cidx (i+1); \
762 bool ja_lt_max = ja < ja_max; \
763 \
764 octave_idx_type jb = m2.cidx (i); \
765 octave_idx_type jb_max = m2.cidx (i+1); \
766 bool jb_lt_max = jb < jb_max; \
767 \
768 while (ja_lt_max || jb_lt_max) \
769 { \
770 octave_quit (); \
771 if ((! jb_lt_max) || \
772 (ja_lt_max && (m1.ridx (ja) < m2.ridx (jb)))) \
773 { \
774 ja++; ja_lt_max= ja < ja_max; \
775 } \
776 else if ((! ja_lt_max) || \
777 (jb_lt_max && (m2.ridx (jb) < m1.ridx (ja)))) \
778 { \
779 jb++; jb_lt_max= jb < jb_max; \
780 } \
781 else \
782 { \
783 if ((m1.data (ja) OP m2.data (jb)) != 0.) \
784 { \
785 r.data (jx) = m1.data (ja) OP m2.data (jb); \
786 r.ridx (jx) = m1.ridx (ja); \
787 jx++; \
788 } \
789 ja++; ja_lt_max= ja < ja_max; \
790 jb++; jb_lt_max= jb < jb_max; \
791 } \
792 } \
793 r.cidx (i+1) = jx; \
794 } \
795 \
796 r.maybe_compress (); \
797 } \
798 else if (m1_nr == m2_nr && (m1_nc == 1 || m2_nc == 1) && \
799 (m1.any_element_is_inf_or_nan () || \
800 m2.any_element_is_inf_or_nan ())) \
801 { \
802 octave_idx_type r_nc = (m1_nc < m2_nc ? m2_nc : m1_nc); \
803 r = R (m1_nr, r_nc, m1_nr * r_nc); \
804 \
805 if (m1_nc == 1) \
806 { \
807 for (octave_idx_type j = 0 ; j < r_nc ; j++) \
808 { \
809 octave_quit (); \
810 for (octave_idx_type i = 0 ; i < m1_nr ; i++) \
811 r.elem (i, j) = m1.elem (i, 0) * m2.elem (i, j); \
812 } \
813 } \
814 else \
815 { \
816 for (octave_idx_type j = 0 ; j < r_nc ; j++) \
817 { \
818 octave_quit (); \
819 for (octave_idx_type i = 0 ; i < m1_nr ; i++) \
820 r.elem (i, j) = m1.elem (i, j) * m2.elem (i, 0); \
821 } \
822 } \
823 r.maybe_compress (true); \
824 } \
825 else if (m1_nr == m2_nr && (m1_nc == 1 || m2_nc == 1)) \
826 { \
827 octave_idx_type r_nc = (m1_nc < m2_nc ? m2_nc : m1_nc); \
828 octave_idx_type rnnz = (m1_nc < m2_nc ? \
829 m1.nnz () * m2_nc + m2.nnz () : \
830 m1.nnz () + m1_nc * m2.nnz ()); \
831 r = R (m1_nr, r_nc, rnnz); \
832 \
833 octave_idx_type jx = 0; \
834 r.cidx (0) = 0; \
835 for (octave_idx_type i = 0 ; i < r_nc ; i++) \
836 { \
837 octave_idx_type ja; \
838 octave_idx_type ja_max; \
839 octave_idx_type jb; \
840 octave_idx_type jb_max; \
841 if (m1_nc == 1) \
842 { \
843 ja = m1.cidx(0); \
844 ja_max = m1.cidx(1); \
845 jb = m2.cidx (i); \
846 jb_max = m2.cidx (i+1); \
847 } \
848 else \
849 { \
850 ja = m1.cidx(i); \
851 ja_max = m1.cidx(i+1); \
852 jb = m2.cidx (0); \
853 jb_max = m2.cidx (1); \
854 } \
855 bool ja_lt_max = ja < ja_max; \
856 bool jb_lt_max = jb < jb_max; \
857 \
858 while (ja_lt_max || jb_lt_max) \
859 { \
860 octave_quit (); \
861 if ((! jb_lt_max) || \
862 (ja_lt_max && (m1.ridx (ja) < m2.ridx (jb)))) \
863 { \
864 ja++; \
865 ja_lt_max = ja < ja_max; \
866 } \
867 else if ((! ja_lt_max) || \
868 (m2.ridx (jb) < m1.ridx (ja))) \
869 { \
870 jb++; \
871 jb_lt_max = jb < jb_max; \
872 } \
873 else \
874 { \
875 if (m1.data (ja) OP m2.data (jb) != 0.) \
876 { \
877 r.data (jx) = m1.data (ja) OP m2.data (jb); \
878 r.ridx (jx) = m1.ridx (ja); \
879 jx++; \
880 } \
881 ja++; \
882 ja_lt_max = ja < ja_max; \
883 jb++; \
884 jb_lt_max = jb < jb_max; \
885 } \
886 } \
887 r.cidx (i+1) = jx; \
888 } \
889 r.maybe_compress (); \
890 } \
891 else if (m1_nc == m2_nc && (m1_nr == 1 || m2_nr == 1)) \
892 r = F (m1.transpose (), m2.transpose ()).transpose (); \
893 else if (m1_nr == 1 && m2_nc == 1 && \
894 (m1.any_element_is_inf_or_nan () || \
895 m2.any_element_is_inf_or_nan ())) \
896 { \
897 r = R (m2_nr, m1_nc, m2_nr * m1_nc); \
898 \
899 for (octave_idx_type j = 0 ; j < m1_nc ; j++) \
900 { \
901 octave_quit (); \
902 for (octave_idx_type i = 0 ; i < m2_nr ; i++) \
903 r.elem (i, j) = m1.elem (0, j) * m2.elem (i, 0); \
904 } \
905 r.maybe_compress (true); \
906 } \
907 else if (m1_nr == 1 && m2_nc == 1) \
908 { \
909 r = R (m2_nr, m1_nc, (m1.nnz () * m2_nr + m2.nnz () * m1_nc)); \
910 \
911 octave_idx_type jx = 0; \
912 r.cidx (0) = 0; \
913 for (octave_idx_type i = 0 ; i < m1_nc ; i++) \
914 { \
915 octave_idx_type ja = m1.cidx (i); \
916 octave_idx_type ja_max = m1.cidx (i+1); \
917 bool ja_lt_max = ja < ja_max; \
918 \
919 octave_idx_type jb = m2.cidx (0); \
920 octave_idx_type jb_max = m2.cidx (1); \
921 bool jb_lt_max = jb < jb_max; \
922 \
923 while (ja_lt_max || jb_lt_max) \
924 { \
925 octave_quit (); \
926 if (! ja_lt_max && jb_lt_max) \
927 { \
928 jb++; \
929 jb_lt_max = jb < jb_max; \
930 } \
931 else if (ja_lt_max && ! jb_lt_max) \
932 { \
933 ja_lt_max = false; \
934 } \
935 else if (ja_lt_max && jb_lt_max) \
936 { \
937 if ((m1.data (ja) OP m2.data (jb)) != 0.) \
938 { \
939 r.data (jx) = m1.data (ja) OP m2.data (jb); \
940 r.ridx (jx) = m2.ridx (jb); \
941 jx++; \
942 } \
943 jb++; \
944 jb_lt_max = jb < jb_max; \
945 ja_lt_max = jb_lt_max; \
946 } \
947 } \
948 r.cidx (i+1) = jx; \
949 } \
950 r.maybe_compress (); \
951 } \
952 else if (m1_nc == 1 && m2_nr == 1) \
953 r = F (m1.transpose (), m2.transpose ()).transpose (); \
954 else \
955 octave::err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \
956 \
957 return r; \
958 }
959
960#define SPARSE_SMSM_BIN_OP_3(R, F, OP, M1, M2) \
961 R \
962 F (const M1& m1, const M2& m2) \
963 { \
964 R r; \
965 Complex Zero = Complex (); \
966 \
967 octave_idx_type m1_nr = m1.rows (); \
968 octave_idx_type m1_nc = m1.cols (); \
969 \
970 octave_idx_type m2_nr = m2.rows (); \
971 octave_idx_type m2_nc = m2.cols (); \
972 \
973 if (m1_nr == 1 && m1_nc == 1) \
974 { \
975 if ((m1.elem (0,0) OP Zero) == Zero) \
976 { \
977 octave_idx_type m2_nnz = m2.nnz (); \
978 r = R (m2); \
979 for (octave_idx_type i = 0 ; i < m2_nnz ; i++) \
980 r.data (i) = m1.elem (0,0) OP r.data (i); \
981 r.maybe_compress (); \
982 } \
983 else \
984 { \
985 r = R (m2_nr, m2_nc, m1.elem (0,0) OP Zero); \
986 for (octave_idx_type j = 0 ; j < m2_nc ; j++) \
987 { \
988 octave_quit (); \
989 octave_idx_type idxj = j * m2_nr; \
990 for (octave_idx_type i = m2.cidx (j) ; i < m2.cidx (j+1) ; i++) \
991 { \
992 octave_quit (); \
993 r.data (idxj + m2.ridx (i)) = m1.elem (0,0) OP m2.data (i); \
994 } \
995 } \
996 r.maybe_compress (); \
997 } \
998 } \
999 else if (m2_nr == 1 && m2_nc == 1) \
1000 { \
1001 if ((Zero OP m1.elem (0,0)) == Zero) \
1002 { \
1003 octave_idx_type m1_nnz = m1.nnz (); \
1004 r = R (m1); \
1005 for (octave_idx_type i = 0 ; i < m1_nnz ; i++) \
1006 r.data (i) = r.data (i) OP m2.elem (0,0); \
1007 r.maybe_compress (); \
1008 } \
1009 else \
1010 { \
1011 r = R (m1_nr, m1_nc, Zero OP m2.elem (0,0)); \
1012 for (octave_idx_type j = 0 ; j < m1_nc ; j++) \
1013 { \
1014 octave_quit (); \
1015 octave_idx_type idxj = j * m1_nr; \
1016 for (octave_idx_type i = m1.cidx (j) ; i < m1.cidx (j+1) ; i++) \
1017 { \
1018 octave_quit (); \
1019 r.data (idxj + m1.ridx (i)) = m1.data (i) OP m2.elem (0,0); \
1020 } \
1021 } \
1022 r.maybe_compress (); \
1023 } \
1024 } \
1025 else if (m1_nr == 0 && (m2_nr == 0 || m2_nr == 1)) \
1026 { \
1027 if (m1_nc == 1 || m2_nc == 1 || m1_nc == m2_nc) \
1028 r.resize (m1_nr, std::max (m1_nc, m2_nc)); \
1029 else \
1030 octave::err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \
1031 } \
1032 else if (m1_nc == 0 && (m2_nc == 0 || m2_nc == 1)) \
1033 { \
1034 if (m1_nr == 1 || m2_nr == 1 || m1_nr == m2_nr) \
1035 r.resize (std::max (m1_nr, m2_nr), m1_nc); \
1036 else \
1037 octave::err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \
1038 } \
1039 else if (m2_nr == 0 && (m1_nr == 0 || m1_nr == 1)) \
1040 { \
1041 if (m2_nc == 1 || m1_nc == 1 || m2_nc == m1_nc) \
1042 r.resize (m2_nr, std::max (m1_nc, m2_nc)); \
1043 else \
1044 octave::err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \
1045 } \
1046 else if (m2_nc == 0 && (m1_nc == 0 || m1_nc == 1)) \
1047 { \
1048 if (m2_nr == 1 || m1_nr == 1 || m2_nr == m1_nr) \
1049 r.resize (std::max (m1_nr, m2_nr), m2_nc); \
1050 else \
1051 octave::err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \
1052 } \
1053 else if (m1_nr == m2_nr && m1_nc == m2_nc) \
1054 { \
1055 \
1056 /* FIXME: Kludge... Always double/Complex, so Complex () */ \
1057 r = R (m1_nr, m1_nc, (Zero OP Zero)); \
1058 \
1059 for (octave_idx_type i = 0 ; i < m1_nc ; i++) \
1060 { \
1061 octave_idx_type ja = m1.cidx (i); \
1062 octave_idx_type ja_max = m1.cidx (i+1); \
1063 bool ja_lt_max = ja < ja_max; \
1064 \
1065 octave_idx_type jb = m2.cidx (i); \
1066 octave_idx_type jb_max = m2.cidx (i+1); \
1067 bool jb_lt_max = jb < jb_max; \
1068 \
1069 while (ja_lt_max || jb_lt_max) \
1070 { \
1071 octave_quit (); \
1072 if ((! jb_lt_max) || \
1073 (ja_lt_max && (m1.ridx (ja) < m2.ridx (jb)))) \
1074 { \
1075 /* keep those kludges coming */ \
1076 r.elem (m1.ridx (ja),i) = m1.data (ja) OP Zero; \
1077 ja++; \
1078 ja_lt_max= ja < ja_max; \
1079 } \
1080 else if ((! ja_lt_max) || \
1081 (jb_lt_max && (m2.ridx (jb) < m1.ridx (ja)))) \
1082 { \
1083 /* keep those kludges coming */ \
1084 r.elem (m2.ridx (jb),i) = Zero OP m2.data (jb); \
1085 jb++; \
1086 jb_lt_max= jb < jb_max; \
1087 } \
1088 else \
1089 { \
1090 r.elem (m1.ridx (ja),i) = m1.data (ja) OP \
1091 m2.data (jb); \
1092 ja++; \
1093 ja_lt_max= ja < ja_max; \
1094 jb++; \
1095 jb_lt_max= jb < jb_max; \
1096 } \
1097 } \
1098 } \
1099 r.maybe_compress (true); \
1100 } \
1101 else if (m1_nr == m2_nr && (m1_nc == 1 || m2_nc == 1)) \
1102 { \
1103 octave_idx_type r_nc = (m1_nc < m2_nc ? m2_nc : m1_nc); \
1104 r = R (m1_nr, r_nc, (Zero / Zero)); \
1105 \
1106 for (octave_idx_type i = 0 ; i < r_nc ; i++) \
1107 { \
1108 octave_idx_type ja; \
1109 octave_idx_type ja_max; \
1110 octave_idx_type jb; \
1111 octave_idx_type jb_max; \
1112 if (m1_nc == 1) \
1113 { \
1114 ja = m1.cidx(0); \
1115 ja_max = m1.cidx(1); \
1116 jb = m2.cidx (i); \
1117 jb_max = m2.cidx (i+1); \
1118 } \
1119 else \
1120 { \
1121 ja = m1.cidx(i); \
1122 ja_max = m1.cidx(i+1); \
1123 jb = m2.cidx (0); \
1124 jb_max = m2.cidx (1); \
1125 } \
1126 bool ja_lt_max = ja < ja_max; \
1127 bool jb_lt_max = jb < jb_max; \
1128 \
1129 while (ja_lt_max || jb_lt_max) \
1130 { \
1131 octave_quit (); \
1132 if ((! jb_lt_max) || \
1133 (ja_lt_max && (m1.ridx (ja) < m2.ridx (jb)))) \
1134 { \
1135 r.elem (m1.ridx (ja), i) = m1.data (ja) OP Zero; \
1136 ja++; \
1137 ja_lt_max = ja < ja_max; \
1138 } \
1139 else if ((! ja_lt_max) || \
1140 (m2.ridx (jb) < m1.ridx (ja))) \
1141 { \
1142 r.elem (m2.ridx (jb), i) = Zero OP m2.data (jb); \
1143 jb++; \
1144 jb_lt_max = jb < jb_max; \
1145 } \
1146 else \
1147 { \
1148 r.elem (m1.ridx (ja), i) = m1.data (ja) OP \
1149 m2.data (jb); \
1150 ja++; \
1151 ja_lt_max = ja < ja_max; \
1152 jb++; \
1153 jb_lt_max = jb < jb_max; \
1154 } \
1155 } \
1156 } \
1157 r.maybe_compress (true); \
1158 } \
1159 else if (m1_nc == m2_nc && (m1_nr == 1 || m2_nr == 1)) \
1160 r = quotient (m1.transpose (), m2.transpose ()).transpose (); \
1161 else if (m1_nr == 1 && m2_nc == 1) \
1162 { \
1163 r = R (m2_nr, m1_nc, (m1_nc * m2_nr)); \
1164 \
1165 for (octave_idx_type j = 0; j < m1_nc; j++) \
1166 { \
1167 const Complex m1_val = m1.elem (0, j); \
1168 \
1169 for (octave_idx_type i = 0; i < m2_nr; i++) \
1170 { \
1171 octave_quit (); \
1172 \
1173 const Complex m2_val = m2.elem (i, 0); \
1174 const Complex val = m1_val OP m2_val; \
1175 \
1176 if (val != Zero || val != val) \
1177 r.elem (i, j) = val; \
1178 } \
1179 } \
1180 r.maybe_compress (true); \
1181 } \
1182 else if (m1_nc == 1 && m2_nr == 1) \
1183 r = quotient (m1.transpose (), m2.transpose ()).transpose (); \
1184 else \
1185 octave::err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \
1186 \
1187 return r; \
1188 }
1189
1190// Note that SM ./ SM needs to take into account the NaN and Inf values
1191// implied by the division by zero.
1192// FIXME: Are the NaNs double(NaN) or Complex(NaN,Nan) in the complex case?
1193#define SPARSE_SMSM_BIN_OPS(R1, R2, M1, M2) \
1194 SPARSE_SMSM_BIN_OP_1 (R1, operator +, +, M1, M2) \
1195 SPARSE_SMSM_BIN_OP_1 (R1, operator -, -, M1, M2) \
1196 SPARSE_SMSM_BIN_OP_2 (R2, product, *, M1, M2) \
1197 SPARSE_SMSM_BIN_OP_3 (R2, quotient, /, M1, M2)
1198
1199// FIXME: this macro duplicates the bodies of the template functions
1200// defined in the SPARSE_SSM_CMP_OP and SPARSE_SMS_CMP_OP macros.
1201
1202#define SPARSE_SMSM_CMP_OP(F, OP, M1, M2) \
1203 SparseBoolMatrix \
1204 F (const M1& m1, const M2& m2) \
1205 { \
1206 SparseBoolMatrix r; \
1207 \
1208 octave_idx_type m1_nr = m1.rows (); \
1209 octave_idx_type m1_nc = m1.cols (); \
1210 \
1211 octave_idx_type m2_nr = m2.rows (); \
1212 octave_idx_type m2_nc = m2.cols (); \
1213 \
1214 M1::element_type Z1 = M1::element_type (); \
1215 M2::element_type Z2 = M2::element_type (); \
1216 \
1217 if (m1_nr == 1 && m1_nc == 1) \
1218 { \
1219 if (m1.elem (0,0) OP Z2) \
1220 { \
1221 r = SparseBoolMatrix (m2_nr, m2_nc, true); \
1222 for (octave_idx_type j = 0; j < m2_nc; j++) \
1223 for (octave_idx_type i = m2.cidx (j); i < m2.cidx (j+1); i++) \
1224 if (! (m1.elem (0,0) OP m2.data (i))) \
1225 r.data (m2.ridx (i) + j * m2_nr) = false; \
1226 r.maybe_compress (true); \
1227 } \
1228 else \
1229 { \
1230 r = SparseBoolMatrix (m2_nr, m2_nc, m2.nnz ()); \
1231 r.cidx (0) = static_cast<octave_idx_type> (0); \
1232 octave_idx_type nel = 0; \
1233 for (octave_idx_type j = 0; j < m2_nc; j++) \
1234 { \
1235 for (octave_idx_type i = m2.cidx (j); i < m2.cidx (j+1); i++) \
1236 if (m1.elem (0,0) OP m2.data (i)) \
1237 { \
1238 r.ridx (nel) = m2.ridx (i); \
1239 r.data (nel++) = true; \
1240 } \
1241 r.cidx (j + 1) = nel; \
1242 } \
1243 r.maybe_compress (false); \
1244 } \
1245 } \
1246 else if (m2_nr == 1 && m2_nc == 1) \
1247 { \
1248 if (Z1 OP m2.elem (0,0)) \
1249 { \
1250 r = SparseBoolMatrix (m1_nr, m1_nc, true); \
1251 for (octave_idx_type j = 0; j < m1_nc; j++) \
1252 for (octave_idx_type i = m1.cidx (j); i < m1.cidx (j+1); i++) \
1253 if (! (m1.data (i) OP m2.elem (0,0))) \
1254 r.data (m1.ridx (i) + j * m1_nr) = false; \
1255 r.maybe_compress (true); \
1256 } \
1257 else \
1258 { \
1259 r = SparseBoolMatrix (m1_nr, m1_nc, m1.nnz ()); \
1260 r.cidx (0) = static_cast<octave_idx_type> (0); \
1261 octave_idx_type nel = 0; \
1262 for (octave_idx_type j = 0; j < m1_nc; j++) \
1263 { \
1264 for (octave_idx_type i = m1.cidx (j); i < m1.cidx (j+1); i++) \
1265 if (m1.data (i) OP m2.elem (0,0)) \
1266 { \
1267 r.ridx (nel) = m1.ridx (i); \
1268 r.data (nel++) = true; \
1269 } \
1270 r.cidx (j + 1) = nel; \
1271 } \
1272 r.maybe_compress (false); \
1273 } \
1274 } \
1275 else if (m1_nr == 0 && (m2_nr == 0 || m2_nr == 1)) \
1276 { \
1277 if (m1_nc == 1 || m2_nc == 1 || m1_nc == m2_nc) \
1278 r.resize (m1_nr, std::max (m1_nc, m2_nc)); \
1279 else \
1280 octave::err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \
1281 } \
1282 else if (m1_nc == 0 && (m2_nc == 0 || m2_nc == 1)) \
1283 { \
1284 if (m1_nr == 1 || m2_nr == 1 || m1_nr == m2_nr) \
1285 r.resize (std::max (m1_nr, m2_nr), m1_nc); \
1286 else \
1287 octave::err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \
1288 } \
1289 else if (m2_nr == 0 && (m1_nr == 0 || m1_nr == 1)) \
1290 { \
1291 if (m2_nc == 1 || m1_nc == 1 || m2_nc == m1_nc) \
1292 r.resize (m2_nr, std::max (m1_nc, m2_nc)); \
1293 else \
1294 octave::err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \
1295 } \
1296 else if (m2_nc == 0 && (m1_nc == 0 || m1_nc == 1)) \
1297 { \
1298 if (m2_nr == 1 || m1_nr == 1 || m2_nr == m1_nr) \
1299 r.resize (std::max (m1_nr, m2_nr), m2_nc); \
1300 else \
1301 octave::err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \
1302 } \
1303 else if (m1_nr == m2_nr && m1_nc == m2_nc) \
1304 { \
1305 if (Z1 OP Z2) \
1306 { \
1307 r = SparseBoolMatrix (m1_nr, m1_nc, true); \
1308 for (octave_idx_type j = 0; j < m1_nc; j++) \
1309 { \
1310 octave_idx_type i1 = m1.cidx (j); \
1311 octave_idx_type e1 = m1.cidx (j+1); \
1312 octave_idx_type i2 = m2.cidx (j); \
1313 octave_idx_type e2 = m2.cidx (j+1); \
1314 while (i1 < e1 || i2 < e2) \
1315 { \
1316 octave_quit (); \
1317 \
1318 if (i1 == e1 || \
1319 (i2 < e2 && m1.ridx (i1) > m2.ridx (i2))) \
1320 { \
1321 if (! (Z1 OP m2.data (i2))) \
1322 r.data (m2.ridx (i2) + j * m1_nr) = false; \
1323 i2++; \
1324 } \
1325 else if (i2 == e2 || m1.ridx (i1) < m2.ridx (i2)) \
1326 { \
1327 if (! (m1.data (i1) OP Z2)) \
1328 r.data (m1.ridx (i1) + j * m1_nr) = false; \
1329 i1++; \
1330 } \
1331 else \
1332 { \
1333 if (! (m1.data (i1) OP m2.data (i2))) \
1334 r.data (m1.ridx (i1) + j * m1_nr) = false; \
1335 i1++; \
1336 i2++; \
1337 } \
1338 } \
1339 } \
1340 r.maybe_compress (true); \
1341 } \
1342 else \
1343 { \
1344 r = SparseBoolMatrix (m1_nr, m1_nc, m1.nnz () + m2.nnz ()); \
1345 r.cidx (0) = static_cast<octave_idx_type> (0); \
1346 octave_idx_type nel = 0; \
1347 for (octave_idx_type j = 0; j < m1_nc; j++) \
1348 { \
1349 octave_idx_type i1 = m1.cidx (j); \
1350 octave_idx_type e1 = m1.cidx (j+1); \
1351 octave_idx_type i2 = m2.cidx (j); \
1352 octave_idx_type e2 = m2.cidx (j+1); \
1353 while (i1 < e1 || i2 < e2) \
1354 { \
1355 octave_quit (); \
1356 \
1357 if (i1 == e1 || \
1358 (i2 < e2 && m1.ridx (i1) > m2.ridx (i2))) \
1359 { \
1360 if (Z1 OP m2.data (i2)) \
1361 { \
1362 r.ridx (nel) = m2.ridx (i2); \
1363 r.data (nel++) = true; \
1364 } \
1365 i2++; \
1366 } \
1367 else if (i2 == e2 || m1.ridx (i1) < m2.ridx (i2)) \
1368 { \
1369 if (m1.data (i1) OP Z2) \
1370 { \
1371 r.ridx (nel) = m1.ridx (i1); \
1372 r.data (nel++) = true; \
1373 } \
1374 i1++; \
1375 } \
1376 else \
1377 { \
1378 if (m1.data (i1) OP m2.data (i2)) \
1379 { \
1380 r.ridx (nel) = m1.ridx (i1); \
1381 r.data (nel++) = true; \
1382 } \
1383 i1++; \
1384 i2++; \
1385 } \
1386 } \
1387 r.cidx (j + 1) = nel; \
1388 } \
1389 r.maybe_compress (false); \
1390 } \
1391 } \
1392 else if (m1_nr == m2_nr && (m1_nc == 1 || m2_nc == 1)) \
1393 { \
1394 if (Z1 OP Z2) \
1395 { \
1396 octave_idx_type r_nc = (m1_nc < m2_nc ? m2_nc : m1_nc); \
1397 r = SparseBoolMatrix (m1_nr, r_nc, true); \
1398 for (octave_idx_type j = 0; j < r_nc; j++) \
1399 { \
1400 octave_idx_type i1; \
1401 octave_idx_type e1; \
1402 octave_idx_type i2; \
1403 octave_idx_type e2; \
1404 if (m1_nc == 1) \
1405 { \
1406 i1 = m1.cidx(0); \
1407 e1 = m1.cidx(1); \
1408 i2 = m2.cidx (j); \
1409 e2 = m2.cidx (j+1); \
1410 } \
1411 else \
1412 { \
1413 i1 = m1.cidx(j); \
1414 e1 = m1.cidx(j+1); \
1415 i2 = m2.cidx (0); \
1416 e2 = m2.cidx (1); \
1417 } \
1418 while (i1 < e1 || i2 < e2) \
1419 { \
1420 octave_quit (); \
1421 \
1422 if (i1 == e1 || \
1423 (i2 < e2 && m1.ridx (i1) > m2.ridx (i2))) \
1424 { \
1425 if (! (Z1 OP m2.data (i2))) \
1426 r.data (m2.ridx (i2) + j * m1_nr) = false; \
1427 i2++; \
1428 } \
1429 else if (i2 == e2 || m1.ridx (i1) < m2.ridx (i2)) \
1430 { \
1431 if (! (m1.data (i1) OP Z2)) \
1432 r.data (m1.ridx (i1) + j * m1_nr) = false; \
1433 i1++; \
1434 } \
1435 else \
1436 { \
1437 if (! (m1.data (i1) OP m2.data (i2))) \
1438 r.data (m1.ridx (i1) + j * m1_nr) = false; \
1439 i1++; \
1440 i2++; \
1441 } \
1442 } \
1443 } \
1444 r.maybe_compress (true); \
1445 } \
1446 else \
1447 { \
1448 octave_idx_type r_nc = (m1_nc < m2_nc ? m2_nc : m1_nc); \
1449 octave_idx_type rnnz = (m1_nc < m2_nc ? \
1450 m1.nnz () * m2_nc + m2.nnz () : \
1451 m1.nnz () + m1_nc * m2.nnz ()); \
1452 r = SparseBoolMatrix (m1_nr, r_nc, rnnz); \
1453 \
1454 r.cidx (0) = static_cast<octave_idx_type> (0); \
1455 octave_idx_type nel = 0; \
1456 for (octave_idx_type j = 0; j < r_nc; j++) \
1457 { \
1458 octave_idx_type i1; \
1459 octave_idx_type e1; \
1460 octave_idx_type i2; \
1461 octave_idx_type e2; \
1462 if (m1_nc == 1) \
1463 { \
1464 i1 = m1.cidx(0); \
1465 e1 = m1.cidx(1); \
1466 i2 = m2.cidx (j); \
1467 e2 = m2.cidx (j+1); \
1468 } \
1469 else \
1470 { \
1471 i1 = m1.cidx(j); \
1472 e1 = m1.cidx(j+1); \
1473 i2 = m2.cidx (0); \
1474 e2 = m2.cidx (1); \
1475 } \
1476 while (i1 < e1 || i2 < e2) \
1477 { \
1478 octave_quit (); \
1479 \
1480 if (i1 == e1 || \
1481 (i2 < e2 && m1.ridx (i1) > m2.ridx (i2))) \
1482 { \
1483 if (Z1 OP m2.data (i2)) \
1484 { \
1485 r.ridx (nel) = m2.ridx (i2); \
1486 r.data (nel++) = true; \
1487 } \
1488 i2++; \
1489 } \
1490 else if (i2 == e2 || m1.ridx (i1) < m2.ridx (i2)) \
1491 { \
1492 if (m1.data (i1) OP Z2) \
1493 { \
1494 r.ridx (nel) = m1.ridx (i1); \
1495 r.data (nel++) = true; \
1496 } \
1497 i1++; \
1498 } \
1499 else \
1500 { \
1501 if (m1.data (i1) OP m2.data (i2)) \
1502 { \
1503 r.ridx (nel) = m1.ridx (i1); \
1504 r.data (nel++) = true; \
1505 } \
1506 i1++; \
1507 i2++; \
1508 } \
1509 } \
1510 r.cidx (j + 1) = nel; \
1511 } \
1512 r.maybe_compress (false); \
1513 } \
1514 } \
1515 else if (m1_nc == m2_nc && (m1_nr == 1 || m2_nr == 1)) \
1516 r = F (m1.transpose (), m2.transpose ()).transpose (); \
1517 else if (m1_nr == 1 && m2_nc == 1) \
1518 { \
1519 if (Z1 OP Z2) \
1520 { \
1521 r = SparseBoolMatrix (m2_nr, m1_nc, true); \
1522 \
1523 for (octave_idx_type j = 0; j < m1_nc; j++) \
1524 { \
1525 const Complex m1_val = m1.elem (0, j); \
1526 \
1527 for (octave_idx_type i = 0; i < m2_nr; i++) \
1528 { \
1529 octave_quit (); \
1530 \
1531 const Complex m2_val = m2.elem (i, 0); \
1532 \
1533 if (! (m1_val OP m2_val)) \
1534 r.elem (i, j) = false; \
1535 } \
1536 } \
1537 r.maybe_compress (true); \
1538 } \
1539 else \
1540 { \
1541 octave_idx_type nz = m1.nnz () * m2_nr + m2.nnz () * m1_nc; \
1542 r = SparseBoolMatrix (m2_nr, m1_nc, nz); \
1543 \
1544 for (octave_idx_type j = 0; j < m1_nc; j++) \
1545 { \
1546 const Complex m1_val = m1.elem (0, j); \
1547 \
1548 for (octave_idx_type i = 0; i < m2_nr; i++) \
1549 { \
1550 octave_quit (); \
1551 \
1552 const Complex m2_val = m2.elem (i, 0); \
1553 \
1554 if (m1_val OP m2_val) \
1555 r.elem (i, j) = true; \
1556 } \
1557 } \
1558 r.maybe_compress (false); \
1559 } \
1560 } \
1561 else if (m1_nc == 1 && m2_nr == 1) \
1562 r = F (m1.transpose (), m2.transpose ()).transpose (); \
1563 else \
1564 octave::err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \
1565 \
1566 return r; \
1567 }
1568
1569#define SPARSE_SMSM_CMP_OPS(M1, M2) \
1570 SPARSE_SMSM_CMP_OP (mx_el_lt, <, M1, M2) \
1571 SPARSE_SMSM_CMP_OP (mx_el_le, <=, M1, M2) \
1572 SPARSE_SMSM_CMP_OP (mx_el_ge, >=, M1, M2) \
1573 SPARSE_SMSM_CMP_OP (mx_el_gt, >, M1, M2) \
1574 SPARSE_SMSM_CMP_OP (mx_el_eq, ==, M1, M2) \
1575 SPARSE_SMSM_CMP_OP (mx_el_ne, !=, M1, M2)
1576
1577#define SPARSE_SMSM_EQNE_OPS(M1, M2) \
1578 SPARSE_SMSM_CMP_OP (mx_el_eq, ==, M1, M2) \
1579 SPARSE_SMSM_CMP_OP (mx_el_ne, !=, M1, M2)
1580
1581#define SPARSE_SMSM_BOOL_AND_OP(M1, M2) \
1582 extern SparseBoolMatrix mx_el_and (const M1&, const M2::element_type&); \
1583 extern SparseBoolMatrix mx_el_and (const M1::element_type&, const M2&); \
1584 SparseBoolMatrix \
1585 mx_el_and (const M1& m1, const M2& m2) \
1586 { \
1587 SparseBoolMatrix r; \
1588 \
1589 octave_idx_type m1_nr = m1.rows (); \
1590 octave_idx_type m1_nc = m1.cols (); \
1591 \
1592 octave_idx_type m2_nr = m2.rows (); \
1593 octave_idx_type m2_nc = m2.cols (); \
1594 \
1595 M1::element_type lhs_zero = M1::element_type (); \
1596 M2::element_type rhs_zero = M2::element_type (); \
1597 \
1598 if (m1_nr == 1 && m1_nc == 1) \
1599 return mx_el_and (m1.elem (0,0), m2); \
1600 else if (m2_nr == 1 && m2_nc == 1) \
1601 return mx_el_and (m1, m2.elem (0,0)); \
1602 else if (m1_nr == 0 && (m2_nr == 0 || m2_nr == 1)) \
1603 { \
1604 if (m1_nc == 1 || m2_nc == 1 || m1_nc == m2_nc) \
1605 r.resize (m1_nr, std::max (m1_nc, m2_nc)); \
1606 else \
1607 octave::err_nonconformant ("mx_el_and", m1_nr, m1_nc, m2_nr, m2_nc); \
1608 } \
1609 else if (m1_nc == 0 && (m2_nc == 0 || m2_nc == 1)) \
1610 { \
1611 if (m1_nr == 1 || m2_nr == 1 || m1_nr == m2_nr) \
1612 r.resize (std::max (m1_nr, m2_nr), m1_nc); \
1613 else \
1614 octave::err_nonconformant ("mx_el_and", m1_nr, m1_nc, m2_nr, m2_nc); \
1615 } \
1616 else if (m2_nr == 0 && (m1_nr == 0 || m1_nr == 1)) \
1617 { \
1618 if (m2_nc == 1 || m1_nc == 1 || m2_nc == m1_nc) \
1619 r.resize (m2_nr, std::max (m1_nc, m2_nc)); \
1620 else \
1621 octave::err_nonconformant ("mx_el_and", m1_nr, m1_nc, m2_nr, m2_nc); \
1622 } \
1623 else if (m2_nc == 0 && (m1_nc == 0 || m1_nc == 1)) \
1624 { \
1625 if (m2_nr == 1 || m1_nr == 1 || m2_nr == m1_nr) \
1626 r.resize (std::max (m1_nr, m2_nr), m2_nc); \
1627 else \
1628 octave::err_nonconformant ("mx_el_and", m1_nr, m1_nc, m2_nr, m2_nc); \
1629 } \
1630 else if (m1_nr == m2_nr && m1_nc == m2_nc) \
1631 { \
1632 r = SparseBoolMatrix (m1_nr, m1_nc, m1.nnz () + m2.nnz ()); \
1633 r.cidx (0) = static_cast<octave_idx_type> (0); \
1634 octave_idx_type nel = 0; \
1635 for (octave_idx_type j = 0; j < m1_nc; j++) \
1636 { \
1637 octave_idx_type i1 = m1.cidx (j); \
1638 octave_idx_type e1 = m1.cidx (j+1); \
1639 octave_idx_type i2 = m2.cidx (j); \
1640 octave_idx_type e2 = m2.cidx (j+1); \
1641 while (i1 < e1 || i2 < e2) \
1642 { \
1643 octave_quit (); \
1644 \
1645 if (i1 == e1 || \
1646 (i2 < e2 && m1.ridx (i1) > m2.ridx (i2))) \
1647 i2++; \
1648 else if (i2 == e2 || m1.ridx (i1) < m2.ridx (i2)) \
1649 i1++; \
1650 else \
1651 { \
1652 if ((m1.data (i1) != lhs_zero) && \
1653 (m2.data (i2) != rhs_zero)) \
1654 { \
1655 r.ridx (nel) = m1.ridx (i1); \
1656 r.data (nel++) = true; \
1657 } \
1658 i1++; \
1659 i2++; \
1660 } \
1661 } \
1662 r.cidx (j + 1) = nel; \
1663 } \
1664 r.maybe_compress (false); \
1665 } \
1666 else if (m1_nr == m2_nr && (m1_nc == 1 || m2_nc == 1)) \
1667 { \
1668 octave_idx_type r_nc = (m1_nc < m2_nc ? m2_nc : m1_nc); \
1669 octave_idx_type rnnz = (m1_nc < m2_nc ? \
1670 m1.nnz () * m2_nc + m2.nnz () : \
1671 m1.nnz () + m1_nc * m2.nnz ()); \
1672 r = SparseBoolMatrix (m1_nr, r_nc, rnnz); \
1673 \
1674 r.cidx (0) = static_cast<octave_idx_type> (0); \
1675 octave_idx_type nel = 0; \
1676 for (octave_idx_type j = 0; j < r_nc; j++) \
1677 { \
1678 octave_idx_type i1; \
1679 octave_idx_type e1; \
1680 octave_idx_type i2; \
1681 octave_idx_type e2; \
1682 if (m1_nc == 1) \
1683 { \
1684 i1 = m1.cidx(0); \
1685 e1 = m1.cidx(1); \
1686 i2 = m2.cidx (j); \
1687 e2 = m2.cidx (j+1); \
1688 } \
1689 else \
1690 { \
1691 i1 = m1.cidx(j); \
1692 e1 = m1.cidx(j+1); \
1693 i2 = m2.cidx (0); \
1694 e2 = m2.cidx (1); \
1695 } \
1696 while (i1 < e1 || i2 < e2) \
1697 { \
1698 octave_quit (); \
1699 \
1700 if (i1 == e1 || \
1701 (i2 < e2 && m1.ridx (i1) > m2.ridx (i2))) \
1702 i2++; \
1703 else if (i2 == e2 || m1.ridx (i1) < m2.ridx (i2)) \
1704 i1++; \
1705 else \
1706 { \
1707 if ((m1.data (i1) != lhs_zero) && \
1708 (m2.data (i2) != rhs_zero)) \
1709 { \
1710 r.ridx (nel) = m1.ridx (i1); \
1711 r.data (nel++) = true; \
1712 } \
1713 i1++; \
1714 i2++; \
1715 } \
1716 } \
1717 r.cidx (j + 1) = nel; \
1718 } \
1719 r.maybe_compress (false); \
1720 } \
1721 else if (m1_nc == m2_nc && (m1_nr == 1 || m2_nr == 1)) \
1722 r = mx_el_and (m1.transpose (), m2.transpose ()).transpose (); \
1723 else if (m1_nr == 1 && m2_nc == 1) \
1724 { \
1725 /* Count num of nonzero elements */ \
1726 octave_idx_type nel = 0; \
1727 for (octave_idx_type j = 0; j < m1_nc; j++) \
1728 for (octave_idx_type i = 0; i < m2_nr; i++) \
1729 if ((m1.elem (0, j) != lhs_zero) && \
1730 (m2.elem (i, 0) != rhs_zero)) \
1731 nel++; \
1732 \
1733 r = SparseBoolMatrix (m2_nr, m1_nc, nel); \
1734 \
1735 for (octave_idx_type j = 0; j < m1_nc; j++) \
1736 { \
1737 const bool m1_val = m1.elem (0, j) != lhs_zero; \
1738 \
1739 for (octave_idx_type i = 0; i < m2_nr; i++) \
1740 { \
1741 octave_quit (); \
1742 \
1743 const bool m2_val = m2.elem (i, 0) != rhs_zero; \
1744 \
1745 if (m1_val && m2_val) \
1746 r.elem (i, j) = true; \
1747 } \
1748 } \
1749 r.maybe_compress (true); \
1750 } \
1751 else if (m1_nc == 1 && m2_nr == 1) \
1752 r = mx_el_and (m1.transpose (), m2.transpose ()).transpose (); \
1753 else \
1754 octave::err_nonconformant ("mx_el_and", m1_nr, m1_nc, m2_nr, m2_nc); \
1755 \
1756 return r; \
1757 }
1758
1759#define SPARSE_SMSM_BOOL_OR_OP(M1, M2) \
1760 extern SparseBoolMatrix mx_el_or (const M1&, const M2&); \
1761 SparseBoolMatrix \
1762 mx_el_or (const M1& m1, const M2& m2) \
1763 { \
1764 SparseBoolMatrix r; \
1765 \
1766 octave_idx_type m1_nr = m1.rows (); \
1767 octave_idx_type m1_nc = m1.cols (); \
1768 \
1769 octave_idx_type m2_nr = m2.rows (); \
1770 octave_idx_type m2_nc = m2.cols (); \
1771 \
1772 M1::element_type lhs_zero = M1::element_type (); \
1773 M2::element_type rhs_zero = M2::element_type (); \
1774 \
1775 if (m1_nr == 1 && m1_nc == 1) \
1776 { \
1777 if (m2_nr > 0 && m2_nc > 0) \
1778 { \
1779 if (m1.elem(0, 0) != rhs_zero) \
1780 r = SparseBoolMatrix (m2_nr, m2_nc, true); \
1781 else \
1782 { \
1783 r = SparseBoolMatrix (m2_nr, m2_nc, m2.nnz ()); \
1784 \
1785 for (octave_idx_type j = 0; j < m2_nc; j++) \
1786 for (octave_idx_type i = 0; i < m2_nr; i++) \
1787 if (m2.elem (i, j) != lhs_zero) \
1788 r.elem (i, j) = true; \
1789 } \
1790 } \
1791 else \
1792 r = SparseBoolMatrix (m2_nr, m2_nc); \
1793 return r; \
1794 } \
1795 else if (m2_nr == 1 && m2_nc == 1) \
1796 { \
1797 if (m1_nr > 0 && m1_nc > 0) \
1798 { \
1799 if (m2.elem(0, 0) != rhs_zero) \
1800 r = SparseBoolMatrix (m1_nr, m1_nc, true); \
1801 else \
1802 { \
1803 r = SparseBoolMatrix (m1_nr, m1_nc, m1.nnz ()); \
1804 \
1805 for (octave_idx_type j = 0; j < m1_nc; j++) \
1806 for (octave_idx_type i = 0; i < m1_nr; i++) \
1807 if (m1.elem (i, j) != lhs_zero) \
1808 r.elem (i, j) = true; \
1809 } \
1810 } \
1811 else \
1812 r = SparseBoolMatrix (m1_nr, m1_nc); \
1813 return r; \
1814 } \
1815 else if (m1_nr == 0 && (m2_nr == 0 || m2_nr == 1)) \
1816 { \
1817 if (m1_nc == 1 || m2_nc == 1 || m1_nc == m2_nc) \
1818 r.resize (m1_nr, std::max (m1_nc, m2_nc)); \
1819 else \
1820 octave::err_nonconformant ("mx_el_or", m1_nr, m1_nc, m2_nr, m2_nc); \
1821 } \
1822 else if (m1_nc == 0 && (m2_nc == 0 || m2_nc == 1)) \
1823 { \
1824 if (m1_nr == 1 || m2_nr == 1 || m1_nr == m2_nr) \
1825 r.resize (std::max (m1_nr, m2_nr), m1_nc); \
1826 else \
1827 octave::err_nonconformant ("mx_el_or", m1_nr, m1_nc, m2_nr, m2_nc); \
1828 } \
1829 else if (m2_nr == 0 && (m1_nr == 0 || m1_nr == 1)) \
1830 { \
1831 if (m2_nc == 1 || m1_nc == 1 || m2_nc == m1_nc) \
1832 r.resize (m2_nr, std::max (m1_nc, m2_nc)); \
1833 else \
1834 octave::err_nonconformant ("mx_el_or", m1_nr, m1_nc, m2_nr, m2_nc); \
1835 } \
1836 else if (m2_nc == 0 && (m1_nc == 0 || m1_nc == 1)) \
1837 { \
1838 if (m2_nr == 1 || m1_nr == 1 || m2_nr == m1_nr) \
1839 r.resize (std::max (m1_nr, m2_nr), m2_nc); \
1840 else \
1841 octave::err_nonconformant ("mx_el_or", m1_nr, m1_nc, m2_nr, m2_nc); \
1842 } \
1843 else if (m1_nr == m2_nr && m1_nc == m2_nc) \
1844 { \
1845 r = SparseBoolMatrix (m1_nr, m1_nc, m1.nnz () + m2.nnz ()); \
1846 r.cidx (0) = static_cast<octave_idx_type> (0); \
1847 octave_idx_type nel = 0; \
1848 for (octave_idx_type j = 0; j < m1_nc; j++) \
1849 { \
1850 octave_idx_type i1 = m1.cidx (j); \
1851 octave_idx_type e1 = m1.cidx (j+1); \
1852 octave_idx_type i2 = m2.cidx (j); \
1853 octave_idx_type e2 = m2.cidx (j+1); \
1854 while (i1 < e1 || i2 < e2) \
1855 { \
1856 octave_quit (); \
1857 \
1858 if (i1 == e1 || (i2 < e2 && m1.ridx (i1) > m2.ridx (i2))) \
1859 { \
1860 if (m2.data (i2) != rhs_zero) \
1861 { \
1862 r.ridx (nel) = m2.ridx (i2); \
1863 r.data (nel++) = true; \
1864 } \
1865 i2++; \
1866 } \
1867 else if (i2 == e2 || m1.ridx (i1) < m2.ridx (i2)) \
1868 { \
1869 if (m1.data (i1) != lhs_zero) \
1870 { \
1871 r.ridx (nel) = m1.ridx (i1); \
1872 r.data (nel++) = true; \
1873 } \
1874 i1++; \
1875 } \
1876 else \
1877 { \
1878 if (m1.data (i1) != lhs_zero || m2.data (i2) != rhs_zero) \
1879 { \
1880 r.ridx (nel) = m1.ridx (i1); \
1881 r.data (nel++) = true; \
1882 } \
1883 i1++; \
1884 i2++; \
1885 } \
1886 } \
1887 r.cidx (j + 1) = nel; \
1888 } \
1889 r.maybe_compress (false); \
1890 } \
1891 else if (m1_nr == m2_nr && (m1_nc == 1 || m2_nc == 1)) \
1892 { \
1893 octave_idx_type r_nc = (m1_nc < m2_nc ? m2_nc : m1_nc); \
1894 octave_idx_type rnnz = (m1_nc < m2_nc ? \
1895 m1.nnz () * m2_nc + m2.nnz () : \
1896 m1.nnz () + m1_nc * m2.nnz ()); \
1897 r = SparseBoolMatrix (m1_nr, r_nc, rnnz); \
1898 \
1899 r.cidx (0) = static_cast<octave_idx_type> (0); \
1900 octave_idx_type nel = 0; \
1901 for (octave_idx_type j = 0; j < r_nc; j++) \
1902 { \
1903 octave_idx_type i1; \
1904 octave_idx_type e1; \
1905 octave_idx_type i2; \
1906 octave_idx_type e2; \
1907 if (m1_nc == 1) \
1908 { \
1909 i1 = m1.cidx(0); \
1910 e1 = m1.cidx(1); \
1911 i2 = m2.cidx (j); \
1912 e2 = m2.cidx (j+1); \
1913 } \
1914 else \
1915 { \
1916 i1 = m1.cidx(j); \
1917 e1 = m1.cidx(j+1); \
1918 i2 = m2.cidx (0); \
1919 e2 = m2.cidx (1); \
1920 } \
1921 while (i1 < e1 || i2 < e2) \
1922 { \
1923 octave_quit (); \
1924 \
1925 if (i1 == e1 || (i2 < e2 && m1.ridx (i1) > m2.ridx (i2))) \
1926 { \
1927 if (m2.data (i2) != rhs_zero) \
1928 { \
1929 r.ridx (nel) = m2.ridx (i2); \
1930 r.data (nel++) = true; \
1931 } \
1932 i2++; \
1933 } \
1934 else if (i2 == e2 || m1.ridx (i1) < m2.ridx (i2)) \
1935 { \
1936 if (m1.data (i1) != lhs_zero) \
1937 { \
1938 r.ridx (nel) = m1.ridx (i1); \
1939 r.data (nel++) = true; \
1940 } \
1941 i1++; \
1942 } \
1943 else \
1944 { \
1945 if (m1.data (i1) != lhs_zero || m2.data (i2) != rhs_zero) \
1946 { \
1947 r.ridx (nel) = m1.ridx (i1); \
1948 r.data (nel++) = true; \
1949 } \
1950 i1++; \
1951 i2++; \
1952 } \
1953 } \
1954 r.cidx (j + 1) = nel; \
1955 } \
1956 r.maybe_compress (false); \
1957 } \
1958 else if (m1_nc == m2_nc && (m1_nr == 1 || m2_nr == 1)) \
1959 r = mx_el_or (m1.transpose (), m2.transpose ()).transpose (); \
1960 else if (m1_nr == 1 && m2_nc == 1) \
1961 { \
1962 /* Count num of nonzero elements */ \
1963 octave_idx_type nel = 0; \
1964 for (octave_idx_type j = 0; j < m1_nc; j++) \
1965 for (octave_idx_type i = 0; i < m2_nr; i++) \
1966 if ((m1.elem (0, j) != lhs_zero) || \
1967 (m2.elem (i, 0) != rhs_zero)) \
1968 nel++; \
1969 \
1970 r = SparseBoolMatrix (m2_nr, m1_nc, nel); \
1971 \
1972 for (octave_idx_type j = 0; j < m1_nc; j++) \
1973 { \
1974 const bool m1_val = m1.elem (0, j) != lhs_zero; \
1975 \
1976 for (octave_idx_type i = 0; i < m2_nr; i++) \
1977 { \
1978 octave_quit (); \
1979 \
1980 const bool m2_val = m2.elem (i, 0) != rhs_zero; \
1981 \
1982 if (m1_val || m2_val) \
1983 r.elem (i, j) = true; \
1984 } \
1985 } \
1986 r.maybe_compress (false); \
1987 } \
1988 else if (m1_nc == 1 && m2_nr == 1) \
1989 r = mx_el_or (m1.transpose (), m2.transpose ()).transpose (); \
1990 else \
1991 octave::err_nonconformant ("mx_el_or", m1_nr, m1_nc, m2_nr, m2_nc); \
1992 \
1993 return r; \
1994 }
1995
1996#define SPARSE_SMSM_BOOL_OPS(M1, M2) \
1997 SPARSE_SMSM_BOOL_AND_OP (M1, M2) \
1998 SPARSE_SMSM_BOOL_OR_OP (M1, M2)
1999
2000// matrix by sparse matrix operations.
2001
2002#define SPARSE_MSM_BIN_OP(R, F, OP, M1, M2) \
2003 R \
2004 F (const M1& m1, const M2& m2) \
2005 { \
2006 R r; \
2007 \
2008 octave_idx_type m1_nr = m1.rows (); \
2009 octave_idx_type m1_nc = m1.cols (); \
2010 \
2011 octave_idx_type m2_nr = m2.rows (); \
2012 octave_idx_type m2_nc = m2.cols (); \
2013 \
2014 if (m2_nr == 1 && m2_nc == 1) \
2015 r = R (m1 OP m2.elem (0,0)); \
2016 else if ((m1_nr == m2_nr && m1_nc == m2_nc) || \
2017 (m1_nr == m2_nr && (m1_nc == 1 || m2_nc == 1)) || \
2018 (m1_nc == m2_nc && (m1_nr == 1 || m2_nr == 1)) || \
2019 (m1_nr == 1 && m2_nc == 1) || \
2020 (m1_nc == 1 && m2_nr == 1)) \
2021 r = R (F (m1, m2.matrix_value ())); \
2022 else \
2023 octave::err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \
2024 \
2025 return r; \
2026 }
2027
2028#define SPARSE_MSM_BIN_OPS(R1, R2, M1, M2) \
2029 SPARSE_MSM_BIN_OP (R1, operator +, +, M1, M2) \
2030 SPARSE_MSM_BIN_OP (R1, operator -, -, M1, M2) \
2031 SPARSE_MSM_BIN_OP (R2, product, *, M1, M2) \
2032 SPARSE_MSM_BIN_OP (R1, quotient, /, M1, M2)
2033
2034#define SPARSE_MSM_CMP_OP(F, OP, M1, M2) \
2035 SparseBoolMatrix \
2036 F (const M1& m1, const M2& m2) \
2037 { \
2038 SparseBoolMatrix r; \
2039 \
2040 octave_idx_type m1_nr = m1.rows (); \
2041 octave_idx_type m1_nc = m1.cols (); \
2042 \
2043 octave_idx_type m2_nr = m2.rows (); \
2044 octave_idx_type m2_nc = m2.cols (); \
2045 \
2046 if (m2_nr == 1 && m2_nc == 1) \
2047 r = SparseBoolMatrix (F (m1, m2.elem (0,0))); \
2048 else if (m1_nr == 0 && (m2_nr == 0 || m2_nr == 1)) \
2049 { \
2050 if (m1_nc == 1 || m2_nc == 1 || m1_nc == m2_nc) \
2051 r.resize (m1_nr, std::max (m1_nc, m2_nc)); \
2052 else \
2053 octave::err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \
2054 } \
2055 else if (m1_nc == 0 && (m2_nc == 0 || m2_nc == 1)) \
2056 { \
2057 if (m1_nr == 1 || m2_nr == 1 || m1_nr == m2_nr) \
2058 r.resize (std::max (m1_nr, m2_nr), m1_nc); \
2059 else \
2060 octave::err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \
2061 } \
2062 else if (m2_nr == 0 && (m1_nr == 0 || m1_nr == 1)) \
2063 { \
2064 if (m2_nc == 1 || m1_nc == 1 || m2_nc == m1_nc) \
2065 r.resize (m2_nr, std::max (m1_nc, m2_nc)); \
2066 else \
2067 octave::err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \
2068 } \
2069 else if (m2_nc == 0 && (m1_nc == 0 || m1_nc == 1)) \
2070 { \
2071 if (m2_nr == 1 || m1_nr == 1 || m2_nr == m1_nr) \
2072 r.resize (std::max (m1_nr, m2_nr), m2_nc); \
2073 else \
2074 octave::err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \
2075 } \
2076 else if (m1_nr == m2_nr && m1_nc == m2_nc) \
2077 { \
2078 /* Count num of nonzero elements */ \
2079 octave_idx_type nel = 0; \
2080 for (octave_idx_type j = 0; j < m1_nc; j++) \
2081 for (octave_idx_type i = 0; i < m1_nr; i++) \
2082 if (m1.elem (i, j) OP m2.elem (i, j)) \
2083 nel++; \
2084 \
2085 r = SparseBoolMatrix (m1_nr, m1_nc, nel); \
2086 \
2087 octave_idx_type ii = 0; \
2088 r.cidx (0) = 0; \
2089 for (octave_idx_type j = 0; j < m1_nc; j++) \
2090 { \
2091 for (octave_idx_type i = 0; i < m1_nr; i++) \
2092 { \
2093 bool el = m1.elem (i, j) OP m2.elem (i, j); \
2094 if (el) \
2095 { \
2096 r.data (ii) = el; \
2097 r.ridx (ii++) = i; \
2098 } \
2099 } \
2100 r.cidx (j+1) = ii; \
2101 } \
2102 } \
2103 else if (m1_nr == m2_nr && (m1_nc == 1 || m2_nc == 1)) \
2104 { \
2105 octave_idx_type r_nc = (m1_nc < m2_nc ? m2_nc : m1_nc); \
2106 \
2107 /* Count num of nonzero elements */ \
2108 octave_idx_type nel = 0; \
2109 if (m1_nc == 1) \
2110 { \
2111 for (octave_idx_type j = 0; j < m2_nc; j++) \
2112 for (octave_idx_type i = 0; i < m1_nr; i++) \
2113 if (m1.elem (i, 0) OP m2.elem (i, j)) \
2114 nel++; \
2115 } \
2116 else \
2117 { \
2118 for (octave_idx_type j = 0; j < m1_nc; j++) \
2119 for (octave_idx_type i = 0; i < m1_nr; i++) \
2120 if (m1.elem (i, j) OP m2.elem (i, 0)) \
2121 nel++; \
2122 } \
2123 \
2124 r = SparseBoolMatrix (m1_nr, r_nc, nel); \
2125 \
2126 octave_idx_type ii = 0; \
2127 r.cidx (0) = 0; \
2128 if (m1_nc == 1) \
2129 { \
2130 for (octave_idx_type j = 0; j < m2_nc; j++) \
2131 { \
2132 for (octave_idx_type i = 0; i < m1_nr; i++) \
2133 { \
2134 bool el = m1.elem (i, 0) OP m2.elem (i, j); \
2135 if (el) \
2136 { \
2137 r.data (ii) = el; \
2138 r.ridx (ii++) = i; \
2139 } \
2140 } \
2141 r.cidx (j+1) = ii; \
2142 } \
2143 } \
2144 else \
2145 { \
2146 for (octave_idx_type j = 0; j < m1_nc; j++) \
2147 { \
2148 for (octave_idx_type i = 0; i < m1_nr; i++) \
2149 { \
2150 bool el = m1.elem (i, j) OP m2.elem (i, 0); \
2151 if (el) \
2152 { \
2153 r.data (ii) = el; \
2154 r.ridx (ii++) = i; \
2155 } \
2156 } \
2157 r.cidx (j+1) = ii; \
2158 } \
2159 } \
2160 r.maybe_compress (false); \
2161 } \
2162 else if (m1_nc == m2_nc && (m1_nr == 1 || m2_nr == 1)) \
2163 r = F (m1.transpose (), m2.transpose ()).transpose (); \
2164 else if (m1_nr == 1 && m2_nc == 1) \
2165 { \
2166 /* Count num of nonzero elements */ \
2167 octave_idx_type nel = 0; \
2168 for (octave_idx_type j = 0; j < m1_nc; j++) \
2169 for (octave_idx_type i = 0; i < m2_nr; i++) \
2170 if (m1.elem (0, j) OP m2.elem (i, 0)) \
2171 nel++; \
2172 \
2173 r = SparseBoolMatrix (m2_nr, m1_nc, nel); \
2174 \
2175 for (octave_idx_type j = 0; j < m1_nc; j++) \
2176 { \
2177 const Complex m1_val = m1.elem (0, j); \
2178 \
2179 for (octave_idx_type i = 0; i < m2_nr; i++) \
2180 { \
2181 octave_quit (); \
2182 \
2183 const Complex m2_val = m2.elem (i, 0); \
2184 \
2185 if (m1_val OP m2_val) \
2186 r.elem (i, j) = true; \
2187 } \
2188 } \
2189 r.maybe_compress (false); \
2190 } \
2191 else if (m1_nc == 1 && m2_nr == 1) \
2192 r = F (m1.transpose (), m2.transpose ()).transpose (); \
2193 else \
2194 octave::err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \
2195 \
2196 return r; \
2197 }
2198
2199#define SPARSE_MSM_CMP_OPS(M1, M2) \
2200 SPARSE_MSM_CMP_OP (mx_el_lt, <, M1, M2) \
2201 SPARSE_MSM_CMP_OP (mx_el_le, <=, M1, M2) \
2202 SPARSE_MSM_CMP_OP (mx_el_ge, >=, M1, M2) \
2203 SPARSE_MSM_CMP_OP (mx_el_gt, >, M1, M2) \
2204 SPARSE_MSM_CMP_OP (mx_el_eq, ==, M1, M2) \
2205 SPARSE_MSM_CMP_OP (mx_el_ne, !=, M1, M2)
2206
2207#define SPARSE_MSM_EQNE_OPS(M1, M2) \
2208 SPARSE_MSM_CMP_OP (mx_el_eq, ==, M1, M2) \
2209 SPARSE_MSM_CMP_OP (mx_el_ne, !=, M1, M2)
2210
2211#define SPARSE_MSM_BOOL_AND_OP(F, OP, M1, M2) \
2212 SparseBoolMatrix \
2213 F (const M1& m1, const M2& m2) \
2214 { \
2215 SparseBoolMatrix r; \
2216 \
2217 octave_idx_type m1_nr = m1.rows (); \
2218 octave_idx_type m1_nc = m1.cols (); \
2219 \
2220 octave_idx_type m2_nr = m2.rows (); \
2221 octave_idx_type m2_nc = m2.cols (); \
2222 \
2223 M1::element_type lhs_zero = M1::element_type (); \
2224 M2::element_type rhs_zero = M2::element_type (); \
2225 \
2226 if (m2_nr == 1 && m2_nc == 1) \
2227 r = SparseBoolMatrix (F (m1, m2.elem (0,0))); \
2228 else if (m1_nr == 0 && (m2_nr == 0 || m2_nr == 1)) \
2229 { \
2230 if (m1_nc == 1 || m2_nc == 1 || m1_nc == m2_nc) \
2231 r.resize (m1_nr, std::max (m1_nc, m2_nc)); \
2232 else \
2233 octave::err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \
2234 } \
2235 else if (m1_nc == 0 && (m2_nc == 0 || m2_nc == 1)) \
2236 { \
2237 if (m1_nr == 1 || m2_nr == 1 || m1_nr == m2_nr) \
2238 r.resize (std::max (m1_nr, m2_nr), m1_nc); \
2239 else \
2240 octave::err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \
2241 } \
2242 else if (m2_nr == 0 && (m1_nr == 0 || m1_nr == 1)) \
2243 { \
2244 if (m2_nc == 1 || m1_nc == 1 || m2_nc == m1_nc) \
2245 r.resize (m2_nr, std::max (m1_nc, m2_nc)); \
2246 else \
2247 octave::err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \
2248 } \
2249 else if (m2_nc == 0 && (m1_nc == 0 || m1_nc == 1)) \
2250 { \
2251 if (m2_nr == 1 || m1_nr == 1 || m2_nr == m1_nr) \
2252 r.resize (std::max (m1_nr, m2_nr), m2_nc); \
2253 else \
2254 octave::err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \
2255 } \
2256 else if (m1_nr == m2_nr && m1_nc == m2_nc) \
2257 { \
2258 /* Count num of nonzero elements */ \
2259 octave_idx_type nel = 0; \
2260 for (octave_idx_type j = 0; j < m1_nc; j++) \
2261 for (octave_idx_type i = 0; i < m1_nr; i++) \
2262 if ((m1.elem (i, j) != lhs_zero) OP \
2263 (m2.elem (i, j) != rhs_zero)) \
2264 nel++; \
2265 \
2266 r = SparseBoolMatrix (m1_nr, m1_nc, nel); \
2267 \
2268 octave_idx_type ii = 0; \
2269 r.cidx (0) = 0; \
2270 for (octave_idx_type j = 0; j < m1_nc; j++) \
2271 { \
2272 for (octave_idx_type i = 0; i < m1_nr; i++) \
2273 { \
2274 bool el = (m1.elem (i, j) != lhs_zero) OP \
2275 (m2.elem (i, j) != rhs_zero); \
2276 if (el) \
2277 { \
2278 r.data (ii) = el; \
2279 r.ridx (ii++) = i; \
2280 } \
2281 } \
2282 r.cidx (j+1) = ii; \
2283 } \
2284 } \
2285 else if (m1_nr == m2_nr && (m1_nc == 1 || m2_nc == 1)) \
2286 { \
2287 octave_idx_type r_nc = (m1_nc < m2_nc ? m2_nc : m1_nc); \
2288 \
2289 /* Count num of nonzero elements */ \
2290 octave_idx_type nel = 0; \
2291 if (m1_nc == 1) \
2292 { \
2293 for (octave_idx_type j = 0; j < m2_nc; j++) \
2294 for (octave_idx_type i = 0; i < m1_nr; i++) \
2295 if ((m1.elem (i, 0) != lhs_zero) OP \
2296 (m2.elem (i, j) != rhs_zero)) \
2297 nel++; \
2298 } \
2299 else \
2300 { \
2301 for (octave_idx_type j = 0; j < m1_nc; j++) \
2302 for (octave_idx_type i = 0; i < m1_nr; i++) \
2303 if ((m1.elem (i, j) != lhs_zero) OP \
2304 (m2.elem (i, 0) != rhs_zero)) \
2305 nel++; \
2306 } \
2307 \
2308 r = SparseBoolMatrix (m1_nr, r_nc, nel); \
2309 \
2310 octave_idx_type ii = 0; \
2311 r.cidx (0) = 0; \
2312 if (m1_nc == 1) \
2313 { \
2314 for (octave_idx_type j = 0; j < m2_nc; j++) \
2315 { \
2316 for (octave_idx_type i = 0; i < m1_nr; i++) \
2317 { \
2318 bool el = (m1.elem (i, 0) != lhs_zero) OP \
2319 (m2.elem (i, j) != rhs_zero); \
2320 if (el) \
2321 { \
2322 r.data (ii) = el; \
2323 r.ridx (ii++) = i; \
2324 } \
2325 } \
2326 r.cidx (j+1) = ii; \
2327 } \
2328 } \
2329 else \
2330 { \
2331 for (octave_idx_type j = 0; j < m1_nc; j++) \
2332 { \
2333 for (octave_idx_type i = 0; i < m1_nr; i++) \
2334 { \
2335 bool el = (m1.elem (i, j) != lhs_zero) OP \
2336 (m2.elem (i, 0) != rhs_zero); \
2337 if (el) \
2338 { \
2339 r.data (ii) = el; \
2340 r.ridx (ii++) = i; \
2341 } \
2342 } \
2343 r.cidx (j+1) = ii; \
2344 } \
2345 } \
2346 r.maybe_compress (false); \
2347 } \
2348 else if (m1_nc == m2_nc && (m1_nr == 1 || m2_nr == 1)) \
2349 r = F (m1.transpose (), m2.transpose ()).transpose (); \
2350 else if (m1_nr == 1 && m2_nc == 1) \
2351 { \
2352 /* Count num of nonzero elements */ \
2353 octave_idx_type nel = 0; \
2354 for (octave_idx_type j = 0; j < m1_nc; j++) \
2355 for (octave_idx_type i = 0; i < m2_nr; i++) \
2356 if ((m1.elem (0, j) != lhs_zero) OP \
2357 (m2.elem (i, 0) != rhs_zero)) \
2358 nel++; \
2359 \
2360 r = SparseBoolMatrix (m2_nr, m1_nc, nel); \
2361 \
2362 for (octave_idx_type j = 0; j < m1_nc; j++) \
2363 { \
2364 const bool m1_val = m1.elem (0, j) != lhs_zero; \
2365 \
2366 for (octave_idx_type i = 0; i < m2_nr; i++) \
2367 { \
2368 octave_quit (); \
2369 \
2370 const bool m2_val = m2.elem (i, 0) != lhs_zero; \
2371 \
2372 if (m1_val OP m2_val) \
2373 r.elem (i, j) = true; \
2374 } \
2375 } \
2376 r.maybe_compress (false); \
2377 } \
2378 else if (m1_nc == 1 && m2_nr == 1) \
2379 r = F (m1.transpose (), m2.transpose ()).transpose (); \
2380 else \
2381 octave::err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \
2382 \
2383 return r; \
2384 }
2385
2386#define SPARSE_MSM_BOOL_OR_OP(F, OP, M1, M2) \
2387 boolMatrix \
2388 F (const M1& m1, const M2& m2) \
2389 { \
2390 boolMatrix r; \
2391 \
2392 octave_idx_type m1_nr = m1.rows (); \
2393 octave_idx_type m1_nc = m1.cols (); \
2394 \
2395 octave_idx_type m2_nr = m2.rows (); \
2396 octave_idx_type m2_nc = m2.cols (); \
2397 \
2398 M1::element_type lhs_zero = M1::element_type (); \
2399 M2::element_type rhs_zero = M2::element_type (); \
2400 \
2401 if (m2_nr == 1 && m2_nc == 1) \
2402 r = boolMatrix (F (m1, m2.elem (0,0))); \
2403 else if (m1_nr == 0 && (m2_nr == 0 || m2_nr == 1)) \
2404 { \
2405 if (m1_nc == 1 || m2_nc == 1 || m1_nc == m2_nc) \
2406 r.resize (m1_nr, std::max (m1_nc, m2_nc)); \
2407 else \
2408 octave::err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \
2409 } \
2410 else if (m1_nc == 0 && (m2_nc == 0 || m2_nc == 1)) \
2411 { \
2412 if (m1_nr == 1 || m2_nr == 1 || m1_nr == m2_nr) \
2413 r.resize (std::max (m1_nr, m2_nr), m1_nc); \
2414 else \
2415 octave::err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \
2416 } \
2417 else if (m2_nr == 0 && (m1_nr == 0 || m1_nr == 1)) \
2418 { \
2419 if (m2_nc == 1 || m1_nc == 1 || m2_nc == m1_nc) \
2420 r.resize (m2_nr, std::max (m1_nc, m2_nc)); \
2421 else \
2422 octave::err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \
2423 } \
2424 else if (m2_nc == 0 && (m1_nc == 0 || m1_nc == 1)) \
2425 { \
2426 if (m2_nr == 1 || m1_nr == 1 || m2_nr == m1_nr) \
2427 r.resize (std::max (m1_nr, m2_nr), m2_nc); \
2428 else \
2429 octave::err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \
2430 } \
2431 else if (m1_nr == m2_nr && m1_nc == m2_nc) \
2432 { \
2433 r = boolMatrix (m1_nr, m1_nc, false); \
2434 \
2435 for (octave_idx_type j = 0; j < m1_nc; j++) \
2436 for (octave_idx_type i = 0; i < m1_nr; i++) \
2437 { \
2438 bool el = (m1.elem (i, j) != lhs_zero) OP \
2439 (m2.elem (i, j) != rhs_zero); \
2440 if (el) \
2441 r.elem (i, j) = el; \
2442 } \
2443 } \
2444 else if (m1_nr == m2_nr && (m1_nc == 1 || m2_nc == 1)) \
2445 { \
2446 octave_idx_type r_nc = (m1_nc < m2_nc ? m2_nc : m1_nc); \
2447 \
2448 r = boolMatrix (m1_nr, r_nc, false); \
2449 \
2450 if (m1_nc == 1) \
2451 { \
2452 for (octave_idx_type j = 0; j < m2_nc; j++) \
2453 for (octave_idx_type i = 0; i < m1_nr; i++) \
2454 { \
2455 bool el = (m1.elem (i, 0) != lhs_zero) OP \
2456 (m2.elem (i, j) != rhs_zero); \
2457 if (el) \
2458 r.elem (i, j) = el; \
2459 } \
2460 } \
2461 else \
2462 { \
2463 for (octave_idx_type j = 0; j < m1_nc; j++) \
2464 for (octave_idx_type i = 0; i < m1_nr; i++) \
2465 { \
2466 bool el = (m1.elem (i, j) != lhs_zero) OP \
2467 (m2.elem (i, 0) != rhs_zero); \
2468 if (el) \
2469 r.elem (i, j) = el; \
2470 } \
2471 } \
2472 } \
2473 else if (m1_nc == m2_nc && (m1_nr == 1 || m2_nr == 1)) \
2474 r = F (m1.transpose (), m2.transpose ()).transpose (); \
2475 else if (m1_nr == 1 && m2_nc == 1) \
2476 { \
2477 r = boolMatrix (m2_nr, m1_nc, false); \
2478 \
2479 for (octave_idx_type j = 0; j < m1_nc; j++) \
2480 { \
2481 const bool m1_val = m1.elem (0, j) != lhs_zero; \
2482 \
2483 for (octave_idx_type i = 0; i < m2_nr; i++) \
2484 { \
2485 octave_quit (); \
2486 \
2487 const bool m2_val = m2.elem (i, 0) != lhs_zero; \
2488 \
2489 if (m1_val OP m2_val) \
2490 r.elem (i, j) = true; \
2491 } \
2492 } \
2493 } \
2494 else if (m1_nc == 1 && m2_nr == 1) \
2495 r = F (m1.transpose (), m2.transpose ()).transpose (); \
2496 else \
2497 octave::err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \
2498 \
2499 return r; \
2500 }
2501
2502#define SPARSE_MSM_BOOL_OPS(M1, M2) \
2503 SPARSE_MSM_BOOL_AND_OP (mx_el_and, &&, M1, M2) \
2504 SPARSE_MSM_BOOL_OR_OP (mx_el_or, ||, M1, M2)
2505
2506// sparse matrix by matrix operations.
2507
2508#define SPARSE_SMM_BIN_OP(R, F, OP, M1, M2) \
2509 R \
2510 F (const M1& m1, const M2& m2) \
2511 { \
2512 R r; \
2513 \
2514 octave_idx_type m1_nr = m1.rows (); \
2515 octave_idx_type m1_nc = m1.cols (); \
2516 \
2517 octave_idx_type m2_nr = m2.rows (); \
2518 octave_idx_type m2_nc = m2.cols (); \
2519 \
2520 if (m1_nr == 1 && m1_nc == 1) \
2521 r = R (m1.elem (0,0) OP m2); \
2522 else if ((m1_nr == m2_nr && m1_nc == m2_nc) || \
2523 (m1_nr == m2_nr && (m1_nc == 1 || m2_nc == 1)) || \
2524 (m1_nc == m2_nc && (m1_nr == 1 || m2_nr == 1)) || \
2525 (m1_nr == 1 && m2_nc == 1) || \
2526 (m1_nc == 1 && m2_nr == 1)) \
2527 r = R (F (m1.matrix_value (), m2)); \
2528 else \
2529 octave::err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \
2530 \
2531 return r; \
2532 }
2533
2534#define SPARSE_SMM_BIN_OPS(R1, R2, M1, M2) \
2535 SPARSE_SMM_BIN_OP (R1, operator +, +, M1, M2) \
2536 SPARSE_SMM_BIN_OP (R1, operator -, -, M1, M2) \
2537 SPARSE_SMM_BIN_OP (R2, product, *, M1, M2) \
2538 SPARSE_SMM_BIN_OP (R2, quotient, /, M1, M2)
2539
2540#define SPARSE_SMM_CMP_OP(F, OP, M1, M2) \
2541 SparseBoolMatrix \
2542 F (const M1& m1, const M2& m2) \
2543 { \
2544 SparseBoolMatrix r; \
2545 \
2546 octave_idx_type m1_nr = m1.rows (); \
2547 octave_idx_type m1_nc = m1.cols (); \
2548 \
2549 octave_idx_type m2_nr = m2.rows (); \
2550 octave_idx_type m2_nc = m2.cols (); \
2551 \
2552 if (m1_nr == 1 && m1_nc == 1) \
2553 r = SparseBoolMatrix (F (m1.elem (0,0), m2)); \
2554 else if (m1_nr == 0 && (m2_nr == 0 || m2_nr == 1)) \
2555 { \
2556 if (m1_nc == 1 || m2_nc == 1 || m1_nc == m2_nc) \
2557 r.resize (m1_nr, std::max (m1_nc, m2_nc)); \
2558 else \
2559 octave::err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \
2560 } \
2561 else if (m1_nc == 0 && (m2_nc == 0 || m2_nc == 1)) \
2562 { \
2563 if (m1_nr == 1 || m2_nr == 1 || m1_nr == m2_nr) \
2564 r.resize (std::max (m1_nr, m2_nr), m1_nc); \
2565 else \
2566 octave::err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \
2567 } \
2568 else if (m2_nr == 0 && (m1_nr == 0 || m1_nr == 1)) \
2569 { \
2570 if (m2_nc == 1 || m1_nc == 1 || m2_nc == m1_nc) \
2571 r.resize (m2_nr, std::max (m1_nc, m2_nc)); \
2572 else \
2573 octave::err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \
2574 } \
2575 else if (m2_nc == 0 && (m1_nc == 0 || m1_nc == 1)) \
2576 { \
2577 if (m2_nr == 1 || m1_nr == 1 || m2_nr == m1_nr) \
2578 r.resize (std::max (m1_nr, m2_nr), m2_nc); \
2579 else \
2580 octave::err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \
2581 } \
2582 else if (m1_nr == m2_nr && m1_nc == m2_nc) \
2583 { \
2584 /* Count num of nonzero elements */ \
2585 octave_idx_type nel = 0; \
2586 for (octave_idx_type j = 0; j < m1_nc; j++) \
2587 for (octave_idx_type i = 0; i < m1_nr; i++) \
2588 if (m1.elem (i, j) OP m2.elem (i, j)) \
2589 nel++; \
2590 \
2591 r = SparseBoolMatrix (m1_nr, m1_nc, nel); \
2592 \
2593 octave_idx_type ii = 0; \
2594 r.cidx (0) = 0; \
2595 for (octave_idx_type j = 0; j < m1_nc; j++) \
2596 { \
2597 for (octave_idx_type i = 0; i < m1_nr; i++) \
2598 { \
2599 bool el = m1.elem (i, j) OP m2.elem (i, j); \
2600 if (el) \
2601 { \
2602 r.data (ii) = el; \
2603 r.ridx (ii++) = i; \
2604 } \
2605 } \
2606 r.cidx (j+1) = ii; \
2607 } \
2608 } \
2609 else if (m1_nr == m2_nr && (m1_nc == 1 || m2_nc == 1)) \
2610 { \
2611 octave_idx_type r_nc = (m1_nc < m2_nc ? m2_nc : m1_nc); \
2612 \
2613 /* Count num of nonzero elements */ \
2614 octave_idx_type nel = 0; \
2615 if (m1_nc == 1) \
2616 { \
2617 for (octave_idx_type j = 0; j < m2_nc; j++) \
2618 for (octave_idx_type i = 0; i < m1_nr; i++) \
2619 if (m1.elem (i, 0) OP m2.elem (i, j)) \
2620 nel++; \
2621 } \
2622 else \
2623 { \
2624 for (octave_idx_type j = 0; j < m1_nc; j++) \
2625 for (octave_idx_type i = 0; i < m1_nr; i++) \
2626 if (m1.elem (i, j) OP m2.elem (i, 0)) \
2627 nel++; \
2628 } \
2629 \
2630 r = SparseBoolMatrix (m1_nr, r_nc, nel); \
2631 \
2632 octave_idx_type ii = 0; \
2633 r.cidx (0) = 0; \
2634 if (m1_nc == 1) \
2635 { \
2636 for (octave_idx_type j = 0; j < m2_nc; j++) \
2637 { \
2638 for (octave_idx_type i = 0; i < m1_nr; i++) \
2639 { \
2640 bool el = m1.elem (i, 0) OP m2.elem (i, j); \
2641 if (el) \
2642 { \
2643 r.data (ii) = el; \
2644 r.ridx (ii++) = i; \
2645 } \
2646 } \
2647 r.cidx (j+1) = ii; \
2648 } \
2649 } \
2650 else \
2651 { \
2652 for (octave_idx_type j = 0; j < m1_nc; j++) \
2653 { \
2654 for (octave_idx_type i = 0; i < m1_nr; i++) \
2655 { \
2656 bool el = m1.elem (i, j) OP m2.elem (i, 0); \
2657 if (el) \
2658 { \
2659 r.data (ii) = el; \
2660 r.ridx (ii++) = i; \
2661 } \
2662 } \
2663 r.cidx (j+1) = ii; \
2664 } \
2665 } \
2666 r.maybe_compress (false); \
2667 } \
2668 else if (m1_nc == m2_nc && (m1_nr == 1 || m2_nr == 1)) \
2669 r = F (m1.transpose (), m2.transpose ()).transpose (); \
2670 else if (m1_nr == 1 && m2_nc == 1) \
2671 { \
2672 /* Count num of nonzero elements */ \
2673 octave_idx_type nel = 0; \
2674 for (octave_idx_type j = 0; j < m1_nc; j++) \
2675 for (octave_idx_type i = 0; i < m2_nr; i++) \
2676 if (m1.elem (0, j) OP m2.elem (i, 0)) \
2677 nel++; \
2678 \
2679 r = SparseBoolMatrix (m2_nr, m1_nc, nel); \
2680 \
2681 for (octave_idx_type j = 0; j < m1_nc; j++) \
2682 { \
2683 const Complex m1_val = m1.elem (0, j); \
2684 \
2685 for (octave_idx_type i = 0; i < m2_nr; i++) \
2686 { \
2687 octave_quit (); \
2688 \
2689 const Complex m2_val = m2.elem (i, 0); \
2690 \
2691 if (m1_val OP m2_val) \
2692 r.elem (i, j) = true; \
2693 } \
2694 } \
2695 r.maybe_compress (false); \
2696 } \
2697 else if (m1_nc == 1 && m2_nr == 1) \
2698 r = F (m1.transpose (), m2.transpose ()).transpose (); \
2699 else \
2700 octave::err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \
2701 \
2702 return r; \
2703 }
2704
2705#define SPARSE_SMM_CMP_OPS(M1, M2) \
2706 SPARSE_SMM_CMP_OP (mx_el_lt, <, M1, M2) \
2707 SPARSE_SMM_CMP_OP (mx_el_le, <=, M1, M2) \
2708 SPARSE_SMM_CMP_OP (mx_el_ge, >=, M1, M2) \
2709 SPARSE_SMM_CMP_OP (mx_el_gt, >, M1, M2) \
2710 SPARSE_SMM_CMP_OP (mx_el_eq, ==, M1, M2) \
2711 SPARSE_SMM_CMP_OP (mx_el_ne, !=, M1, M2)
2712
2713#define SPARSE_SMM_EQNE_OPS(M1, M2) \
2714 SPARSE_SMM_CMP_OP (mx_el_eq, ==, M1, M2) \
2715 SPARSE_SMM_CMP_OP (mx_el_ne, !=, M1, M2)
2716
2717#define SPARSE_SMM_BOOL_AND_OP(F, OP, M1, M2) \
2718 SparseBoolMatrix \
2719 F (const M1& m1, const M2& m2) \
2720 { \
2721 SparseBoolMatrix r; \
2722 \
2723 octave_idx_type m1_nr = m1.rows (); \
2724 octave_idx_type m1_nc = m1.cols (); \
2725 \
2726 octave_idx_type m2_nr = m2.rows (); \
2727 octave_idx_type m2_nc = m2.cols (); \
2728 \
2729 M1::element_type lhs_zero = M1::element_type (); \
2730 M2::element_type rhs_zero = M2::element_type (); \
2731 \
2732 if (m1_nr == 1 && m1_nc == 1) \
2733 r = SparseBoolMatrix (F (m1.elem (0,0), m2)); \
2734 else if (m1_nr == 0 && (m2_nr == 0 || m2_nr == 1)) \
2735 { \
2736 if (m1_nc == 1 || m2_nc == 1 || m1_nc == m2_nc) \
2737 r.resize (m1_nr, std::max (m1_nc, m2_nc)); \
2738 else \
2739 octave::err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \
2740 } \
2741 else if (m1_nc == 0 && (m2_nc == 0 || m2_nc == 1)) \
2742 { \
2743 if (m1_nr == 1 || m2_nr == 1 || m1_nr == m2_nr) \
2744 r.resize (std::max (m1_nr, m2_nr), m1_nc); \
2745 else \
2746 octave::err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \
2747 } \
2748 else if (m2_nr == 0 && (m1_nr == 0 || m1_nr == 1)) \
2749 { \
2750 if (m2_nc == 1 || m1_nc == 1 || m2_nc == m1_nc) \
2751 r.resize (m2_nr, std::max (m1_nc, m2_nc)); \
2752 else \
2753 octave::err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \
2754 } \
2755 else if (m2_nc == 0 && (m1_nc == 0 || m1_nc == 1)) \
2756 { \
2757 if (m2_nr == 1 || m1_nr == 1 || m2_nr == m1_nr) \
2758 r.resize (std::max (m1_nr, m2_nr), m2_nc); \
2759 else \
2760 octave::err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \
2761 } \
2762 else if (m1_nr == m2_nr && m1_nc == m2_nc) \
2763 { \
2764 /* Count num of nonzero elements */ \
2765 octave_idx_type nel = 0; \
2766 for (octave_idx_type j = 0; j < m1_nc; j++) \
2767 for (octave_idx_type i = 0; i < m1_nr; i++) \
2768 if ((m1.elem (i, j) != lhs_zero) OP \
2769 (m2.elem (i, j) != rhs_zero)) \
2770 nel++; \
2771 \
2772 r = SparseBoolMatrix (m1_nr, m1_nc, nel); \
2773 \
2774 octave_idx_type ii = 0; \
2775 r.cidx (0) = 0; \
2776 for (octave_idx_type j = 0; j < m1_nc; j++) \
2777 { \
2778 for (octave_idx_type i = 0; i < m1_nr; i++) \
2779 { \
2780 bool el = (m1.elem (i, j) != lhs_zero) OP \
2781 (m2.elem (i, j) != rhs_zero); \
2782 if (el) \
2783 { \
2784 r.data (ii) = el; \
2785 r.ridx (ii++) = i; \
2786 } \
2787 } \
2788 r.cidx (j+1) = ii; \
2789 } \
2790 } \
2791 else if (m1_nr == m2_nr && (m1_nc == 1 || m2_nc == 1)) \
2792 { \
2793 octave_idx_type r_nc = (m1_nc < m2_nc ? m2_nc : m1_nc); \
2794 \
2795 /* Count num of nonzero elements */ \
2796 octave_idx_type nel = 0; \
2797 if (m1_nc == 1) \
2798 { \
2799 for (octave_idx_type j = 0; j < m2_nc; j++) \
2800 for (octave_idx_type i = 0; i < m1_nr; i++) \
2801 if ((m1.elem (i, 0) != lhs_zero) OP \
2802 (m2.elem (i, j) != rhs_zero)) \
2803 nel++; \
2804 } \
2805 else \
2806 { \
2807 for (octave_idx_type j = 0; j < m1_nc; j++) \
2808 for (octave_idx_type i = 0; i < m1_nr; i++) \
2809 if ((m1.elem (i, j) != lhs_zero) OP \
2810 (m2.elem (i, 0) != rhs_zero)) \
2811 nel++; \
2812 } \
2813 \
2814 r = SparseBoolMatrix (m1_nr, r_nc, nel); \
2815 \
2816 octave_idx_type ii = 0; \
2817 r.cidx (0) = 0; \
2818 if (m1_nc == 1) \
2819 { \
2820 for (octave_idx_type j = 0; j < m2_nc; j++) \
2821 { \
2822 for (octave_idx_type i = 0; i < m1_nr; i++) \
2823 { \
2824 bool el = (m1.elem (i, 0) != lhs_zero) OP \
2825 (m2.elem (i, j) != rhs_zero); \
2826 if (el) \
2827 { \
2828 r.data (ii) = el; \
2829 r.ridx (ii++) = i; \
2830 } \
2831 } \
2832 r.cidx (j+1) = ii; \
2833 } \
2834 } \
2835 else \
2836 { \
2837 for (octave_idx_type j = 0; j < m1_nc; j++) \
2838 { \
2839 for (octave_idx_type i = 0; i < m1_nr; i++) \
2840 { \
2841 bool el = (m1.elem (i, j) != lhs_zero) OP \
2842 (m2.elem (i, 0) != rhs_zero); \
2843 if (el) \
2844 { \
2845 r.data (ii) = el; \
2846 r.ridx (ii++) = i; \
2847 } \
2848 } \
2849 r.cidx (j+1) = ii; \
2850 } \
2851 } \
2852 r.maybe_compress (false); \
2853 } \
2854 else if (m1_nc == m2_nc && (m1_nr == 1 || m2_nr == 1)) \
2855 r = F (m1.transpose (), m2.transpose ()).transpose (); \
2856 else if (m1_nr == 1 && m2_nc == 1) \
2857 { \
2858 /* Count num of nonzero elements */ \
2859 octave_idx_type nel = 0; \
2860 for (octave_idx_type j = 0; j < m1_nc; j++) \
2861 for (octave_idx_type i = 0; i < m2_nr; i++) \
2862 if ((m1.elem (0, j) != lhs_zero) OP \
2863 (m2.elem (i, 0) != rhs_zero)) \
2864 nel++; \
2865 \
2866 r = SparseBoolMatrix (m2_nr, m1_nc, nel); \
2867 \
2868 for (octave_idx_type j = 0; j < m1_nc; j++) \
2869 { \
2870 const bool m1_val = m1.elem (0, j) != lhs_zero; \
2871 \
2872 for (octave_idx_type i = 0; i < m2_nr; i++) \
2873 { \
2874 octave_quit (); \
2875 \
2876 const bool m2_val = m2.elem (i, 0) != lhs_zero; \
2877 \
2878 if (m1_val OP m2_val) \
2879 r.elem (i, j) = true; \
2880 } \
2881 } \
2882 r.maybe_compress (false); \
2883 } \
2884 else if (m1_nc == 1 && m2_nr == 1) \
2885 r = F (m1.transpose (), m2.transpose ()).transpose (); \
2886 else \
2887 octave::err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \
2888 \
2889 return r; \
2890 }
2891
2892#define SPARSE_SMM_BOOL_OR_OP(F, OP, M1, M2) \
2893 boolMatrix \
2894 F (const M1& m1, const M2& m2) \
2895 { \
2896 boolMatrix r; \
2897 \
2898 octave_idx_type m1_nr = m1.rows (); \
2899 octave_idx_type m1_nc = m1.cols (); \
2900 \
2901 octave_idx_type m2_nr = m2.rows (); \
2902 octave_idx_type m2_nc = m2.cols (); \
2903 \
2904 M1::element_type lhs_zero = M1::element_type (); \
2905 M2::element_type rhs_zero = M2::element_type (); \
2906 \
2907 if (m1_nr == 1 && m1_nc == 1) \
2908 r = boolMatrix (F (m1.elem (0,0), m2)); \
2909 else if (m1_nr == 0 && (m2_nr == 0 || m2_nr == 1)) \
2910 { \
2911 if (m1_nc == 1 || m2_nc == 1 || m1_nc == m2_nc) \
2912 r.resize (m1_nr, std::max (m1_nc, m2_nc)); \
2913 else \
2914 octave::err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \
2915 } \
2916 else if (m1_nc == 0 && (m2_nc == 0 || m2_nc == 1)) \
2917 { \
2918 if (m1_nr == 1 || m2_nr == 1 || m1_nr == m2_nr) \
2919 r.resize (std::max (m1_nr, m2_nr), m1_nc); \
2920 else \
2921 octave::err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \
2922 } \
2923 else if (m2_nr == 0 && (m1_nr == 0 || m1_nr == 1)) \
2924 { \
2925 if (m2_nc == 1 || m1_nc == 1 || m2_nc == m1_nc) \
2926 r.resize (m2_nr, std::max (m1_nc, m2_nc)); \
2927 else \
2928 octave::err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \
2929 } \
2930 else if (m2_nc == 0 && (m1_nc == 0 || m1_nc == 1)) \
2931 { \
2932 if (m2_nr == 1 || m1_nr == 1 || m2_nr == m1_nr) \
2933 r.resize (std::max (m1_nr, m2_nr), m2_nc); \
2934 else \
2935 octave::err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \
2936 } \
2937 else if (m1_nr == m2_nr && m1_nc == m2_nc) \
2938 { \
2939 r = boolMatrix (m1_nr, m1_nc, false); \
2940 \
2941 for (octave_idx_type j = 0; j < m1_nc; j++) \
2942 for (octave_idx_type i = 0; i < m1_nr; i++) \
2943 { \
2944 bool el = (m1.elem (i, j) != lhs_zero) OP \
2945 (m2.elem (i, j) != rhs_zero); \
2946 if (el) \
2947 r.elem (i, j) = el; \
2948 } \
2949 } \
2950 else if (m1_nr == m2_nr && (m1_nc == 1 || m2_nc == 1)) \
2951 { \
2952 octave_idx_type r_nc = (m1_nc < m2_nc ? m2_nc : m1_nc); \
2953 \
2954 r = boolMatrix (m1_nr, r_nc, false); \
2955 \
2956 if (m1_nc == 1) \
2957 { \
2958 for (octave_idx_type j = 0; j < m2_nc; j++) \
2959 for (octave_idx_type i = 0; i < m1_nr; i++) \
2960 { \
2961 bool el = (m1.elem (i, 0) != lhs_zero) OP \
2962 (m2.elem (i, j) != rhs_zero); \
2963 if (el) \
2964 r.elem (i, j) = el; \
2965 } \
2966 } \
2967 else \
2968 { \
2969 for (octave_idx_type j = 0; j < m1_nc; j++) \
2970 for (octave_idx_type i = 0; i < m1_nr; i++) \
2971 { \
2972 bool el = (m1.elem (i, j) != lhs_zero) OP \
2973 (m2.elem (i, 0) != rhs_zero); \
2974 if (el) \
2975 r.elem (i, j) = el; \
2976 } \
2977 } \
2978 } \
2979 else if (m1_nc == m2_nc && (m1_nr == 1 || m2_nr == 1)) \
2980 r = F (m1.transpose (), m2.transpose ()).transpose (); \
2981 else if (m1_nr == 1 && m2_nc == 1) \
2982 { \
2983 r = boolMatrix (m2_nr, m1_nc, false); \
2984 \
2985 for (octave_idx_type j = 0; j < m1_nc; j++) \
2986 { \
2987 const bool m1_val = m1.elem (0, j) != lhs_zero; \
2988 \
2989 for (octave_idx_type i = 0; i < m2_nr; i++) \
2990 { \
2991 octave_quit (); \
2992 \
2993 const bool m2_val = m2.elem (i, 0) != lhs_zero; \
2994 \
2995 if (m1_val OP m2_val) \
2996 r.elem (i, j) = true; \
2997 } \
2998 } \
2999 } \
3000 else if (m1_nc == 1 && m2_nr == 1) \
3001 r = F (m1.transpose (), m2.transpose ()).transpose (); \
3002 else \
3003 octave::err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \
3004 \
3005 return r; \
3006 }
3007
3008#define SPARSE_SMM_BOOL_OPS(M1, M2) \
3009 SPARSE_SMM_BOOL_AND_OP (mx_el_and, &&, M1, M2) \
3010 SPARSE_SMM_BOOL_OR_OP (mx_el_or, ||, M1, M2)
3011
3012// 'cumsum' and 'cumprod' specializations for sparse matrices
3013
3014#define SPARSE_CUMSUM(RET_TYPE, ELT_TYPE, FCN, NAN_EXPR, EXPR) \
3015 \
3016 octave_idx_type nr = rows (); \
3017 octave_idx_type nc = cols (); \
3018 \
3019 RET_TYPE retval; \
3020 \
3021 if (nr > 0 && nc > 0) \
3022 { \
3023 if ((nr == 1 && dim == -1) || dim == 1) \
3024 /* Ugly!! Is there a better way? */ \
3025 retval = transpose (). FCN (0, nanflag) .transpose (); \
3026 else \
3027 { \
3028 octave_idx_type nel = 0; \
3029 for (octave_idx_type i = 0; i < nc; i++) \
3030 { \
3031 ELT_TYPE t = ELT_TYPE (); \
3032 for (octave_idx_type j = cidx (i); j < cidx (i+1); j++) \
3033 { \
3034 t += data (j); \
3035 if (t != ELT_TYPE ()) \
3036 { \
3037 if (j == cidx (i+1) - 1) \
3038 nel += nr - ridx (j); \
3039 else \
3040 nel += ridx (j+1) - ridx (j); \
3041 } \
3042 } \
3043 } \
3044 retval = RET_TYPE (nr, nc, nel); \
3045 retval.cidx (0) = 0; \
3046 octave_idx_type ii = 0; \
3047 if (nanflag) \
3048 { \
3049 for (octave_idx_type i = 0; i < nc; i++) \
3050 { \
3051 ELT_TYPE t = ELT_TYPE (); \
3052 for (octave_idx_type j = cidx (i); j < cidx (i+1); j++) \
3053 { \
3054 NAN_EXPR; \
3055 if (t != ELT_TYPE ()) \
3056 { \
3057 if (j == cidx (i+1) - 1) \
3058 { \
3059 for (octave_idx_type k = ridx (j); k < nr; k++) \
3060 { \
3061 retval.data (ii) = t; \
3062 retval.ridx (ii++) = k; \
3063 } \
3064 } \
3065 else \
3066 { \
3067 for (octave_idx_type k = ridx (j); \
3068 k < ridx (j+1); k++) \
3069 { \
3070 retval.data (ii) = t; \
3071 retval.ridx (ii++) = k; \
3072 } \
3073 } \
3074 } \
3075 } \
3076 retval.cidx (i+1) = ii; \
3077 } \
3078 } \
3079 else \
3080 { \
3081 for (octave_idx_type i = 0; i < nc; i++) \
3082 { \
3083 ELT_TYPE t = ELT_TYPE (); \
3084 for (octave_idx_type j = cidx (i); j < cidx (i+1); j++) \
3085 { \
3086 EXPR; \
3087 if (t != ELT_TYPE ()) \
3088 { \
3089 if (j == cidx (i+1) - 1) \
3090 { \
3091 for (octave_idx_type k = ridx (j); k < nr; k++) \
3092 { \
3093 retval.data (ii) = t; \
3094 retval.ridx (ii++) = k; \
3095 } \
3096 } \
3097 else \
3098 { \
3099 for (octave_idx_type k = ridx (j); \
3100 k < ridx (j+1); k++) \
3101 { \
3102 retval.data (ii) = t; \
3103 retval.ridx (ii++) = k; \
3104 } \
3105 } \
3106 } \
3107 } \
3108 retval.cidx (i+1) = ii; \
3109 } \
3110 } \
3111 } \
3112 } \
3113 else \
3114 retval = RET_TYPE (nr, nc); \
3115 \
3116 return retval
3117
3118#define SPARSE_CUMPROD(RET_TYPE, ELT_TYPE, FCN) \
3119 \
3120 octave_idx_type nr = rows (); \
3121 octave_idx_type nc = cols (); \
3122 \
3123 RET_TYPE retval; \
3124 \
3125 if (nr > 0 && nc > 0) \
3126 { \
3127 if ((nr == 1 && dim == -1) || dim == 1) \
3128 /* Ugly!! Is there a better way? */ \
3129 retval = transpose (). FCN (0, nanflag) .transpose (); \
3130 else \
3131 { \
3132 octave_idx_type nel = 0; \
3133 for (octave_idx_type i = 0; i < nc; i++) \
3134 { \
3135 octave_idx_type jj = 0; \
3136 for (octave_idx_type j = cidx (i); j < cidx (i+1); j++) \
3137 { \
3138 if (jj == ridx (j)) \
3139 { \
3140 nel++; \
3141 jj++; \
3142 } \
3143 else \
3144 break; \
3145 } \
3146 } \
3147 retval = RET_TYPE (nr, nc, nel); \
3148 retval.cidx (0) = 0; \
3149 octave_idx_type ii = 0; \
3150 if (nanflag) \
3151 { \
3152 for (octave_idx_type i = 0; i < nc; i++) \
3153 { \
3154 ELT_TYPE t = ELT_TYPE (1.); \
3155 octave_idx_type jj = 0; \
3156 for (octave_idx_type j = cidx (i); j < cidx (i+1); j++) \
3157 { \
3158 if (jj == ridx (j)) \
3159 { \
3160 if (! octave::math::isnan (data (j))) \
3161 t *= data (j); \
3162 retval.data (ii) = t; \
3163 retval.ridx (ii++) = jj++; \
3164 } \
3165 else \
3166 break; \
3167 } \
3168 retval.cidx (i+1) = ii; \
3169 } \
3170 } \
3171 else \
3172 { \
3173 for (octave_idx_type i = 0; i < nc; i++) \
3174 { \
3175 ELT_TYPE t = ELT_TYPE (1.); \
3176 octave_idx_type jj = 0; \
3177 for (octave_idx_type j = cidx (i); j < cidx (i+1); j++) \
3178 { \
3179 if (jj == ridx (j)) \
3180 { \
3181 t *= data (j); \
3182 retval.data (ii) = t; \
3183 retval.ridx (ii++) = jj++; \
3184 } \
3185 else \
3186 break; \
3187 } \
3188 retval.cidx (i+1) = ii; \
3189 } \
3190 } \
3191 } \
3192 } \
3193 else \
3194 retval = RET_TYPE (nr, nc); \
3195 \
3196 return retval
3197
3198// Reduction operations for sparse matrices ('prod', 'sum', and 'sumsq')
3199
3200#define SPARSE_SUMSQ_HEADER(RET_TYPE, EXPR) \
3201 if (rows () > 0 && cols () > 0 && dim > 1) \
3202 { \
3203 RET_TYPE r = RET_TYPE (rows (), cols (), nnz ()); \
3204 r.cidx (0) = static_cast<octave_idx_type> (0); \
3205 octave_idx_type nel = 0; \
3206 for (octave_idx_type j = 0; j < cols (); j++) \
3207 { \
3208 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++) \
3209 if (data (i) != 0.0) \
3210 { \
3211 r.ridx (nel) = ridx (i); \
3212 EXPR \
3213 } \
3214 r.cidx (j + 1) = nel; \
3215 } \
3216 r.maybe_compress (false); \
3217 return r; \
3218 } \
3219 if (dim > 1) \
3220 return *this;
3221
3222#define SPARSE_BASE_REDUCTION_OP(RET_TYPE, EL_TYPE, ROW_EXPR, COL_EXPR, \
3223 INIT_VAL, MT_RESULT) \
3224 \
3225 octave_idx_type nr = rows (); \
3226 octave_idx_type nc = cols (); \
3227 \
3228 RET_TYPE retval; \
3229 \
3230 if (nr > 0 && nc > 0) \
3231 { \
3232 if ((nr == 1 && dim == -1) || dim == 1) \
3233 { \
3234 /* Define j here to allow fancy definition for prod method */ \
3235 octave_idx_type j = 0; \
3236 OCTAVE_LOCAL_BUFFER (EL_TYPE, tmp, nr); \
3237 \
3238 for (octave_idx_type i = 0; i < nr; i++) \
3239 tmp[i] = INIT_VAL; \
3240 for (j = 0; j < nc; j++) \
3241 { \
3242 for (octave_idx_type i = cidx (j); i < cidx (j + 1); i++) \
3243 { \
3244 ROW_EXPR; \
3245 } \
3246 } \
3247 octave_idx_type nel = 0; \
3248 for (octave_idx_type i = 0; i < nr; i++) \
3249 if (tmp[i] != EL_TYPE ()) \
3250 nel++; \
3251 retval = RET_TYPE (nr, static_cast<octave_idx_type> (1), nel); \
3252 retval.cidx (0) = 0; \
3253 retval.cidx (1) = nel; \
3254 nel = 0; \
3255 for (octave_idx_type i = 0; i < nr; i++) \
3256 if (tmp[i] != EL_TYPE ()) \
3257 { \
3258 retval.data (nel) = tmp[i]; \
3259 retval.ridx (nel++) = i; \
3260 } \
3261 } \
3262 else \
3263 { \
3264 OCTAVE_LOCAL_BUFFER (EL_TYPE, tmp, nc); \
3265 \
3266 for (octave_idx_type j = 0; j < nc; j++) \
3267 { \
3268 tmp[j] = INIT_VAL; \
3269 for (octave_idx_type i = cidx (j); i < cidx (j + 1); i++) \
3270 { \
3271 COL_EXPR; \
3272 } \
3273 } \
3274 octave_idx_type nel = 0; \
3275 for (octave_idx_type i = 0; i < nc; i++) \
3276 if (tmp[i] != EL_TYPE ()) \
3277 nel++; \
3278 retval = RET_TYPE (static_cast<octave_idx_type> (1), nc, nel); \
3279 retval.cidx (0) = 0; \
3280 nel = 0; \
3281 for (octave_idx_type i = 0; i < nc; i++) \
3282 if (tmp[i] != EL_TYPE ()) \
3283 { \
3284 retval.data (nel) = tmp[i]; \
3285 retval.ridx (nel++) = 0; \
3286 retval.cidx (i+1) = retval.cidx (i) + 1; \
3287 } \
3288 else \
3289 retval.cidx (i+1) = retval.cidx (i); \
3290 } \
3291 } \
3292 else if (nc == 0 && nr == 0 && dim == -1) \
3293 { \
3294 if (MT_RESULT) \
3295 { \
3296 retval = RET_TYPE (static_cast<octave_idx_type> (1), \
3297 static_cast<octave_idx_type> (1), \
3298 static_cast<octave_idx_type> (1)); \
3299 retval.cidx (0) = 0; \
3300 retval.cidx (1) = 1; \
3301 retval.ridx (0) = 0; \
3302 retval.data (0) = MT_RESULT; \
3303 } \
3304 else \
3305 retval = RET_TYPE (static_cast<octave_idx_type> (1), \
3306 static_cast<octave_idx_type> (1), \
3307 static_cast<octave_idx_type> (0)); \
3308 } \
3309 else if (nc == 0 && nr == 0 && dim == 0) \
3310 retval = RET_TYPE (static_cast<octave_idx_type> (1), \
3311 static_cast<octave_idx_type> (0), \
3312 static_cast<octave_idx_type> (0)); \
3313 else if (nc == 0 && nr == 0 && dim == 1) \
3314 retval = RET_TYPE (static_cast<octave_idx_type> (0), \
3315 static_cast<octave_idx_type> (1), \
3316 static_cast<octave_idx_type> (0)); \
3317 else if (nr == 0 && (dim == 0 || dim == -1)) \
3318 { \
3319 if (MT_RESULT) \
3320 { \
3321 retval = RET_TYPE (static_cast<octave_idx_type> (1), nc, nc); \
3322 retval.cidx (0) = 0; \
3323 for (octave_idx_type i = 0; i < nc ; i++) \
3324 { \
3325 retval.ridx (i) = 0; \
3326 retval.cidx (i+1) = i+1; \
3327 retval.data (i) = MT_RESULT; \
3328 } \
3329 } \
3330 else \
3331 retval = RET_TYPE (static_cast<octave_idx_type> (1), nc, \
3332 static_cast<octave_idx_type> (0)); \
3333 } \
3334 else if (nc == 0 && (dim == 1 || dim == -1)) \
3335 { \
3336 if (MT_RESULT) \
3337 { \
3338 retval = RET_TYPE (nr, static_cast<octave_idx_type> (1), nr); \
3339 retval.cidx (0) = 0; \
3340 retval.cidx (1) = nr; \
3341 for (octave_idx_type i = 0; i < nr; i++) \
3342 { \
3343 retval.ridx (i) = i; \
3344 retval.data (i) = MT_RESULT; \
3345 } \
3346 } \
3347 else \
3348 retval = RET_TYPE (nr, static_cast<octave_idx_type> (1), \
3349 static_cast<octave_idx_type> (0)); \
3350 } \
3351 else \
3352 retval.resize (nr > 0, nc > 0); \
3353 \
3354 return retval
3355
3356// High accuracy accumulator for sparse matrices ('xsum')
3357
3358#define SPARSE_XSUM_REDUCTION_OP(RET_TYPE, EL_TYPE) \
3359 \
3360 octave_idx_type nr = rows (); \
3361 octave_idx_type nc = cols (); \
3362 \
3363 EL_TYPE inf = std::numeric_limits<EL_TYPE>::infinity (); \
3364 \
3365 RET_TYPE retval; \
3366 \
3367 if (nr > 0 && nc > 0) \
3368 { \
3369 if ((nr == 1 && dim == -1) || dim == 1) \
3370 { \
3371 OCTAVE_LOCAL_BUFFER (EL_TYPE, tmp, nr); \
3372 \
3373 for (octave_idx_type i = 0; i < nr; i++) \
3374 tmp[i] = 0.0; \
3375 for (octave_idx_type i = 0; i < nr; i++) \
3376 { \
3377 EL_TYPE acc = 0.0; \
3378 EL_TYPE err = 0.0; \
3379 bool posinf = false; \
3380 bool neginf = false; \
3381 for (octave_idx_type j = 0; j < nc; j++) \
3382 { \
3383 EL_TYPE d = elem (i, j); \
3384 if (d == EL_TYPE ()); \
3385 else if (nanflag && octave::math::isnan (d)); \
3386 else if (! octave::math::isinf (d)) \
3387 twosum_accum (acc, err, d); \
3388 else if (d > 0.0) \
3389 posinf = true; \
3390 else \
3391 neginf = true; \
3392 } \
3393 if (posinf && neginf) \
3394 tmp[i] = NAN; \
3395 else if (posinf) \
3396 tmp[i] = acc + err + inf; \
3397 else if (neginf) \
3398 tmp[i] = acc + err - inf; \
3399 else \
3400 tmp[i] = acc + err; \
3401 } \
3402 octave_idx_type nel = 0; \
3403 for (octave_idx_type i = 0; i < nr; i++) \
3404 if (tmp[i] != EL_TYPE ()) \
3405 nel++; \
3406 retval = RET_TYPE (nr, static_cast<octave_idx_type> (1), nel); \
3407 retval.cidx (0) = 0; \
3408 retval.cidx (1) = nel; \
3409 nel = 0; \
3410 for (octave_idx_type i = 0; i < nr; i++) \
3411 if (tmp[i] != EL_TYPE ()) \
3412 { \
3413 retval.data (nel) = tmp[i]; \
3414 retval.ridx (nel++) = i; \
3415 } \
3416 } \
3417 else \
3418 { \
3419 OCTAVE_LOCAL_BUFFER (EL_TYPE, tmp, nc); \
3420 \
3421 for (octave_idx_type j = 0; j < nc; j++) \
3422 { \
3423 EL_TYPE acc = 0.0; \
3424 EL_TYPE err = 0.0; \
3425 bool posinf = false; \
3426 bool neginf = false; \
3427 for (octave_idx_type i = 0; i < nr; i++) \
3428 { \
3429 EL_TYPE d = elem (i, j); \
3430 if (d == EL_TYPE ()); \
3431 else if (nanflag && octave::math::isnan (d)); \
3432 else if (! octave::math::isinf (d)) \
3433 twosum_accum (acc, err, d); \
3434 else if (d > 0.0) \
3435 posinf = true; \
3436 else \
3437 neginf = true; \
3438 } \
3439 if (posinf && neginf) \
3440 tmp[j] = NAN; \
3441 else if (posinf) \
3442 tmp[j] = acc + err + inf; \
3443 else if (neginf) \
3444 tmp[j] = acc + err - inf; \
3445 else \
3446 tmp[j] = acc + err; \
3447 } \
3448 octave_idx_type nel = 0; \
3449 for (octave_idx_type i = 0; i < nc; i++) \
3450 if (tmp[i] != EL_TYPE ()) \
3451 nel++; \
3452 retval = RET_TYPE (static_cast<octave_idx_type> (1), nc, nel); \
3453 retval.cidx (0) = 0; \
3454 nel = 0; \
3455 for (octave_idx_type i = 0; i < nc; i++) \
3456 if (tmp[i] != EL_TYPE ()) \
3457 { \
3458 retval.data (nel) = tmp[i]; \
3459 retval.ridx (nel++) = 0; \
3460 retval.cidx (i+1) = retval.cidx (i) + 1; \
3461 } \
3462 else \
3463 retval.cidx (i+1) = retval.cidx (i); \
3464 } \
3465 } \
3466 else if (nc == 0 && nr == 0 && dim == -1) \
3467 retval = RET_TYPE (static_cast<octave_idx_type> (1), \
3468 static_cast<octave_idx_type> (1), \
3469 static_cast<octave_idx_type> (0)); \
3470 else if (nc == 0 && nr == 0 && dim == 0) \
3471 retval = RET_TYPE (static_cast<octave_idx_type> (1), \
3472 static_cast<octave_idx_type> (0), \
3473 static_cast<octave_idx_type> (0)); \
3474 else if (nc == 0 && nr == 0 && dim == 1) \
3475 retval = RET_TYPE (static_cast<octave_idx_type> (0), \
3476 static_cast<octave_idx_type> (1), \
3477 static_cast<octave_idx_type> (0)); \
3478 else if (nr == 0 && (dim == 0 || dim == -1)) \
3479 retval = RET_TYPE (static_cast<octave_idx_type> (1), nc, \
3480 static_cast<octave_idx_type> (0)); \
3481 else if (nc == 0 && (dim == 1 || dim == -1)) \
3482 retval = RET_TYPE (nr, static_cast<octave_idx_type> (1), \
3483 static_cast<octave_idx_type> (0)); \
3484 else \
3485 retval.resize (nr > 0, nc > 0); \
3486 \
3487 return retval
3488
3489// Specializations for 'any' and 'all' functions
3490
3491#define SPARSE_ANY_ALL_HEADER \
3492 if (rows () > 0 && cols () > 0 && dim > 1) \
3493 { \
3494 SparseBoolMatrix r = SparseBoolMatrix (rows (), cols (), nnz ()); \
3495 r.cidx (0) = static_cast<octave_idx_type> (0); \
3496 octave_idx_type nel = 0; \
3497 for (octave_idx_type j = 0; j < cols (); j++) \
3498 { \
3499 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++) \
3500 if (data (i) != 0.0) \
3501 { \
3502 r.ridx (nel) = ridx (i); \
3503 r.data (nel++) = true; \
3504 } \
3505 r.cidx (j + 1) = nel; \
3506 } \
3507 r.maybe_compress (false); \
3508 return r; \
3509 }
3510
3511// Don't break from this loop if the test succeeds because
3512// we are looping over the rows and not the columns in the inner loop.
3513#define SPARSE_ANY_ALL_OP_ROW_CODE(TEST_OP, TEST_TRUE_VAL) \
3514 if (data (i) TEST_OP 0.0) \
3515 tmp[ridx (i)] = TEST_TRUE_VAL;
3516
3517#define SPARSE_ANY_ALL_OP_COL_CODE(TEST_OP, TEST_TRUE_VAL) \
3518 if (data (i) TEST_OP 0.0) \
3519 { \
3520 tmp[j] = TEST_TRUE_VAL; \
3521 break; \
3522 }
3523
3524#define SPARSE_ANY_ALL_OP(DIM, INIT_VAL, MT_RESULT, TEST_OP, TEST_TRUE_VAL) \
3525 SPARSE_BASE_REDUCTION_OP (SparseBoolMatrix, char, \
3526 SPARSE_ANY_ALL_OP_ROW_CODE (TEST_OP, TEST_TRUE_VAL), \
3527 SPARSE_ANY_ALL_OP_COL_CODE (TEST_OP, TEST_TRUE_VAL), \
3528 INIT_VAL, MT_RESULT)
3529
3530#define SPARSE_ALL_OP(DIM) \
3531 if ((rows () == 1 && dim == -1) || dim == 1) \
3532 return transpose (). all (0). transpose (); \
3533 else \
3534 { \
3535 SPARSE_ANY_ALL_OP (DIM, (cidx (j+1) - cidx (j) < nr ? false : true), \
3536 true, ==, false); \
3537 }
3538
3539#define SPARSE_ANY_OP(DIM) SPARSE_ANY_ALL_OP (DIM, false, false, !=, true)
3540
3541#define SPARSE_SPARSE_MUL(RET_TYPE, RET_EL_TYPE, EL_TYPE) \
3542 octave_idx_type nr = m.rows (); \
3543 octave_idx_type nc = m.cols (); \
3544 \
3545 octave_idx_type a_nr = a.rows (); \
3546 octave_idx_type a_nc = a.cols (); \
3547 \
3548 if (nr == 1 && nc == 1) \
3549 { \
3550 RET_EL_TYPE s = m.elem (0,0); \
3551 if (octave::math::isnan (s)) \
3552 { \
3553 RET_TYPE r (a_nr, a_nc, a_nr * a_nc); \
3554 for (octave_idx_type i = 0 ; i < r.numel () ; i++) \
3555 r.elem(i) = octave::numeric_limits<EL_TYPE>::NaN (); \
3556 return r; \
3557 } \
3558 if (octave::math::isinf (s)) \
3559 { \
3560 RET_TYPE r (a_nr, a_nc, a_nr * a_nc); \
3561 for (octave_idx_type j = 0 ; j < a_nc ; j++) \
3562 { \
3563 octave_quit (); \
3564 for (octave_idx_type i = 0 ; i < a_nr ; i++) \
3565 { \
3566 if (a.elem (i, j) == 0.0) \
3567 r.elem (i, j) = octave::numeric_limits<EL_TYPE>::NaN (); \
3568 else \
3569 r.elem (i, j) = s * a.elem (i, j); \
3570 } \
3571 } \
3572 r.maybe_compress (true); \
3573 return r; \
3574 } \
3575 octave_idx_type nz = a.nnz (); \
3576 RET_TYPE r (a_nr, a_nc, nz); \
3577 \
3578 for (octave_idx_type i = 0; i < nz; i++) \
3579 { \
3580 octave_quit (); \
3581 r.data (i) = s * a.data (i); \
3582 r.ridx (i) = a.ridx (i); \
3583 } \
3584 for (octave_idx_type i = 0; i < a_nc + 1; i++) \
3585 { \
3586 octave_quit (); \
3587 r.cidx (i) = a.cidx (i); \
3588 } \
3589 \
3590 r.maybe_compress (true); \
3591 return r; \
3592 } \
3593 else if (a_nr == 1 && a_nc == 1) \
3594 { \
3595 RET_EL_TYPE s = a.elem (0,0); \
3596 if (octave::math::isnan (s)) \
3597 { \
3598 RET_TYPE r (nr, nc, nr * nc); \
3599 for (octave_idx_type i = 0 ; i < r.numel () ; i++) \
3600 r.elem(i) = octave::numeric_limits<EL_TYPE>::NaN (); \
3601 return r; \
3602 } \
3603 if (octave::math::isinf (s)) \
3604 { \
3605 RET_TYPE r (nr, nc, nr * nc); \
3606 for (octave_idx_type j = 0 ; j < nc ; j++) \
3607 { \
3608 octave_quit (); \
3609 for (octave_idx_type i = 0 ; i < nr ; i++) \
3610 { \
3611 if (m.elem (i, j) == 0.0) \
3612 r.elem (i, j) = octave::numeric_limits<EL_TYPE>::NaN (); \
3613 else \
3614 r.elem (i, j) = s * m.elem (i, j); \
3615 } \
3616 } \
3617 r.maybe_compress (true); \
3618 return r; \
3619 } \
3620 octave_idx_type nz = m.nnz (); \
3621 RET_TYPE r (nr, nc, nz); \
3622 \
3623 for (octave_idx_type i = 0; i < nz; i++) \
3624 { \
3625 octave_quit (); \
3626 r.data (i) = m.data (i) * s; \
3627 r.ridx (i) = m.ridx (i); \
3628 } \
3629 for (octave_idx_type i = 0; i < nc + 1; i++) \
3630 { \
3631 octave_quit (); \
3632 r.cidx (i) = m.cidx (i); \
3633 } \
3634 \
3635 r.maybe_compress (true); \
3636 return r; \
3637 } \
3638 else if (nc != a_nr) \
3639 octave::err_nonconformant ("operator *", nr, nc, a_nr, a_nc); \
3640 else \
3641 { \
3642 OCTAVE_LOCAL_BUFFER (octave_idx_type, w, nr); \
3643 RET_TYPE retval (nr, a_nc, static_cast<octave_idx_type> (0)); \
3644 for (octave_idx_type i = 0; i < nr; i++) \
3645 w[i] = 0; \
3646 retval.xcidx (0) = 0; \
3647 \
3648 octave_idx_type nel = 0; \
3649 \
3650 for (octave_idx_type i = 0; i < a_nc; i++) \
3651 { \
3652 for (octave_idx_type j = a.cidx (i); j < a.cidx (i+1); j++) \
3653 { \
3654 octave_idx_type col = a.ridx (j); \
3655 for (octave_idx_type k = m.cidx (col) ; k < m.cidx (col+1); k++) \
3656 { \
3657 if (w[m.ridx (k)] < i + 1) \
3658 { \
3659 w[m.ridx (k)] = i + 1; \
3660 nel++; \
3661 } \
3662 octave_quit (); \
3663 } \
3664 } \
3665 retval.xcidx (i+1) = nel; \
3666 } \
3667 \
3668 if (nel == 0) \
3669 return RET_TYPE (nr, a_nc); \
3670 else \
3671 { \
3672 for (octave_idx_type i = 0; i < nr; i++) \
3673 w[i] = 0; \
3674 \
3675 OCTAVE_LOCAL_BUFFER (RET_EL_TYPE, Xcol, nr); \
3676 \
3677 retval.change_capacity (nel); \
3678 /* The optimal break-point as estimated from simulations */ \
3679 /* Note that Mergesort is O(nz log(nz)) while searching */ \
3680 /* all values is O(nr), where nz here is nonzero per row */ \
3681 /* of length nr. The test itself was then derived from */ \
3682 /* the simulation with random square matrices and the */ \
3683 /* observation of the number of nonzero elements in the */ \
3684 /* output matrix it was found that the breakpoints were */ \
3685 /* nr: 500 1000 2000 5000 10000 */ \
3686 /* nz: 6 25 97 585 2202 */ \
3687 /* The below is a simplication of the 'polyfit'-ed */ \
3688 /* parameters to these breakpoints */ \
3689 octave_idx_type n_per_col = (a_nc > 43000 ? 43000 : \
3690 (a_nc * a_nc) / 43000); \
3691 octave_idx_type ii = 0; \
3692 octave_idx_type *ri = retval.xridx (); \
3693 octave_sort<octave_idx_type> sort; \
3694 \
3695 for (octave_idx_type i = 0; i < a_nc ; i++) \
3696 { \
3697 if (retval.xcidx (i+1) - retval.xcidx (i) > n_per_col) \
3698 { \
3699 for (octave_idx_type j = a.cidx (i); j < a.cidx (i+1); j++) \
3700 { \
3701 octave_idx_type col = a.ridx (j); \
3702 EL_TYPE tmpval = a.data (j); \
3703 for (octave_idx_type k = m.cidx (col) ; \
3704 k < m.cidx (col+1); k++) \
3705 { \
3706 octave_quit (); \
3707 octave_idx_type row = m.ridx (k); \
3708 if (w[row] < i + 1) \
3709 { \
3710 w[row] = i + 1; \
3711 Xcol[row] = tmpval * m.data (k); \
3712 } \
3713 else \
3714 Xcol[row] += tmpval * m.data (k); \
3715 } \
3716 } \
3717 for (octave_idx_type k = 0; k < nr; k++) \
3718 if (w[k] == i + 1) \
3719 { \
3720 retval.xdata (ii) = Xcol[k]; \
3721 retval.xridx (ii++) = k; \
3722 } \
3723 } \
3724 else \
3725 { \
3726 for (octave_idx_type j = a.cidx (i); j < a.cidx (i+1); j++) \
3727 { \
3728 octave_idx_type col = a.ridx (j); \
3729 EL_TYPE tmpval = a.data (j); \
3730 for (octave_idx_type k = m.cidx (col) ; \
3731 k < m.cidx (col+1); k++) \
3732 { \
3733 octave_quit (); \
3734 octave_idx_type row = m.ridx (k); \
3735 if (w[row] < i + 1) \
3736 { \
3737 w[row] = i + 1; \
3738 retval.xridx (ii++) = row; \
3739 Xcol[row] = tmpval * m.data (k); \
3740 } \
3741 else \
3742 Xcol[row] += tmpval * m.data (k); \
3743 } \
3744 } \
3745 sort.sort (ri + retval.xcidx (i), ii - retval.xcidx (i)); \
3746 for (octave_idx_type k = retval.xcidx (i); k < ii; k++) \
3747 retval.xdata (k) = Xcol[retval.xridx (k)]; \
3748 } \
3749 } \
3750 retval.maybe_compress (true); \
3751 return retval; \
3752 } \
3753 }
3754
3755#define SPARSE_FULL_MUL(RET_TYPE, EL_TYPE) \
3756 octave_idx_type nr = m.rows (); \
3757 octave_idx_type nc = m.cols (); \
3758 \
3759 octave_idx_type a_nr = a.rows (); \
3760 octave_idx_type a_nc = a.cols (); \
3761 \
3762 if (nr == 1 && nc == 1) \
3763 { \
3764 RET_TYPE retval = m.elem (0,0) * a; \
3765 return retval; \
3766 } \
3767 else if (nc != a_nr) \
3768 octave::err_nonconformant ("operator *", nr, nc, a_nr, a_nc); \
3769 else \
3770 { \
3771 RET_TYPE::element_type zero = RET_TYPE::element_type (); \
3772 \
3773 RET_TYPE retval (nr, a_nc, zero); \
3774 \
3775 for (octave_idx_type i = 0; i < a_nc ; i++) \
3776 { \
3777 for (octave_idx_type j = 0; j < a_nr; j++) \
3778 { \
3779 octave_quit (); \
3780 \
3781 EL_TYPE tmpval = a.elem (j,i); \
3782 for (octave_idx_type k = m.cidx (j) ; k < m.cidx (j+1); k++) \
3783 retval.elem (m.ridx (k),i) += tmpval * m.data (k); \
3784 } \
3785 } \
3786 return retval; \
3787 }
3788
3789#define SPARSE_FULL_TRANS_MUL(RET_TYPE, EL_TYPE, CONJ_OP) \
3790 octave_idx_type nr = m.rows (); \
3791 octave_idx_type nc = m.cols (); \
3792 \
3793 octave_idx_type a_nr = a.rows (); \
3794 octave_idx_type a_nc = a.cols (); \
3795 \
3796 if (nr == 1 && nc == 1) \
3797 { \
3798 RET_TYPE retval = CONJ_OP (m.elem (0,0)) * a; \
3799 return retval; \
3800 } \
3801 else if (nr != a_nr) \
3802 octave::err_nonconformant ("operator *", nc, nr, a_nr, a_nc); \
3803 else \
3804 { \
3805 RET_TYPE retval (nc, a_nc); \
3806 \
3807 for (octave_idx_type i = 0; i < a_nc ; i++) \
3808 { \
3809 for (octave_idx_type j = 0; j < nc; j++) \
3810 { \
3811 octave_quit (); \
3812 \
3813 EL_TYPE acc = EL_TYPE (); \
3814 for (octave_idx_type k = m.cidx (j) ; k < m.cidx (j+1); k++) \
3815 acc += a.elem (m.ridx (k),i) * CONJ_OP (m.data (k)); \
3816 retval.xelem (j,i) = acc; \
3817 } \
3818 } \
3819 return retval; \
3820 }
3821
3822#define FULL_SPARSE_MUL(RET_TYPE, EL_TYPE) \
3823 octave_idx_type nr = m.rows (); \
3824 octave_idx_type nc = m.cols (); \
3825 \
3826 octave_idx_type a_nr = a.rows (); \
3827 octave_idx_type a_nc = a.cols (); \
3828 \
3829 if (a_nr == 1 && a_nc == 1) \
3830 { \
3831 RET_TYPE retval = m * a.elem (0,0); \
3832 return retval; \
3833 } \
3834 else if (nc != a_nr) \
3835 octave::err_nonconformant ("operator *", nr, nc, a_nr, a_nc); \
3836 else \
3837 { \
3838 RET_TYPE::element_type zero = RET_TYPE::element_type (); \
3839 \
3840 RET_TYPE retval (nr, a_nc, zero); \
3841 \
3842 for (octave_idx_type i = 0; i < a_nc ; i++) \
3843 { \
3844 octave_quit (); \
3845 for (octave_idx_type j = a.cidx (i); j < a.cidx (i+1); j++) \
3846 { \
3847 octave_idx_type col = a.ridx (j); \
3848 EL_TYPE tmpval = a.data (j); \
3849 \
3850 for (octave_idx_type k = 0 ; k < nr; k++) \
3851 retval.xelem (k,i) += tmpval * m.elem (k,col); \
3852 } \
3853 } \
3854 return retval; \
3855 }
3856
3857#define FULL_SPARSE_MUL_TRANS(RET_TYPE, EL_TYPE, CONJ_OP) \
3858 octave_idx_type nr = m.rows (); \
3859 octave_idx_type nc = m.cols (); \
3860 \
3861 octave_idx_type a_nr = a.rows (); \
3862 octave_idx_type a_nc = a.cols (); \
3863 \
3864 if (a_nr == 1 && a_nc == 1) \
3865 { \
3866 RET_TYPE retval = m * CONJ_OP (a.elem (0,0)); \
3867 return retval; \
3868 } \
3869 else if (nc != a_nc) \
3870 octave::err_nonconformant ("operator *", nr, nc, a_nc, a_nr); \
3871 else \
3872 { \
3873 RET_TYPE::element_type zero = RET_TYPE::element_type (); \
3874 \
3875 RET_TYPE retval (nr, a_nr, zero); \
3876 \
3877 for (octave_idx_type i = 0; i < a_nc ; i++) \
3878 { \
3879 octave_quit (); \
3880 for (octave_idx_type j = a.cidx (i); j < a.cidx (i+1); j++) \
3881 { \
3882 octave_idx_type col = a.ridx (j); \
3883 EL_TYPE tmpval = CONJ_OP (a.data (j)); \
3884 for (octave_idx_type k = 0 ; k < nr; k++) \
3885 retval.xelem (k,col) += tmpval * m.elem (k,i); \
3886 } \
3887 } \
3888 return retval; \
3889 }
3890
3891#endif