Mercurial > octave

comparison liboctave/Sparse-op-defs.h @ 5275:23b37da9fd5b

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2005-04-08 16:07:35 by jwe]
author jwe
date Fri, 08 Apr 2005 16:07:37 +0000
parents 57077d0ddc8e
children 4c8a2e4e0717
comparison
equal deleted inserted replaced
5274:eae7b40388e9 5275:23b37da9fd5b
43 43
44 #define SPARSE_SMS_BIN_OP_1(R, F, OP, M, S) \ 44 #define SPARSE_SMS_BIN_OP_1(R, F, OP, M, S) \
45 R \ 45 R \
46 F (const M& m, const S& s) \ 46 F (const M& m, const S& s) \
47 { \ 47 { \
48 int nr = m.rows (); \ 48 octave_idx_type nr = m.rows (); \
49 int nc = m.cols (); \ 49 octave_idx_type nc = m.cols (); \
50 \ 50 \
51 R r (nr, nc, (0.0 OP s)); \ 51 R r (nr, nc, (0.0 OP s)); \
52 \ 52 \
53 for (int j = 0; j < nc; j++) \ 53 for (octave_idx_type j = 0; j < nc; j++) \
54 for (int i = m.cidx (j); i < m.cidx (j+1); i++) \ 54 for (octave_idx_type i = m.cidx (j); i < m.cidx (j+1); i++) \
55 r.elem (m.ridx (i), j) = m.data (i) OP s; \ 55 r.elem (m.ridx (i), j) = m.data (i) OP s; \
56 return r; \ 56 return r; \
57 } 57 }
58 58
59 #define SPARSE_SMS_BIN_OP_2(R, F, OP, M, S) \ 59 #define SPARSE_SMS_BIN_OP_2(R, F, OP, M, S) \
60 R \ 60 R \
61 F (const M& m, const S& s) \ 61 F (const M& m, const S& s) \
62 { \ 62 { \
63 int nr = m.rows (); \ 63 octave_idx_type nr = m.rows (); \
64 int nc = m.cols (); \ 64 octave_idx_type nc = m.cols (); \
65 int nz = m.nnz (); \ 65 octave_idx_type nz = m.nnz (); \
66 \ 66 \
67 R r (nr, nc, nz); \ 67 R r (nr, nc, nz); \
68 \ 68 \
69 for (int i = 0; i < nz; i++) \ 69 for (octave_idx_type i = 0; i < nz; i++) \
70 { \ 70 { \
71 r.data(i) = m.data(i) OP s; \ 71 r.data(i) = m.data(i) OP s; \
72 r.ridx(i) = m.ridx(i); \ 72 r.ridx(i) = m.ridx(i); \
73 } \ 73 } \
74 for (int i = 0; i < nc + 1; i++) \ 74 for (octave_idx_type i = 0; i < nc + 1; i++) \
75 r.cidx(i) = m.cidx(i); \ 75 r.cidx(i) = m.cidx(i); \
76 \ 76 \
77 r.maybe_compress (true); \ 77 r.maybe_compress (true); \
78 return r; \ 78 return r; \
79 } 79 }
99 #define SPARSE_SMS_CMP_OP(F, OP, M, MZ, MC, S, SZ, SC) \ 99 #define SPARSE_SMS_CMP_OP(F, OP, M, MZ, MC, S, SZ, SC) \
100 SparseBoolMatrix \ 100 SparseBoolMatrix \
101 F (const M& m, const S& s) \ 101 F (const M& m, const S& s) \
102 { \ 102 { \
103 /* Count num of non-zero elements */ \ 103 /* Count num of non-zero elements */ \
104 int nel = 0; \ 104 octave_idx_type nel = 0; \
105 int nz = m.nnz (); \ 105 octave_idx_type nz = m.nnz (); \
106 if (MC (MZ) OP SC (s)) \ 106 if (MC (MZ) OP SC (s)) \
107 nel += m.numel() - nz; \ 107 nel += m.numel() - nz; \
108 for (int i = 0; i < nz; i++) \ 108 for (octave_idx_type i = 0; i < nz; i++) \
109 if (MC (m.data (i)) OP SC (s)) \ 109 if (MC (m.data (i)) OP SC (s)) \
110 nel++; \ 110 nel++; \
111 \ 111 \
112 int nr = m.rows (); \ 112 octave_idx_type nr = m.rows (); \
113 int nc = m.cols (); \ 113 octave_idx_type nc = m.cols (); \
114 SparseBoolMatrix r (nr, nc, nel); \ 114 SparseBoolMatrix r (nr, nc, nel); \
115 \ 115 \
116 if (nr > 0 && nc > 0) \ 116 if (nr > 0 && nc > 0) \
117 { \ 117 { \
118 if (MC (MZ) OP SC (s)) \ 118 if (MC (MZ) OP SC (s)) \
119 { \ 119 { \
120 int ii = 0; \ 120 octave_idx_type ii = 0; \
121 r.cidx (0) = 0; \ 121 r.cidx (0) = 0; \
122 for (int j = 0; j < nc; j++) \ 122 for (octave_idx_type j = 0; j < nc; j++) \
123 { \ 123 { \
124 for (int i = 0; i < nr; i++) \ 124 for (octave_idx_type i = 0; i < nr; i++) \
125 { \ 125 { \
126 bool el = MC (m.elem(i, j)) OP SC (s); \ 126 bool el = MC (m.elem(i, j)) OP SC (s); \
127 if (el) \ 127 if (el) \
128 { \ 128 { \
129 r.data(ii) = el; \ 129 r.data(ii) = el; \
133 r.cidx(j+1) = ii; \ 133 r.cidx(j+1) = ii; \
134 } \ 134 } \
135 } \ 135 } \
136 else \ 136 else \
137 { \ 137 { \
138 int ii = 0; \ 138 octave_idx_type ii = 0; \
139 r.cidx (0) = 0; \ 139 r.cidx (0) = 0; \
140 for (int j = 0; j < nc; j++) \ 140 for (octave_idx_type j = 0; j < nc; j++) \
141 { \ 141 { \
142 for (int i = m.cidx(j); i < m.cidx(j+1); i++) \ 142 for (octave_idx_type i = m.cidx(j); i < m.cidx(j+1); i++) \
143 { \ 143 { \
144 bool el = MC (m.data(i)) OP SC (s); \ 144 bool el = MC (m.data(i)) OP SC (s); \
145 if (el) \ 145 if (el) \
146 { \ 146 { \
147 r.data(ii) = el; \ 147 r.data(ii) = el; \
174 #define SPARSE_SMS_BOOL_OP(F, OP, M, S, LHS_ZERO, RHS_ZERO) \ 174 #define SPARSE_SMS_BOOL_OP(F, OP, M, S, LHS_ZERO, RHS_ZERO) \
175 SparseBoolMatrix \ 175 SparseBoolMatrix \
176 F (const M& m, const S& s) \ 176 F (const M& m, const S& s) \
177 { \ 177 { \
178 /* Count num of non-zero elements */ \ 178 /* Count num of non-zero elements */ \
179 int nel = 0; \ 179 octave_idx_type nel = 0; \
180 int nz = m.nnz (); \ 180 octave_idx_type nz = m.nnz (); \
181 if (LHS_ZERO OP (s != RHS_ZERO)) \ 181 if (LHS_ZERO OP (s != RHS_ZERO)) \
182 nel += m.numel() - nz; \ 182 nel += m.numel() - nz; \
183 for (int i = 0; i < nz; i++) \ 183 for (octave_idx_type i = 0; i < nz; i++) \
184 if ((m.data(i) != LHS_ZERO) OP (s != RHS_ZERO))\ 184 if ((m.data(i) != LHS_ZERO) OP (s != RHS_ZERO))\
185 nel++; \ 185 nel++; \
186 \ 186 \
187 int nr = m.rows (); \ 187 octave_idx_type nr = m.rows (); \
188 int nc = m.cols (); \ 188 octave_idx_type nc = m.cols (); \
189 SparseBoolMatrix r (nr, nc, nel); \ 189 SparseBoolMatrix r (nr, nc, nel); \
190 \ 190 \
191 if (nr > 0 && nc > 0) \ 191 if (nr > 0 && nc > 0) \
192 { \ 192 { \
193 if (LHS_ZERO OP (s != RHS_ZERO)) \ 193 if (LHS_ZERO OP (s != RHS_ZERO)) \
194 { \ 194 { \
195 int ii = 0; \ 195 octave_idx_type ii = 0; \
196 r.cidx (0) = 0; \ 196 r.cidx (0) = 0; \
197 for (int j = 0; j < nc; j++) \ 197 for (octave_idx_type j = 0; j < nc; j++) \
198 { \ 198 { \
199 for (int i = 0; i < nr; i++) \ 199 for (octave_idx_type i = 0; i < nr; i++) \
200 { \ 200 { \
201 bool el = (m.elem(i, j) != LHS_ZERO) OP (s != RHS_ZERO); \ 201 bool el = (m.elem(i, j) != LHS_ZERO) OP (s != RHS_ZERO); \
202 if (el) \ 202 if (el) \
203 { \ 203 { \
204 r.data(ii) = el; \ 204 r.data(ii) = el; \
208 r.cidx(j+1) = ii; \ 208 r.cidx(j+1) = ii; \
209 } \ 209 } \
210 } \ 210 } \
211 else \ 211 else \
212 { \ 212 { \
213 int ii = 0; \ 213 octave_idx_type ii = 0; \
214 r.cidx (0) = 0; \ 214 r.cidx (0) = 0; \
215 for (int j = 0; j < nc; j++) \ 215 for (octave_idx_type j = 0; j < nc; j++) \
216 { \ 216 { \
217 for (int i = m.cidx(j); i < m.cidx(j+1); i++) \ 217 for (octave_idx_type i = m.cidx(j); i < m.cidx(j+1); i++) \
218 { \ 218 { \
219 bool el = (m.data(i) != LHS_ZERO) OP (s != RHS_ZERO); \ 219 bool el = (m.data(i) != LHS_ZERO) OP (s != RHS_ZERO); \
220 if (el) \ 220 if (el) \
221 { \ 221 { \
222 r.data(ii) = el; \ 222 r.data(ii) = el; \
252 252
253 #define SPARSE_SSM_BIN_OP_1(R, F, OP, S, M) \ 253 #define SPARSE_SSM_BIN_OP_1(R, F, OP, S, M) \
254 R \ 254 R \
255 F (const S& s, const M& m) \ 255 F (const S& s, const M& m) \
256 { \ 256 { \
257 int nr = m.rows (); \ 257 octave_idx_type nr = m.rows (); \
258 int nc = m.cols (); \ 258 octave_idx_type nc = m.cols (); \
259 \ 259 \
260 R r (nr, nc, (s OP 0.0)); \ 260 R r (nr, nc, (s OP 0.0)); \
261 \ 261 \
262 for (int j = 0; j < nc; j++) \ 262 for (octave_idx_type j = 0; j < nc; j++) \
263 for (int i = m.cidx (j); i < m.cidx (j+1); i++) \ 263 for (octave_idx_type i = m.cidx (j); i < m.cidx (j+1); i++) \
264 r.elem (m.ridx (i), j) = s OP m.data (i); \ 264 r.elem (m.ridx (i), j) = s OP m.data (i); \
265 \ 265 \
266 return r; \ 266 return r; \
267 } 267 }
268 268
269 #define SPARSE_SSM_BIN_OP_2(R, F, OP, S, M) \ 269 #define SPARSE_SSM_BIN_OP_2(R, F, OP, S, M) \
270 R \ 270 R \
271 F (const S& s, const M& m) \ 271 F (const S& s, const M& m) \
272 { \ 272 { \
273 int nr = m.rows (); \ 273 octave_idx_type nr = m.rows (); \
274 int nc = m.cols (); \ 274 octave_idx_type nc = m.cols (); \
275 int nz = m.nnz (); \ 275 octave_idx_type nz = m.nnz (); \
276 \ 276 \
277 R r (nr, nc, nz); \ 277 R r (nr, nc, nz); \
278 \ 278 \
279 for (int i = 0; i < nz; i++) \ 279 for (octave_idx_type i = 0; i < nz; i++) \
280 { \ 280 { \
281 r.data(i) = s OP m.data(i); \ 281 r.data(i) = s OP m.data(i); \
282 r.ridx(i) = m.ridx(i); \ 282 r.ridx(i) = m.ridx(i); \
283 } \ 283 } \
284 for (int i = 0; i < nc + 1; i++) \ 284 for (octave_idx_type i = 0; i < nc + 1; i++) \
285 r.cidx(i) = m.cidx(i); \ 285 r.cidx(i) = m.cidx(i); \
286 \ 286 \
287 r.maybe_compress(true); \ 287 r.maybe_compress(true); \
288 return r; \ 288 return r; \
289 } 289 }
309 #define SPARSE_SSM_CMP_OP(F, OP, S, SZ, SC, M, MZ, MC) \ 309 #define SPARSE_SSM_CMP_OP(F, OP, S, SZ, SC, M, MZ, MC) \
310 SparseBoolMatrix \ 310 SparseBoolMatrix \
311 F (const S& s, const M& m) \ 311 F (const S& s, const M& m) \
312 { \ 312 { \
313 /* Count num of non-zero elements */ \ 313 /* Count num of non-zero elements */ \
314 int nel = 0; \ 314 octave_idx_type nel = 0; \
315 int nz = m.nnz (); \ 315 octave_idx_type nz = m.nnz (); \
316 if (SC (s) OP MC (MZ)) \ 316 if (SC (s) OP MC (MZ)) \
317 nel += m.numel() - nz; \ 317 nel += m.numel() - nz; \
318 for (int i = 0; i < nz; i++) \ 318 for (octave_idx_type i = 0; i < nz; i++) \
319 if (SC (s) OP MC (m.data (i))) \ 319 if (SC (s) OP MC (m.data (i))) \
320 nel++; \ 320 nel++; \
321 \ 321 \
322 int nr = m.rows (); \ 322 octave_idx_type nr = m.rows (); \
323 int nc = m.cols (); \ 323 octave_idx_type nc = m.cols (); \
324 SparseBoolMatrix r (nr, nc, nel); \ 324 SparseBoolMatrix r (nr, nc, nel); \
325 \ 325 \
326 if (nr > 0 && nc > 0) \ 326 if (nr > 0 && nc > 0) \
327 { \ 327 { \
328 if (SC (s) OP MC (MZ))\ 328 if (SC (s) OP MC (MZ))\
329 { \ 329 { \
330 int ii = 0; \ 330 octave_idx_type ii = 0; \
331 r.cidx (0) = 0; \ 331 r.cidx (0) = 0; \
332 for (int j = 0; j < nc; j++) \ 332 for (octave_idx_type j = 0; j < nc; j++) \
333 { \ 333 { \
334 for (int i = 0; i < nr; i++) \ 334 for (octave_idx_type i = 0; i < nr; i++) \
335 { \ 335 { \
336 bool el = SC (s) OP MC (m.elem(i, j)); \ 336 bool el = SC (s) OP MC (m.elem(i, j)); \
337 if (el) \ 337 if (el) \
338 { \ 338 { \
339 r.data(ii) = el; \ 339 r.data(ii) = el; \
343 r.cidx(j+1) = ii; \ 343 r.cidx(j+1) = ii; \
344 } \ 344 } \
345 } \ 345 } \
346 else \ 346 else \
347 { \ 347 { \
348 int ii = 0; \ 348 octave_idx_type ii = 0; \
349 r.cidx (0) = 0; \ 349 r.cidx (0) = 0; \
350 for (int j = 0; j < nc; j++) \ 350 for (octave_idx_type j = 0; j < nc; j++) \
351 { \ 351 { \
352 for (int i = m.cidx(j); i < m.cidx(j+1); i++) \ 352 for (octave_idx_type i = m.cidx(j); i < m.cidx(j+1); i++) \
353 { \ 353 { \
354 bool el = SC (s) OP MC (m.data(i)); \ 354 bool el = SC (s) OP MC (m.data(i)); \
355 if (el) \ 355 if (el) \
356 { \ 356 { \
357 r.data(ii) = el; \ 357 r.data(ii) = el; \
384 #define SPARSE_SSM_BOOL_OP(F, OP, S, M, LHS_ZERO, RHS_ZERO) \ 384 #define SPARSE_SSM_BOOL_OP(F, OP, S, M, LHS_ZERO, RHS_ZERO) \
385 SparseBoolMatrix \ 385 SparseBoolMatrix \
386 F (const S& s, const M& m) \ 386 F (const S& s, const M& m) \
387 { \ 387 { \
388 /* Count num of non-zero elements */ \ 388 /* Count num of non-zero elements */ \
389 int nel = 0; \ 389 octave_idx_type nel = 0; \
390 int nz = m.nnz (); \ 390 octave_idx_type nz = m.nnz (); \
391 if ((s != LHS_ZERO) OP RHS_ZERO) \ 391 if ((s != LHS_ZERO) OP RHS_ZERO) \
392 nel += m.numel() - nz; \ 392 nel += m.numel() - nz; \
393 for (int i = 0; i < nz; i++) \ 393 for (octave_idx_type i = 0; i < nz; i++) \
394 if ((s != LHS_ZERO) OP m.data(i) != RHS_ZERO) \ 394 if ((s != LHS_ZERO) OP m.data(i) != RHS_ZERO) \
395 nel++; \ 395 nel++; \
396 \ 396 \
397 int nr = m.rows (); \ 397 octave_idx_type nr = m.rows (); \
398 int nc = m.cols (); \ 398 octave_idx_type nc = m.cols (); \
399 SparseBoolMatrix r (nr, nc, nel); \ 399 SparseBoolMatrix r (nr, nc, nel); \
400 \ 400 \
401 if (nr > 0 && nc > 0) \ 401 if (nr > 0 && nc > 0) \
402 { \ 402 { \
403 if ((s != LHS_ZERO) OP RHS_ZERO) \ 403 if ((s != LHS_ZERO) OP RHS_ZERO) \
404 { \ 404 { \
405 int ii = 0; \ 405 octave_idx_type ii = 0; \
406 r.cidx (0) = 0; \ 406 r.cidx (0) = 0; \
407 for (int j = 0; j < nc; j++) \ 407 for (octave_idx_type j = 0; j < nc; j++) \
408 { \ 408 { \
409 for (int i = 0; i < nr; i++) \ 409 for (octave_idx_type i = 0; i < nr; i++) \
410 { \ 410 { \
411 bool el = (s != LHS_ZERO) OP (m.elem(i, j) != RHS_ZERO); \ 411 bool el = (s != LHS_ZERO) OP (m.elem(i, j) != RHS_ZERO); \
412 if (el) \ 412 if (el) \
413 { \ 413 { \
414 r.data(ii) = el; \ 414 r.data(ii) = el; \
418 r.cidx(j+1) = ii; \ 418 r.cidx(j+1) = ii; \
419 } \ 419 } \
420 } \ 420 } \
421 else \ 421 else \
422 { \ 422 { \
423 int ii = 0; \ 423 octave_idx_type ii = 0; \
424 r.cidx (0) = 0; \ 424 r.cidx (0) = 0; \
425 for (int j = 0; j < nc; j++) \ 425 for (octave_idx_type j = 0; j < nc; j++) \
426 { \ 426 { \
427 for (int i = m.cidx(j); i < m.cidx(j+1); i++) \ 427 for (octave_idx_type i = m.cidx(j); i < m.cidx(j+1); i++) \
428 { \ 428 { \
429 bool el = (s != LHS_ZERO) OP (m.data(i) != RHS_ZERO); \ 429 bool el = (s != LHS_ZERO) OP (m.data(i) != RHS_ZERO); \
430 if (el) \ 430 if (el) \
431 { \ 431 { \
432 r.data(ii) = el; \ 432 r.data(ii) = el; \
464 R \ 464 R \
465 F (const M1& m1, const M2& m2) \ 465 F (const M1& m1, const M2& m2) \
466 { \ 466 { \
467 R r; \ 467 R r; \
468 \ 468 \
469 int m1_nr = m1.rows (); \ 469 octave_idx_type m1_nr = m1.rows (); \
470 int m1_nc = m1.cols (); \ 470 octave_idx_type m1_nc = m1.cols (); \
471 \ 471 \
472 int m2_nr = m2.rows (); \ 472 octave_idx_type m2_nr = m2.rows (); \
473 int m2_nc = m2.cols (); \ 473 octave_idx_type m2_nc = m2.cols (); \
474 \ 474 \
475 if (m1_nr != m2_nr || m1_nc != m2_nc) \ 475 if (m1_nr != m2_nr || m1_nc != m2_nc) \
476 gripe_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \ 476 gripe_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \
477 else \ 477 else \
478 { \ 478 { \
479 r = R (m1_nr, m1_nc, (m1.nnz () + m2.nnz ())); \ 479 r = R (m1_nr, m1_nc, (m1.nnz () + m2.nnz ())); \
480 \ 480 \
481 int jx = 0; \ 481 octave_idx_type jx = 0; \
482 r.cidx (0) = 0; \ 482 r.cidx (0) = 0; \
483 for (int i = 0 ; i < m1_nc ; i++) \ 483 for (octave_idx_type i = 0 ; i < m1_nc ; i++) \
484 { \ 484 { \
485 int ja = m1.cidx(i); \ 485 octave_idx_type ja = m1.cidx(i); \
486 int ja_max = m1.cidx(i+1); \ 486 octave_idx_type ja_max = m1.cidx(i+1); \
487 bool ja_lt_max= ja < ja_max; \ 487 bool ja_lt_max= ja < ja_max; \
488 \ 488 \
489 int jb = m2.cidx(i); \ 489 octave_idx_type jb = m2.cidx(i); \
490 int jb_max = m2.cidx(i+1); \ 490 octave_idx_type jb_max = m2.cidx(i+1); \
491 bool jb_lt_max = jb < jb_max; \ 491 bool jb_lt_max = jb < jb_max; \
492 \ 492 \
493 while (ja_lt_max || jb_lt_max ) \ 493 while (ja_lt_max || jb_lt_max ) \
494 { \ 494 { \
495 OCTAVE_QUIT; \ 495 OCTAVE_QUIT; \
538 R \ 538 R \
539 F (const M1& m1, const M2& m2) \ 539 F (const M1& m1, const M2& m2) \
540 { \ 540 { \
541 R r; \ 541 R r; \
542 \ 542 \
543 int m1_nr = m1.rows (); \ 543 octave_idx_type m1_nr = m1.rows (); \
544 int m1_nc = m1.cols (); \ 544 octave_idx_type m1_nc = m1.cols (); \
545 \ 545 \
546 int m2_nr = m2.rows (); \ 546 octave_idx_type m2_nr = m2.rows (); \
547 int m2_nc = m2.cols (); \ 547 octave_idx_type m2_nc = m2.cols (); \
548 \ 548 \
549 if (m1_nr != m2_nr || m1_nc != m2_nc) \ 549 if (m1_nr != m2_nr || m1_nc != m2_nc) \
550 gripe_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \ 550 gripe_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \
551 else \ 551 else \
552 { \ 552 { \
553 r = R (m1_nr, m1_nc, (m1.nnz () > m2.nnz () ? m1.nnz () : m2.nnz ())); \ 553 r = R (m1_nr, m1_nc, (m1.nnz () > m2.nnz () ? m1.nnz () : m2.nnz ())); \
554 \ 554 \
555 int jx = 0; \ 555 octave_idx_type jx = 0; \
556 r.cidx (0) = 0; \ 556 r.cidx (0) = 0; \
557 for (int i = 0 ; i < m1_nc ; i++) \ 557 for (octave_idx_type i = 0 ; i < m1_nc ; i++) \
558 { \ 558 { \
559 int ja = m1.cidx(i); \ 559 octave_idx_type ja = m1.cidx(i); \
560 int ja_max = m1.cidx(i+1); \ 560 octave_idx_type ja_max = m1.cidx(i+1); \
561 bool ja_lt_max= ja < ja_max; \ 561 bool ja_lt_max= ja < ja_max; \
562 \ 562 \
563 int jb = m2.cidx(i); \ 563 octave_idx_type jb = m2.cidx(i); \
564 int jb_max = m2.cidx(i+1); \ 564 octave_idx_type jb_max = m2.cidx(i+1); \
565 bool jb_lt_max = jb < jb_max; \ 565 bool jb_lt_max = jb < jb_max; \
566 \ 566 \
567 while (ja_lt_max || jb_lt_max ) \ 567 while (ja_lt_max || jb_lt_max ) \
568 { \ 568 { \
569 OCTAVE_QUIT; \ 569 OCTAVE_QUIT; \
602 R \ 602 R \
603 F (const M1& m1, const M2& m2) \ 603 F (const M1& m1, const M2& m2) \
604 { \ 604 { \
605 R r; \ 605 R r; \
606 \ 606 \
607 int m1_nr = m1.rows (); \ 607 octave_idx_type m1_nr = m1.rows (); \
608 int m1_nc = m1.cols (); \ 608 octave_idx_type m1_nc = m1.cols (); \
609 \ 609 \
610 int m2_nr = m2.rows (); \ 610 octave_idx_type m2_nr = m2.rows (); \
611 int m2_nc = m2.cols (); \ 611 octave_idx_type m2_nc = m2.cols (); \
612 \ 612 \
613 if (m1_nr != m2_nr || m1_nc != m2_nc) \ 613 if (m1_nr != m2_nr || m1_nc != m2_nc) \
614 gripe_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \ 614 gripe_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \
615 else \ 615 else \
616 { \ 616 { \
617 \ 617 \
618 /* XXX FIXME XXX Kludge... Always double/Complex, so Complex () */ \ 618 /* XXX FIXME XXX Kludge... Always double/Complex, so Complex () */ \
619 r = R (m1_nr, m1_nc, (Complex () OP Complex ())); \ 619 r = R (m1_nr, m1_nc, (Complex () OP Complex ())); \
620 \ 620 \
621 for (int i = 0 ; i < m1_nc ; i++) \ 621 for (octave_idx_type i = 0 ; i < m1_nc ; i++) \
622 { \ 622 { \
623 int ja = m1.cidx(i); \ 623 octave_idx_type ja = m1.cidx(i); \
624 int ja_max = m1.cidx(i+1); \ 624 octave_idx_type ja_max = m1.cidx(i+1); \
625 bool ja_lt_max= ja < ja_max; \ 625 bool ja_lt_max= ja < ja_max; \
626 \ 626 \
627 int jb = m2.cidx(i); \ 627 octave_idx_type jb = m2.cidx(i); \
628 int jb_max = m2.cidx(i+1); \ 628 octave_idx_type jb_max = m2.cidx(i+1); \
629 bool jb_lt_max = jb < jb_max; \ 629 bool jb_lt_max = jb < jb_max; \
630 \ 630 \
631 while (ja_lt_max || jb_lt_max ) \ 631 while (ja_lt_max || jb_lt_max ) \
632 { \ 632 { \
633 OCTAVE_QUIT; \ 633 OCTAVE_QUIT; \
689 SparseBoolMatrix \ 689 SparseBoolMatrix \
690 F (const M1& m1, const M2& m2) \ 690 F (const M1& m1, const M2& m2) \
691 { \ 691 { \
692 SparseBoolMatrix r; \ 692 SparseBoolMatrix r; \
693 \ 693 \
694 int m1_nr = m1.rows (); \ 694 octave_idx_type m1_nr = m1.rows (); \
695 int m1_nc = m1.cols (); \ 695 octave_idx_type m1_nc = m1.cols (); \
696 \ 696 \
697 int m2_nr = m2.rows (); \ 697 octave_idx_type m2_nr = m2.rows (); \
698 int m2_nc = m2.cols (); \ 698 octave_idx_type m2_nc = m2.cols (); \
699 \ 699 \
700 if (m1_nr == m2_nr && m1_nc == m2_nc) \ 700 if (m1_nr == m2_nr && m1_nc == m2_nc) \
701 { \ 701 { \
702 if (m1_nr != 0 || m1_nc != 0) \ 702 if (m1_nr != 0 || m1_nc != 0) \
703 { \ 703 { \
704 /* Count num of non-zero elements */ \ 704 /* Count num of non-zero elements */ \
705 int nel = 0; \ 705 octave_idx_type nel = 0; \
706 for (int j = 0; j < m1_nc; j++) \ 706 for (octave_idx_type j = 0; j < m1_nc; j++) \
707 for (int i = 0; i < m1_nr; i++) \ 707 for (octave_idx_type i = 0; i < m1_nr; i++) \
708 if (C1 (m1.elem(i, j)) OP C2 (m2.elem(i, j))) \ 708 if (C1 (m1.elem(i, j)) OP C2 (m2.elem(i, j))) \
709 nel++; \ 709 nel++; \
710 \ 710 \
711 r = SparseBoolMatrix (m1_nr, m1_nc, nel); \ 711 r = SparseBoolMatrix (m1_nr, m1_nc, nel); \
712 \ 712 \
713 int ii = 0; \ 713 octave_idx_type ii = 0; \
714 r.cidx (0) = 0; \ 714 r.cidx (0) = 0; \
715 for (int j = 0; j < m1_nc; j++) \ 715 for (octave_idx_type j = 0; j < m1_nc; j++) \
716 { \ 716 { \
717 for (int i = 0; i < m1_nr; i++) \ 717 for (octave_idx_type i = 0; i < m1_nr; i++) \
718 { \ 718 { \
719 bool el = C1 (m1.elem(i, j)) OP C2 (m2.elem(i, j)); \ 719 bool el = C1 (m1.elem(i, j)) OP C2 (m2.elem(i, j)); \
720 if (el) \ 720 if (el) \
721 { \ 721 { \
722 r.data(ii) = el; \ 722 r.data(ii) = el; \
755 SparseBoolMatrix \ 755 SparseBoolMatrix \
756 F (const M1& m1, const M2& m2) \ 756 F (const M1& m1, const M2& m2) \
757 { \ 757 { \
758 SparseBoolMatrix r; \ 758 SparseBoolMatrix r; \
759 \ 759 \
760 int m1_nr = m1.rows (); \ 760 octave_idx_type m1_nr = m1.rows (); \
761 int m1_nc = m1.cols (); \ 761 octave_idx_type m1_nc = m1.cols (); \
762 \ 762 \
763 int m2_nr = m2.rows (); \ 763 octave_idx_type m2_nr = m2.rows (); \
764 int m2_nc = m2.cols (); \ 764 octave_idx_type m2_nc = m2.cols (); \
765 \ 765 \
766 if (m1_nr == m2_nr && m1_nc == m2_nc) \ 766 if (m1_nr == m2_nr && m1_nc == m2_nc) \
767 { \ 767 { \
768 if (m1_nr != 0 || m1_nc != 0) \ 768 if (m1_nr != 0 || m1_nc != 0) \
769 { \ 769 { \
770 /* Count num of non-zero elements */ \ 770 /* Count num of non-zero elements */ \
771 int nel = 0; \ 771 octave_idx_type nel = 0; \
772 for (int j = 0; j < m1_nc; j++) \ 772 for (octave_idx_type j = 0; j < m1_nc; j++) \
773 for (int i = 0; i < m1_nr; i++) \ 773 for (octave_idx_type i = 0; i < m1_nr; i++) \
774 if ((m1.elem(i, j) != LHS_ZERO) \ 774 if ((m1.elem(i, j) != LHS_ZERO) \
775 OP (m2.elem(i, j) != RHS_ZERO)) \ 775 OP (m2.elem(i, j) != RHS_ZERO)) \
776 nel++; \ 776 nel++; \
777 \ 777 \
778 r = SparseBoolMatrix (m1_nr, m1_nc, nel); \ 778 r = SparseBoolMatrix (m1_nr, m1_nc, nel); \
779 \ 779 \
780 int ii = 0; \ 780 octave_idx_type ii = 0; \
781 r.cidx (0) = 0; \ 781 r.cidx (0) = 0; \
782 for (int j = 0; j < m1_nc; j++) \ 782 for (octave_idx_type j = 0; j < m1_nc; j++) \
783 { \ 783 { \
784 for (int i = 0; i < m1_nr; i++) \ 784 for (octave_idx_type i = 0; i < m1_nr; i++) \
785 { \ 785 { \
786 bool el = (m1.elem(i, j) != LHS_ZERO) \ 786 bool el = (m1.elem(i, j) != LHS_ZERO) \
787 OP (m2.elem(i, j) != RHS_ZERO); \ 787 OP (m2.elem(i, j) != RHS_ZERO); \
788 if (el) \ 788 if (el) \
789 { \ 789 { \
827 R \ 827 R \
828 F (const M1& m1, const M2& m2) \ 828 F (const M1& m1, const M2& m2) \
829 { \ 829 { \
830 R r; \ 830 R r; \
831 \ 831 \
832 int m1_nr = m1.rows (); \ 832 octave_idx_type m1_nr = m1.rows (); \
833 int m1_nc = m1.cols (); \ 833 octave_idx_type m1_nc = m1.cols (); \
834 \ 834 \
835 int m2_nr = m2.rows (); \ 835 octave_idx_type m2_nr = m2.rows (); \
836 int m2_nc = m2.cols (); \ 836 octave_idx_type m2_nc = m2.cols (); \
837 \ 837 \
838 if (m1_nr != m2_nr || m1_nc != m2_nc) \ 838 if (m1_nr != m2_nr || m1_nc != m2_nc) \
839 gripe_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \ 839 gripe_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \
840 else \ 840 else \
841 { \ 841 { \
842 r = R (m1_nr, m1_nc); \ 842 r = R (m1_nr, m1_nc); \
843 \ 843 \
844 for (int j = 0; j < m1_nc; j++) \ 844 for (octave_idx_type j = 0; j < m1_nc; j++) \
845 for (int i = 0; i < m1_nr; i++) \ 845 for (octave_idx_type i = 0; i < m1_nr; i++) \
846 r.elem (i, j) = m1.elem (i, j) OP m2.elem (i, j); \ 846 r.elem (i, j) = m1.elem (i, j) OP m2.elem (i, j); \
847 } \ 847 } \
848 return r; \ 848 return r; \
849 } 849 }
850 850
852 R \ 852 R \
853 F (const M1& m1, const M2& m2) \ 853 F (const M1& m1, const M2& m2) \
854 { \ 854 { \
855 R r; \ 855 R r; \
856 \ 856 \
857 int m1_nr = m1.rows (); \ 857 octave_idx_type m1_nr = m1.rows (); \
858 int m1_nc = m1.cols (); \ 858 octave_idx_type m1_nc = m1.cols (); \
859 \ 859 \
860 int m2_nr = m2.rows (); \ 860 octave_idx_type m2_nr = m2.rows (); \
861 int m2_nc = m2.cols (); \ 861 octave_idx_type m2_nc = m2.cols (); \
862 \ 862 \
863 if (m1_nr != m2_nr || m1_nc != m2_nc) \ 863 if (m1_nr != m2_nr || m1_nc != m2_nc) \
864 gripe_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \ 864 gripe_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \
865 else \ 865 else \
866 { \ 866 { \
867 /* Count num of non-zero elements */ \ 867 /* Count num of non-zero elements */ \
868 int nel = 0; \ 868 octave_idx_type nel = 0; \
869 for (int j = 0; j < m1_nc; j++) \ 869 for (octave_idx_type j = 0; j < m1_nc; j++) \
870 for (int i = 0; i < m1_nr; i++) \ 870 for (octave_idx_type i = 0; i < m1_nr; i++) \
871 if ((m1.elem(i, j) OP m2.elem(i, j)) != ZERO) \ 871 if ((m1.elem(i, j) OP m2.elem(i, j)) != ZERO) \
872 nel++; \ 872 nel++; \
873 \ 873 \
874 r = R (m1_nr, m1_nc, nel); \ 874 r = R (m1_nr, m1_nc, nel); \
875 \ 875 \
876 int ii = 0; \ 876 octave_idx_type ii = 0; \
877 r.cidx (0) = 0; \ 877 r.cidx (0) = 0; \
878 for (int j = 0 ; j < m1_nc ; j++) \ 878 for (octave_idx_type j = 0 ; j < m1_nc ; j++) \
879 { \ 879 { \
880 for (int i = 0 ; i < m1_nr ; i++) \ 880 for (octave_idx_type i = 0 ; i < m1_nr ; i++) \
881 { \ 881 { \
882 if ((m1.elem(i, j) OP m2.elem(i, j)) != ZERO) \ 882 if ((m1.elem(i, j) OP m2.elem(i, j)) != ZERO) \
883 { \ 883 { \
884 r.data (ii) = m1.elem(i, j) OP m2.elem(i,j); \ 884 r.data (ii) = m1.elem(i, j) OP m2.elem(i,j); \
885 r.ridx (ii++) = i; \ 885 r.ridx (ii++) = i; \
915 SparseBoolMatrix \ 915 SparseBoolMatrix \
916 F (const M1& m1, const M2& m2) \ 916 F (const M1& m1, const M2& m2) \
917 { \ 917 { \
918 SparseBoolMatrix r; \ 918 SparseBoolMatrix r; \
919 \ 919 \
920 int m1_nr = m1.rows (); \ 920 octave_idx_type m1_nr = m1.rows (); \
921 int m1_nc = m1.cols (); \ 921 octave_idx_type m1_nc = m1.cols (); \
922 \ 922 \
923 int m2_nr = m2.rows (); \ 923 octave_idx_type m2_nr = m2.rows (); \
924 int m2_nc = m2.cols (); \ 924 octave_idx_type m2_nc = m2.cols (); \
925 \ 925 \
926 if (m1_nr == m2_nr && m1_nc == m2_nc) \ 926 if (m1_nr == m2_nr && m1_nc == m2_nc) \
927 { \ 927 { \
928 if (m1_nr != 0 || m1_nc != 0) \ 928 if (m1_nr != 0 || m1_nc != 0) \
929 { \ 929 { \
930 /* Count num of non-zero elements */ \ 930 /* Count num of non-zero elements */ \
931 int nel = 0; \ 931 octave_idx_type nel = 0; \
932 for (int j = 0; j < m1_nc; j++) \ 932 for (octave_idx_type j = 0; j < m1_nc; j++) \
933 for (int i = 0; i < m1_nr; i++) \ 933 for (octave_idx_type i = 0; i < m1_nr; i++) \
934 if (C1 (m1.elem(i, j)) OP C2 (m2.elem(i, j))) \ 934 if (C1 (m1.elem(i, j)) OP C2 (m2.elem(i, j))) \
935 nel++; \ 935 nel++; \
936 \ 936 \
937 r = SparseBoolMatrix (m1_nr, m1_nc, nel); \ 937 r = SparseBoolMatrix (m1_nr, m1_nc, nel); \
938 \ 938 \
939 int ii = 0; \ 939 octave_idx_type ii = 0; \
940 r.cidx (0) = 0; \ 940 r.cidx (0) = 0; \
941 for (int j = 0; j < m1_nc; j++) \ 941 for (octave_idx_type j = 0; j < m1_nc; j++) \
942 { \ 942 { \
943 for (int i = 0; i < m1_nr; i++) \ 943 for (octave_idx_type i = 0; i < m1_nr; i++) \
944 { \ 944 { \
945 bool el = C1 (m1.elem(i, j)) OP C2 (m2.elem(i, j)); \ 945 bool el = C1 (m1.elem(i, j)) OP C2 (m2.elem(i, j)); \
946 if (el) \ 946 if (el) \
947 { \ 947 { \
948 r.data(ii) = el; \ 948 r.data(ii) = el; \
981 SparseBoolMatrix \ 981 SparseBoolMatrix \
982 F (const M1& m1, const M2& m2) \ 982 F (const M1& m1, const M2& m2) \
983 { \ 983 { \
984 SparseBoolMatrix r; \ 984 SparseBoolMatrix r; \
985 \ 985 \
986 int m1_nr = m1.rows (); \ 986 octave_idx_type m1_nr = m1.rows (); \
987 int m1_nc = m1.cols (); \ 987 octave_idx_type m1_nc = m1.cols (); \
988 \ 988 \
989 int m2_nr = m2.rows (); \ 989 octave_idx_type m2_nr = m2.rows (); \
990 int m2_nc = m2.cols (); \ 990 octave_idx_type m2_nc = m2.cols (); \
991 \ 991 \
992 if (m1_nr == m2_nr && m1_nc == m2_nc) \ 992 if (m1_nr == m2_nr && m1_nc == m2_nc) \
993 { \ 993 { \
994 if (m1_nr != 0 || m1_nc != 0) \ 994 if (m1_nr != 0 || m1_nc != 0) \
995 { \ 995 { \
996 /* Count num of non-zero elements */ \ 996 /* Count num of non-zero elements */ \
997 int nel = 0; \ 997 octave_idx_type nel = 0; \
998 for (int j = 0; j < m1_nc; j++) \ 998 for (octave_idx_type j = 0; j < m1_nc; j++) \
999 for (int i = 0; i < m1_nr; i++) \ 999 for (octave_idx_type i = 0; i < m1_nr; i++) \
1000 if ((m1.elem(i, j) != LHS_ZERO) \ 1000 if ((m1.elem(i, j) != LHS_ZERO) \
1001 OP (m2.elem(i, j) != RHS_ZERO)) \ 1001 OP (m2.elem(i, j) != RHS_ZERO)) \
1002 nel++; \ 1002 nel++; \
1003 \ 1003 \
1004 r = SparseBoolMatrix (m1_nr, m1_nc, nel); \ 1004 r = SparseBoolMatrix (m1_nr, m1_nc, nel); \
1005 \ 1005 \
1006 int ii = 0; \ 1006 octave_idx_type ii = 0; \
1007 r.cidx (0) = 0; \ 1007 r.cidx (0) = 0; \
1008 for (int j = 0; j < m1_nc; j++) \ 1008 for (octave_idx_type j = 0; j < m1_nc; j++) \
1009 { \ 1009 { \
1010 for (int i = 0; i < m1_nr; i++) \ 1010 for (octave_idx_type i = 0; i < m1_nr; i++) \
1011 { \ 1011 { \
1012 bool el = (m1.elem(i, j) != LHS_ZERO) \ 1012 bool el = (m1.elem(i, j) != LHS_ZERO) \
1013 OP (m2.elem(i, j) != RHS_ZERO); \ 1013 OP (m2.elem(i, j) != RHS_ZERO); \
1014 if (el) \ 1014 if (el) \
1015 { \ 1015 { \
1053 R \ 1053 R \
1054 F (const M1& m1, const M2& m2) \ 1054 F (const M1& m1, const M2& m2) \
1055 { \ 1055 { \
1056 R r; \ 1056 R r; \
1057 \ 1057 \
1058 int m1_nr = m1.rows (); \ 1058 octave_idx_type m1_nr = m1.rows (); \
1059 int m1_nc = m1.cols (); \ 1059 octave_idx_type m1_nc = m1.cols (); \
1060 \ 1060 \
1061 int m2_nr = m2.rows (); \ 1061 octave_idx_type m2_nr = m2.rows (); \
1062 int m2_nc = m2.cols (); \ 1062 octave_idx_type m2_nc = m2.cols (); \
1063 \ 1063 \
1064 if (m1_nr != m2_nr || m1_nc != m2_nc) \ 1064 if (m1_nr != m2_nr || m1_nc != m2_nc) \
1065 gripe_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \ 1065 gripe_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \
1066 else \ 1066 else \
1067 { \ 1067 { \
1068 r = R (m1_nr, m1_nc); \ 1068 r = R (m1_nr, m1_nc); \
1069 \ 1069 \
1070 for (int j = 0; j < m1_nc; j++) \ 1070 for (octave_idx_type j = 0; j < m1_nc; j++) \
1071 for (int i = 0; i < m1_nr; i++) \ 1071 for (octave_idx_type i = 0; i < m1_nr; i++) \
1072 r.elem (i, j) = m1.elem (i, j) OP m2.elem (i, j); \ 1072 r.elem (i, j) = m1.elem (i, j) OP m2.elem (i, j); \
1073 } \ 1073 } \
1074 return r; \ 1074 return r; \
1075 } 1075 }
1076 1076
1078 R \ 1078 R \
1079 F (const M1& m1, const M2& m2) \ 1079 F (const M1& m1, const M2& m2) \
1080 { \ 1080 { \
1081 R r; \ 1081 R r; \
1082 \ 1082 \
1083 int m1_nr = m1.rows (); \ 1083 octave_idx_type m1_nr = m1.rows (); \
1084 int m1_nc = m1.cols (); \ 1084 octave_idx_type m1_nc = m1.cols (); \
1085 \ 1085 \
1086 int m2_nr = m2.rows (); \ 1086 octave_idx_type m2_nr = m2.rows (); \
1087 int m2_nc = m2.cols (); \ 1087 octave_idx_type m2_nc = m2.cols (); \
1088 \ 1088 \
1089 if (m1_nr != m2_nr || m1_nc != m2_nc) \ 1089 if (m1_nr != m2_nr || m1_nc != m2_nc) \
1090 gripe_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \ 1090 gripe_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \
1091 else \ 1091 else \
1092 { \ 1092 { \
1093 /* Count num of non-zero elements */ \ 1093 /* Count num of non-zero elements */ \
1094 int nel = 0; \ 1094 octave_idx_type nel = 0; \
1095 for (int j = 0; j < m1_nc; j++) \ 1095 for (octave_idx_type j = 0; j < m1_nc; j++) \
1096 for (int i = 0; i < m1_nr; i++) \ 1096 for (octave_idx_type i = 0; i < m1_nr; i++) \
1097 if ((m1.elem(i, j) OP m2.elem(i, j)) != ZERO) \ 1097 if ((m1.elem(i, j) OP m2.elem(i, j)) != ZERO) \
1098 nel++; \ 1098 nel++; \
1099 \ 1099 \
1100 r = R (m1_nr, m1_nc, nel); \ 1100 r = R (m1_nr, m1_nc, nel); \
1101 \ 1101 \
1102 int ii = 0; \ 1102 octave_idx_type ii = 0; \
1103 r.cidx (0) = 0; \ 1103 r.cidx (0) = 0; \
1104 for (int j = 0 ; j < m1_nc ; j++) \ 1104 for (octave_idx_type j = 0 ; j < m1_nc ; j++) \
1105 { \ 1105 { \
1106 for (int i = 0 ; i < m1_nr ; i++) \ 1106 for (octave_idx_type i = 0 ; i < m1_nr ; i++) \
1107 { \ 1107 { \
1108 if ((m1.elem(i, j) OP m2.elem(i, j)) != ZERO) \ 1108 if ((m1.elem(i, j) OP m2.elem(i, j)) != ZERO) \
1109 { \ 1109 { \
1110 r.data (ii) = m1.elem(i, j) OP m2.elem(i,j); \ 1110 r.data (ii) = m1.elem(i, j) OP m2.elem(i,j); \
1111 r.ridx (ii++) = i; \ 1111 r.ridx (ii++) = i; \
1141 SparseBoolMatrix \ 1141 SparseBoolMatrix \
1142 F (const M1& m1, const M2& m2) \ 1142 F (const M1& m1, const M2& m2) \
1143 { \ 1143 { \
1144 SparseBoolMatrix r; \ 1144 SparseBoolMatrix r; \
1145 \ 1145 \
1146 int m1_nr = m1.rows (); \ 1146 octave_idx_type m1_nr = m1.rows (); \
1147 int m1_nc = m1.cols (); \ 1147 octave_idx_type m1_nc = m1.cols (); \
1148 \ 1148 \
1149 int m2_nr = m2.rows (); \ 1149 octave_idx_type m2_nr = m2.rows (); \
1150 int m2_nc = m2.cols (); \ 1150 octave_idx_type m2_nc = m2.cols (); \
1151 \ 1151 \
1152 if (m1_nr == m2_nr && m1_nc == m2_nc) \ 1152 if (m1_nr == m2_nr && m1_nc == m2_nc) \
1153 { \ 1153 { \
1154 if (m1_nr != 0 || m1_nc != 0) \ 1154 if (m1_nr != 0 || m1_nc != 0) \
1155 { \ 1155 { \
1156 /* Count num of non-zero elements */ \ 1156 /* Count num of non-zero elements */ \
1157 int nel = 0; \ 1157 octave_idx_type nel = 0; \
1158 for (int j = 0; j < m1_nc; j++) \ 1158 for (octave_idx_type j = 0; j < m1_nc; j++) \
1159 for (int i = 0; i < m1_nr; i++) \ 1159 for (octave_idx_type i = 0; i < m1_nr; i++) \
1160 if (C1 (m1.elem(i, j)) OP C2 (m2.elem(i, j))) \ 1160 if (C1 (m1.elem(i, j)) OP C2 (m2.elem(i, j))) \
1161 nel++; \ 1161 nel++; \
1162 \ 1162 \
1163 r = SparseBoolMatrix (m1_nr, m1_nc, nel); \ 1163 r = SparseBoolMatrix (m1_nr, m1_nc, nel); \
1164 \ 1164 \
1165 int ii = 0; \ 1165 octave_idx_type ii = 0; \
1166 r.cidx (0) = 0; \ 1166 r.cidx (0) = 0; \
1167 for (int j = 0; j < m1_nc; j++) \ 1167 for (octave_idx_type j = 0; j < m1_nc; j++) \
1168 { \ 1168 { \
1169 for (int i = 0; i < m1_nr; i++) \ 1169 for (octave_idx_type i = 0; i < m1_nr; i++) \
1170 { \ 1170 { \
1171 bool el = C1 (m1.elem(i, j)) OP C2 (m2.elem(i, j)); \ 1171 bool el = C1 (m1.elem(i, j)) OP C2 (m2.elem(i, j)); \
1172 if (el) \ 1172 if (el) \
1173 { \ 1173 { \
1174 r.data(ii) = el; \ 1174 r.data(ii) = el; \
1207 SparseBoolMatrix \ 1207 SparseBoolMatrix \
1208 F (const M1& m1, const M2& m2) \ 1208 F (const M1& m1, const M2& m2) \
1209 { \ 1209 { \
1210 SparseBoolMatrix r; \ 1210 SparseBoolMatrix r; \
1211 \ 1211 \
1212 int m1_nr = m1.rows (); \ 1212 octave_idx_type m1_nr = m1.rows (); \
1213 int m1_nc = m1.cols (); \ 1213 octave_idx_type m1_nc = m1.cols (); \
1214 \ 1214 \
1215 int m2_nr = m2.rows (); \ 1215 octave_idx_type m2_nr = m2.rows (); \
1216 int m2_nc = m2.cols (); \ 1216 octave_idx_type m2_nc = m2.cols (); \
1217 \ 1217 \
1218 if (m1_nr == m2_nr && m1_nc == m2_nc) \ 1218 if (m1_nr == m2_nr && m1_nc == m2_nc) \
1219 { \ 1219 { \
1220 if (m1_nr != 0 || m1_nc != 0) \ 1220 if (m1_nr != 0 || m1_nc != 0) \
1221 { \ 1221 { \
1222 /* Count num of non-zero elements */ \ 1222 /* Count num of non-zero elements */ \
1223 int nel = 0; \ 1223 octave_idx_type nel = 0; \
1224 for (int j = 0; j < m1_nc; j++) \ 1224 for (octave_idx_type j = 0; j < m1_nc; j++) \
1225 for (int i = 0; i < m1_nr; i++) \ 1225 for (octave_idx_type i = 0; i < m1_nr; i++) \
1226 if ((m1.elem(i, j) != LHS_ZERO) \ 1226 if ((m1.elem(i, j) != LHS_ZERO) \
1227 OP (m2.elem(i, j) != RHS_ZERO)) \ 1227 OP (m2.elem(i, j) != RHS_ZERO)) \
1228 nel++; \ 1228 nel++; \
1229 \ 1229 \
1230 r = SparseBoolMatrix (m1_nr, m1_nc, nel); \ 1230 r = SparseBoolMatrix (m1_nr, m1_nc, nel); \
1231 \ 1231 \
1232 int ii = 0; \ 1232 octave_idx_type ii = 0; \
1233 r.cidx (0) = 0; \ 1233 r.cidx (0) = 0; \
1234 for (int j = 0; j < m1_nc; j++) \ 1234 for (octave_idx_type j = 0; j < m1_nc; j++) \
1235 { \ 1235 { \
1236 for (int i = 0; i < m1_nr; i++) \ 1236 for (octave_idx_type i = 0; i < m1_nr; i++) \
1237 { \ 1237 { \
1238 bool el = (m1.elem(i, j) != LHS_ZERO) \ 1238 bool el = (m1.elem(i, j) != LHS_ZERO) \
1239 OP (m2.elem(i, j) != RHS_ZERO); \ 1239 OP (m2.elem(i, j) != RHS_ZERO); \
1240 if (el) \ 1240 if (el) \
1241 { \ 1241 { \
1269 1269
1270 // Avoid some code duplication. Maybe we should use templates. 1270 // Avoid some code duplication. Maybe we should use templates.
1271 1271
1272 #define SPARSE_CUMSUM(RET_TYPE, ELT_TYPE, FCN) \ 1272 #define SPARSE_CUMSUM(RET_TYPE, ELT_TYPE, FCN) \
1273 \ 1273 \
1274 int nr = rows (); \ 1274 octave_idx_type nr = rows (); \
1275 int nc = cols (); \ 1275 octave_idx_type nc = cols (); \
1276 \ 1276 \
1277 RET_TYPE retval; \ 1277 RET_TYPE retval; \
1278 \ 1278 \
1279 if (nr > 0 && nc > 0) \ 1279 if (nr > 0 && nc > 0) \
1280 { \ 1280 { \
1281 if ((nr == 1 && dim == -1) || dim == 1) \ 1281 if ((nr == 1 && dim == -1) || dim == 1) \
1282 /* Ugly!! Is there a better way? */ \ 1282 /* Ugly!! Is there a better way? */ \
1283 retval = transpose (). FCN (0) .transpose (); \ 1283 retval = transpose (). FCN (0) .transpose (); \
1284 else \ 1284 else \
1285 { \ 1285 { \
1286 int nel = 0; \ 1286 octave_idx_type nel = 0; \
1287 for (int i = 0; i < nc; i++) \ 1287 for (octave_idx_type i = 0; i < nc; i++) \
1288 { \ 1288 { \
1289 ELT_TYPE t = ELT_TYPE (); \ 1289 ELT_TYPE t = ELT_TYPE (); \
1290 for (int j = cidx (i); j < cidx (i+1); j++) \ 1290 for (octave_idx_type j = cidx (i); j < cidx (i+1); j++) \
1291 { \ 1291 { \
1292 t += data(j); \ 1292 t += data(j); \
1293 if (t != ELT_TYPE ()) \ 1293 if (t != ELT_TYPE ()) \
1294 if (j == cidx(i+1) - 1) \ 1294 if (j == cidx(i+1) - 1) \
1295 nel += nr - ridx(j); \ 1295 nel += nr - ridx(j); \
1297 nel += ridx(j+1) - ridx(j); \ 1297 nel += ridx(j+1) - ridx(j); \
1298 } \ 1298 } \
1299 } \ 1299 } \
1300 retval = RET_TYPE (nr, nc, nel); \ 1300 retval = RET_TYPE (nr, nc, nel); \
1301 retval.cidx(0) = 0; \ 1301 retval.cidx(0) = 0; \
1302 int ii = 0; \ 1302 octave_idx_type ii = 0; \
1303 for (int i = 0; i < nc; i++) \ 1303 for (octave_idx_type i = 0; i < nc; i++) \
1304 { \ 1304 { \
1305 ELT_TYPE t = ELT_TYPE (); \ 1305 ELT_TYPE t = ELT_TYPE (); \
1306 for (int j = cidx (i); j < cidx (i+1); j++) \ 1306 for (octave_idx_type j = cidx (i); j < cidx (i+1); j++) \
1307 { \ 1307 { \
1308 t += data(j); \ 1308 t += data(j); \
1309 if (t != ELT_TYPE ()) \ 1309 if (t != ELT_TYPE ()) \
1310 { \ 1310 { \
1311 if (j == cidx(i+1) - 1) \ 1311 if (j == cidx(i+1) - 1) \
1312 { \ 1312 { \
1313 for (int k = ridx(j); k < nr; k++) \ 1313 for (octave_idx_type k = ridx(j); k < nr; k++) \
1314 { \ 1314 { \
1315 retval.data (ii) = t; \ 1315 retval.data (ii) = t; \
1316 retval.ridx (ii++) = k; \ 1316 retval.ridx (ii++) = k; \
1317 } \ 1317 } \
1318 } \ 1318 } \
1319 else \ 1319 else \
1320 { \ 1320 { \
1321 for (int k = ridx(j); k < ridx(j+1); k++) \ 1321 for (octave_idx_type k = ridx(j); k < ridx(j+1); k++) \
1322 { \ 1322 { \
1323 retval.data (ii) = t; \ 1323 retval.data (ii) = t; \
1324 retval.ridx (ii++) = k; \ 1324 retval.ridx (ii++) = k; \
1325 } \ 1325 } \
1326 } \ 1326 } \
1336 return retval 1336 return retval
1337 1337
1338 1338
1339 #define SPARSE_CUMPROD(RET_TYPE, ELT_TYPE, FCN) \ 1339 #define SPARSE_CUMPROD(RET_TYPE, ELT_TYPE, FCN) \
1340 \ 1340 \
1341 int nr = rows (); \ 1341 octave_idx_type nr = rows (); \
1342 int nc = cols (); \ 1342 octave_idx_type nc = cols (); \
1343 \ 1343 \
1344 RET_TYPE retval; \ 1344 RET_TYPE retval; \
1345 \ 1345 \
1346 if (nr > 0 && nc > 0) \ 1346 if (nr > 0 && nc > 0) \
1347 { \ 1347 { \
1348 if ((nr == 1 && dim == -1) || dim == 1) \ 1348 if ((nr == 1 && dim == -1) || dim == 1) \
1349 /* Ugly!! Is there a better way? */ \ 1349 /* Ugly!! Is there a better way? */ \
1350 retval = transpose (). FCN (0) .transpose (); \ 1350 retval = transpose (). FCN (0) .transpose (); \
1351 else \ 1351 else \
1352 { \ 1352 { \
1353 int nel = 0; \ 1353 octave_idx_type nel = 0; \
1354 for (int i = 0; i < nc; i++) \ 1354 for (octave_idx_type i = 0; i < nc; i++) \
1355 { \ 1355 { \
1356 int jj = 0; \ 1356 octave_idx_type jj = 0; \
1357 for (int j = cidx (i); j < cidx (i+1); j++) \ 1357 for (octave_idx_type j = cidx (i); j < cidx (i+1); j++) \
1358 { \ 1358 { \
1359 if (jj == ridx(j)) \ 1359 if (jj == ridx(j)) \
1360 { \ 1360 { \
1361 nel++; \ 1361 nel++; \
1362 jj++; \ 1362 jj++; \
1365 break; \ 1365 break; \
1366 } \ 1366 } \
1367 } \ 1367 } \
1368 retval = RET_TYPE (nr, nc, nel); \ 1368 retval = RET_TYPE (nr, nc, nel); \
1369 retval.cidx(0) = 0; \ 1369 retval.cidx(0) = 0; \
1370 int ii = 0; \ 1370 octave_idx_type ii = 0; \
1371 for (int i = 0; i < nc; i++) \ 1371 for (octave_idx_type i = 0; i < nc; i++) \
1372 { \ 1372 { \
1373 ELT_TYPE t = ELT_TYPE (1.); \ 1373 ELT_TYPE t = ELT_TYPE (1.); \
1374 int jj = 0; \ 1374 octave_idx_type jj = 0; \
1375 for (int j = cidx (i); j < cidx (i+1); j++) \ 1375 for (octave_idx_type j = cidx (i); j < cidx (i+1); j++) \
1376 { \ 1376 { \
1377 if (jj == ridx(j)) \ 1377 if (jj == ridx(j)) \
1378 { \ 1378 { \
1379 t *= data(j); \ 1379 t *= data(j); \
1380 retval.data(ii) = t; \ 1380 retval.data(ii) = t; \
1393 return retval 1393 return retval
1394 1394
1395 #define SPARSE_BASE_REDUCTION_OP(RET_TYPE, EL_TYPE, ROW_EXPR, COL_EXPR, \ 1395 #define SPARSE_BASE_REDUCTION_OP(RET_TYPE, EL_TYPE, ROW_EXPR, COL_EXPR, \
1396 INIT_VAL, MT_RESULT) \ 1396 INIT_VAL, MT_RESULT) \
1397 \ 1397 \
1398 int nr = rows (); \ 1398 octave_idx_type nr = rows (); \
1399 int nc = cols (); \ 1399 octave_idx_type nc = cols (); \
1400 \ 1400 \
1401 RET_TYPE retval; \ 1401 RET_TYPE retval; \
1402 \ 1402 \
1403 if (nr > 0 && nc > 0) \ 1403 if (nr > 0 && nc > 0) \
1404 { \ 1404 { \
1405 if ((nr == 1 && dim == -1) || dim == 1) \ 1405 if ((nr == 1 && dim == -1) || dim == 1) \
1406 { \ 1406 { \
1407 OCTAVE_LOCAL_BUFFER (EL_TYPE, tmp, nr); \ 1407 OCTAVE_LOCAL_BUFFER (EL_TYPE, tmp, nr); \
1408 \ 1408 \
1409 for (int i = 0; i < nr; i++) \ 1409 for (octave_idx_type i = 0; i < nr; i++) \
1410 { \ 1410 { \
1411 tmp[i] = INIT_VAL; \ 1411 tmp[i] = INIT_VAL; \
1412 for (int j = 0; j < nc; j++) \ 1412 for (octave_idx_type j = 0; j < nc; j++) \
1413 { \ 1413 { \
1414 ROW_EXPR; \ 1414 ROW_EXPR; \
1415 } \ 1415 } \
1416 } \ 1416 } \
1417 int nel = 0; \ 1417 octave_idx_type nel = 0; \
1418 for (int i = 0; i < nr; i++) \ 1418 for (octave_idx_type i = 0; i < nr; i++) \
1419 if (tmp[i] != EL_TYPE ()) \ 1419 if (tmp[i] != EL_TYPE ()) \
1420 nel++ ; \ 1420 nel++ ; \
1421 retval = RET_TYPE (nr, 1, nel); \ 1421 retval = RET_TYPE (nr, static_cast<octave_idx_type> (1), nel); \
1422 retval.cidx(0) = 0; \ 1422 retval.cidx(0) = 0; \
1423 retval.cidx(1) = nel; \ 1423 retval.cidx(1) = nel; \
1424 nel = 0; \ 1424 nel = 0; \
1425 for (int i = 0; i < nr; i++) \ 1425 for (octave_idx_type i = 0; i < nr; i++) \
1426 if (tmp[i] != EL_TYPE ()) \ 1426 if (tmp[i] != EL_TYPE ()) \
1427 { \ 1427 { \
1428 retval.data(nel) = tmp[i]; \ 1428 retval.data(nel) = tmp[i]; \
1429 retval.ridx(nel++) = i; \ 1429 retval.ridx(nel++) = i; \
1430 } \ 1430 } \
1431 } \ 1431 } \
1432 else \ 1432 else \
1433 { \ 1433 { \
1434 OCTAVE_LOCAL_BUFFER (EL_TYPE, tmp, nc); \ 1434 OCTAVE_LOCAL_BUFFER (EL_TYPE, tmp, nc); \
1435 \ 1435 \
1436 for (int j = 0; j < nc; j++) \ 1436 for (octave_idx_type j = 0; j < nc; j++) \
1437 { \ 1437 { \
1438 tmp[j] = INIT_VAL; \ 1438 tmp[j] = INIT_VAL; \
1439 for (int i = 0; i < nr; i++) \ 1439 for (octave_idx_type i = 0; i < nr; i++) \
1440 { \ 1440 { \
1441 COL_EXPR; \ 1441 COL_EXPR; \
1442 } \ 1442 } \
1443 } \ 1443 } \
1444 int nel = 0; \ 1444 octave_idx_type nel = 0; \
1445 for (int i = 0; i < nc; i++) \ 1445 for (octave_idx_type i = 0; i < nc; i++) \
1446 if (tmp[i] != EL_TYPE ()) \ 1446 if (tmp[i] != EL_TYPE ()) \
1447 nel++ ; \ 1447 nel++ ; \
1448 retval = RET_TYPE (1, nc, nel); \ 1448 retval = RET_TYPE (static_cast<octave_idx_type> (1), nc, nel); \
1449 retval.cidx(0) = 0; \ 1449 retval.cidx(0) = 0; \
1450 nel = 0; \ 1450 nel = 0; \
1451 for (int i = 0; i < nc; i++) \ 1451 for (octave_idx_type i = 0; i < nc; i++) \
1452 if (tmp[i] != EL_TYPE ()) \ 1452 if (tmp[i] != EL_TYPE ()) \
1453 { \ 1453 { \
1454 retval.data(nel) = tmp[i]; \ 1454 retval.data(nel) = tmp[i]; \
1455 retval.ridx(nel++) = 0; \ 1455 retval.ridx(nel++) = 0; \
1456 retval.cidx(i+1) = retval.cidx(i) + 1; \ 1456 retval.cidx(i+1) = retval.cidx(i) + 1; \
1459 retval.cidx(i+1) = retval.cidx(i); \ 1459 retval.cidx(i+1) = retval.cidx(i); \
1460 } \ 1460 } \
1461 } \ 1461 } \
1462 else if (nc == 0 && (nr == 0 || (nr == 1 && dim == -1))) \ 1462 else if (nc == 0 && (nr == 0 || (nr == 1 && dim == -1))) \
1463 { \ 1463 { \
1464 retval = RET_TYPE (1, 1, 1); \ 1464 retval = RET_TYPE (static_cast<octave_idx_type> (1), \
1465 static_cast<octave_idx_type> (1), \
1466 static_cast<octave_idx_type> (1)); \
1465 retval.cidx(0) = 0; \ 1467 retval.cidx(0) = 0; \
1466 retval.cidx(1) = 1; \ 1468 retval.cidx(1) = 1; \
1467 retval.ridx(0) = 0; \ 1469 retval.ridx(0) = 0; \
1468 retval.data(0) = MT_RESULT; \ 1470 retval.data(0) = MT_RESULT; \
1469 } \ 1471 } \
1470 else if (nr == 0 && (dim == 0 || dim == -1)) \ 1472 else if (nr == 0 && (dim == 0 || dim == -1)) \
1471 { \ 1473 { \
1472 retval = RET_TYPE (1, nc, nc); \ 1474 retval = RET_TYPE (static_cast<octave_idx_type> (1), nc, nc); \
1473 retval.cidx (0) = 0; \ 1475 retval.cidx (0) = 0; \
1474 for (int i = 0; i < nc ; i++) \ 1476 for (octave_idx_type i = 0; i < nc ; i++) \
1475 { \ 1477 { \
1476 retval.ridx (i) = 0; \ 1478 retval.ridx (i) = 0; \
1477 retval.cidx (i+1) = i; \ 1479 retval.cidx (i+1) = i; \
1478 retval.data (i) = MT_RESULT; \ 1480 retval.data (i) = MT_RESULT; \
1479 } \ 1481 } \
1480 } \ 1482 } \
1481 else if (nc == 0 && dim == 1) \ 1483 else if (nc == 0 && dim == 1) \
1482 { \ 1484 { \
1483 retval = RET_TYPE (nr, 1, nr); \ 1485 retval = RET_TYPE (nr, static_cast<octave_idx_type> (1), nr); \
1484 retval.cidx(0) = 0; \ 1486 retval.cidx(0) = 0; \
1485 retval.cidx(1) = nr; \ 1487 retval.cidx(1) = nr; \
1486 for (int i = 0; i < nr; i++) \ 1488 for (octave_idx_type i = 0; i < nr; i++) \
1487 { \ 1489 { \
1488 retval.ridx(i) = i; \ 1490 retval.ridx(i) = i; \
1489 retval.data(i) = MT_RESULT; \ 1491 retval.data(i) = MT_RESULT; \
1490 } \ 1492 } \
1491 } \ 1493 } \
1529 #define SPARSE_ALL_OP(DIM) SPARSE_ANY_ALL_OP (DIM, true, ==, false) 1531 #define SPARSE_ALL_OP(DIM) SPARSE_ANY_ALL_OP (DIM, true, ==, false)
1530 1532
1531 #define SPARSE_ANY_OP(DIM) SPARSE_ANY_ALL_OP (DIM, false, !=, true) 1533 #define SPARSE_ANY_OP(DIM) SPARSE_ANY_ALL_OP (DIM, false, !=, true)
1532 1534
1533 #define SPARSE_SPARSE_MUL( RET_TYPE, EL_TYPE ) \ 1535 #define SPARSE_SPARSE_MUL( RET_TYPE, EL_TYPE ) \
1534 int nr = m.rows (); \ 1536 octave_idx_type nr = m.rows (); \
1535 int nc = m.cols (); \ 1537 octave_idx_type nc = m.cols (); \
1536 \ 1538 \
1537 int a_nr = a.rows (); \ 1539 octave_idx_type a_nr = a.rows (); \
1538 int a_nc = a.cols (); \ 1540 octave_idx_type a_nc = a.cols (); \
1539 \ 1541 \
1540 if (nc != a_nr) \ 1542 if (nc != a_nr) \
1541 { \ 1543 { \
1542 gripe_nonconformant ("operator *", nr, nc, a_nr, a_nc); \ 1544 gripe_nonconformant ("operator *", nr, nc, a_nr, a_nc); \
1543 return RET_TYPE (); \ 1545 return RET_TYPE (); \
1544 } \ 1546 } \
1545 else \ 1547 else \
1546 { \ 1548 { \
1547 OCTAVE_LOCAL_BUFFER (EL_TYPE, Xcol, nr); \ 1549 OCTAVE_LOCAL_BUFFER (EL_TYPE, Xcol, nr); \
1548 \ 1550 \
1549 int nel = 0; \ 1551 octave_idx_type nel = 0; \
1550 \ 1552 \
1551 for (int i = 0; i < a_nc; i++) \ 1553 for (octave_idx_type i = 0; i < a_nc; i++) \
1552 { \ 1554 { \
1553 OCTAVE_QUIT; \ 1555 OCTAVE_QUIT; \
1554 for (int k = 0; k < nr; k++) \ 1556 for (octave_idx_type k = 0; k < nr; k++) \
1555 Xcol[k]= 0.; \ 1557 Xcol[k]= 0.; \
1556 for (int j = a.cidx(i); j < a.cidx(i+1); j++) \ 1558 for (octave_idx_type j = a.cidx(i); j < a.cidx(i+1); j++) \
1557 { \ 1559 { \
1558 int col = a.ridx(j); \ 1560 octave_idx_type col = a.ridx(j); \
1559 for (int k = m.cidx(col) ; k < m.cidx(col+1); k++) \ 1561 for (octave_idx_type k = m.cidx(col) ; k < m.cidx(col+1); k++) \
1560 if (Xcol[m.ridx(k)] == 0.) \ 1562 if (Xcol[m.ridx(k)] == 0.) \
1561 { \ 1563 { \
1562 Xcol[m.ridx(k)] = 1.; \ 1564 Xcol[m.ridx(k)] = 1.; \
1563 nel++; \ 1565 nel++; \
1564 } \ 1566 } \
1569 return RET_TYPE (nr, a_nc); \ 1571 return RET_TYPE (nr, a_nc); \
1570 else \ 1572 else \
1571 { \ 1573 { \
1572 RET_TYPE retval (nr, a_nc, nel); \ 1574 RET_TYPE retval (nr, a_nc, nel); \
1573 \ 1575 \
1574 int ii = 0; \ 1576 octave_idx_type ii = 0; \
1575 \ 1577 \
1576 retval.cidx(0) = 0; \ 1578 retval.cidx(0) = 0; \
1577 for (int i = 0; i < a_nc ; i++) \ 1579 for (octave_idx_type i = 0; i < a_nc ; i++) \
1578 { \ 1580 { \
1579 OCTAVE_QUIT; \ 1581 OCTAVE_QUIT; \
1580 for (int k = 0; k < nr; k++) \ 1582 for (octave_idx_type k = 0; k < nr; k++) \
1581 Xcol[k]= 0.; \ 1583 Xcol[k]= 0.; \
1582 for (int j = a.cidx(i); j < a.cidx(i+1); j++) \ 1584 for (octave_idx_type j = a.cidx(i); j < a.cidx(i+1); j++) \
1583 { \ 1585 { \
1584 int col = a.ridx(j); \ 1586 octave_idx_type col = a.ridx(j); \
1585 EL_TYPE tmpval = a.data(j); \ 1587 EL_TYPE tmpval = a.data(j); \
1586 for (int k = m.cidx(col) ; k < m.cidx(col+1); k++) \ 1588 for (octave_idx_type k = m.cidx(col) ; k < m.cidx(col+1); k++) \
1587 Xcol[m.ridx(k)] += tmpval * m.data(k); \ 1589 Xcol[m.ridx(k)] += tmpval * m.data(k); \
1588 } \ 1590 } \
1589 for (int k = 0; k < nr; k++) \ 1591 for (octave_idx_type k = 0; k < nr; k++) \
1590 { \ 1592 { \
1591 if (Xcol[k] !=0. ) \ 1593 if (Xcol[k] !=0. ) \
1592 { \ 1594 { \
1593 retval.ridx (ii) = k; \ 1595 retval.ridx (ii) = k; \
1594 retval.data (ii++) = Xcol[k]; \ 1596 retval.data (ii++) = Xcol[k]; \