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