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
lo-array-errwarn.h
oct-locbuf.h
Array-util.h
mx-inlines.cc
liboctave
operators
Sparse-op-defs.h
Generated on Sat Sep 1 2018 14:01:46 for GNU Octave by
1.8.14