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