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