Mercurial > hg > octave-nkf
annotate liboctave/MSparse.cc @ 10396:a0b51ac0f88a
optimize accumdim with summation
| author | Jaroslav Hajek <highegg@gmail.com> |
|---|---|
| date | Fri, 05 Mar 2010 12:31:30 +0100 |
| parents | 12884915a8e4 |
| children | fd0a3ac60b0e |
| 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 |
| 5 | |
| 6 This file is part of Octave. | |
| 5164 | 7 |
| 8 Octave is free software; you can redistribute it and/or modify it | |
| 9 under the terms of the GNU General Public License as published by the | |
| 7016 | 10 Free Software Foundation; either version 3 of the License, or (at your |
| 11 option) any later version. | |
| 5164 | 12 |
| 13 Octave is distributed in the hope that it will be useful, but WITHOUT | |
| 14 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or | |
| 15 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License | |
| 16 for more details. | |
| 17 | |
| 18 You should have received a copy of the GNU General Public License | |
| 7016 | 19 along with Octave; see the file COPYING. If not, see |
| 20 <http://www.gnu.org/licenses/>. | |
| 5164 | 21 |
| 22 */ | |
| 23 | |
| 24 #ifdef HAVE_CONFIG_H | |
| 25 #include <config.h> | |
| 26 #endif | |
| 27 | |
| 28 #include "quit.h" | |
| 29 #include "lo-error.h" | |
|
10350
12884915a8e4
merge MArray classes & improve Array interface
Jaroslav Hajek <highegg@gmail.com>
parents:
10314
diff
changeset
|
30 #include "MArray.h" |
| 5164 | 31 #include "Array-util.h" |
| 32 | |
| 33 #include "MSparse.h" | |
| 34 #include "MSparse-defs.h" | |
| 35 | |
| 36 // sparse array with math ops. | |
| 37 | |
| 38 // Element by element MSparse by MSparse ops. | |
| 39 | |
| 40 template <class T> | |
| 41 MSparse<T>& | |
| 42 operator += (MSparse<T>& a, const MSparse<T>& b) | |
| 43 { | |
| 44 MSparse<T> r; | |
| 45 | |
| 5275 | 46 octave_idx_type a_nr = a.rows (); |
| 47 octave_idx_type a_nc = a.cols (); | |
| 5164 | 48 |
| 5275 | 49 octave_idx_type b_nr = b.rows (); |
| 50 octave_idx_type b_nc = b.cols (); | |
| 5164 | 51 |
| 52 if (a_nr != b_nr || a_nc != b_nc) | |
| 53 gripe_nonconformant ("operator +=" , a_nr, a_nc, b_nr, b_nc); | |
| 54 else | |
| 55 { | |
| 5681 | 56 r = MSparse<T> (a_nr, a_nc, (a.nnz () + b.nnz ())); |
| 5164 | 57 |
| 5275 | 58 octave_idx_type jx = 0; |
| 59 for (octave_idx_type i = 0 ; i < a_nc ; i++) | |
| 5164 | 60 { |
| 5275 | 61 octave_idx_type ja = a.cidx(i); |
| 62 octave_idx_type ja_max = a.cidx(i+1); | |
| 5164 | 63 bool ja_lt_max= ja < ja_max; |
| 64 | |
| 5275 | 65 octave_idx_type jb = b.cidx(i); |
| 66 octave_idx_type jb_max = b.cidx(i+1); | |
| 5164 | 67 bool jb_lt_max = jb < jb_max; |
| 68 | |
| 69 while (ja_lt_max || jb_lt_max ) | |
| 70 { | |
|
10142
829e69ec3110
make OCTAVE_QUIT a function
Jaroslav Hajek <highegg@gmail.com>
parents:
7344
diff
changeset
|
71 octave_quit (); |
| 5164 | 72 if ((! jb_lt_max) || |
| 73 (ja_lt_max && (a.ridx(ja) < b.ridx(jb)))) | |
| 74 { | |
| 75 r.ridx(jx) = a.ridx(ja); | |
| 76 r.data(jx) = a.data(ja) + 0.; | |
| 77 jx++; | |
| 78 ja++; | |
| 79 ja_lt_max= ja < ja_max; | |
| 80 } | |
| 81 else if (( !ja_lt_max ) || | |
| 82 (jb_lt_max && (b.ridx(jb) < a.ridx(ja)) ) ) | |
| 83 { | |
|
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
84 r.ridx(jx) = b.ridx(jb); |
|
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
85 r.data(jx) = 0. + b.data(jb); |
|
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
86 jx++; |
| 5164 | 87 jb++; |
| 88 jb_lt_max= jb < jb_max; | |
| 89 } | |
| 90 else | |
| 91 { | |
|
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
92 if ((a.data(ja) + b.data(jb)) != 0.) |
|
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
93 { |
| 5164 | 94 r.data(jx) = a.data(ja) + b.data(jb); |
| 95 r.ridx(jx) = a.ridx(ja); | |
| 96 jx++; | |
| 97 } | |
| 98 ja++; | |
| 99 ja_lt_max= ja < ja_max; | |
| 100 jb++; | |
| 101 jb_lt_max= jb < jb_max; | |
| 102 } | |
| 103 } | |
| 104 r.cidx(i+1) = jx; | |
| 105 } | |
| 106 | |
|
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
107 a = r.maybe_compress (); |
| 5164 | 108 } |
| 109 | |
| 110 return a; | |
| 111 } | |
| 112 | |
| 113 template <class T> | |
| 114 MSparse<T>& | |
| 115 operator -= (MSparse<T>& a, const MSparse<T>& b) | |
| 116 { | |
| 117 MSparse<T> r; | |
| 118 | |
| 5275 | 119 octave_idx_type a_nr = a.rows (); |
| 120 octave_idx_type a_nc = a.cols (); | |
| 5164 | 121 |
| 5275 | 122 octave_idx_type b_nr = b.rows (); |
| 123 octave_idx_type b_nc = b.cols (); | |
| 5164 | 124 |
| 125 if (a_nr != b_nr || a_nc != b_nc) | |
| 126 gripe_nonconformant ("operator -=" , a_nr, a_nc, b_nr, b_nc); | |
| 127 else | |
| 128 { | |
| 5681 | 129 r = MSparse<T> (a_nr, a_nc, (a.nnz () + b.nnz ())); |
| 5164 | 130 |
| 5275 | 131 octave_idx_type jx = 0; |
| 132 for (octave_idx_type i = 0 ; i < a_nc ; i++) | |
| 5164 | 133 { |
| 5275 | 134 octave_idx_type ja = a.cidx(i); |
| 135 octave_idx_type ja_max = a.cidx(i+1); | |
| 5164 | 136 bool ja_lt_max= ja < ja_max; |
| 137 | |
| 5275 | 138 octave_idx_type jb = b.cidx(i); |
| 139 octave_idx_type jb_max = b.cidx(i+1); | |
| 5164 | 140 bool jb_lt_max = jb < jb_max; |
| 141 | |
| 142 while (ja_lt_max || jb_lt_max ) | |
| 143 { | |
|
10142
829e69ec3110
make OCTAVE_QUIT a function
Jaroslav Hajek <highegg@gmail.com>
parents:
7344
diff
changeset
|
144 octave_quit (); |
| 5164 | 145 if ((! jb_lt_max) || |
| 146 (ja_lt_max && (a.ridx(ja) < b.ridx(jb)))) | |
| 147 { | |
| 148 r.ridx(jx) = a.ridx(ja); | |
| 149 r.data(jx) = a.data(ja) - 0.; | |
| 150 jx++; | |
| 151 ja++; | |
| 152 ja_lt_max= ja < ja_max; | |
| 153 } | |
| 154 else if (( !ja_lt_max ) || | |
| 155 (jb_lt_max && (b.ridx(jb) < a.ridx(ja)) ) ) | |
| 156 { | |
|
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
157 r.ridx(jx) = b.ridx(jb); |
|
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
158 r.data(jx) = 0. - b.data(jb); |
|
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
159 jx++; |
| 5164 | 160 jb++; |
| 161 jb_lt_max= jb < jb_max; | |
| 162 } | |
| 163 else | |
| 164 { | |
|
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
165 if ((a.data(ja) - b.data(jb)) != 0.) |
|
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
166 { |
| 5164 | 167 r.data(jx) = a.data(ja) - b.data(jb); |
| 168 r.ridx(jx) = a.ridx(ja); | |
| 169 jx++; | |
| 170 } | |
| 171 ja++; | |
| 172 ja_lt_max= ja < ja_max; | |
| 173 jb++; | |
| 174 jb_lt_max= jb < jb_max; | |
| 175 } | |
| 176 } | |
| 177 r.cidx(i+1) = jx; | |
| 178 } | |
| 179 | |
|
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
180 a = r.maybe_compress (); |
| 5164 | 181 } |
| 182 | |
| 183 return a; | |
| 184 } | |
| 185 | |
| 186 // Element by element MSparse by scalar ops. | |
| 187 | |
| 188 #define SPARSE_A2S_OP_1(OP) \ | |
| 189 template <class T> \ | |
|
10350
12884915a8e4
merge MArray classes & improve Array interface
Jaroslav Hajek <highegg@gmail.com>
parents:
10314
diff
changeset
|
190 MArray<T> \ |
| 5164 | 191 operator OP (const MSparse<T>& a, const T& s) \ |
| 192 { \ | |
| 5275 | 193 octave_idx_type nr = a.rows (); \ |
| 194 octave_idx_type nc = a.cols (); \ | |
| 5164 | 195 \ |
|
10350
12884915a8e4
merge MArray classes & improve Array interface
Jaroslav Hajek <highegg@gmail.com>
parents:
10314
diff
changeset
|
196 MArray<T> r (nr, nc, (0.0 OP s)); \ |
| 5164 | 197 \ |
| 5275 | 198 for (octave_idx_type j = 0; j < nc; j++) \ |
|
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
199 for (octave_idx_type i = a.cidx(j); i < a.cidx(j+1); i++) \ |
|
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
200 r.elem (a.ridx (i), j) = a.data (i) OP s; \ |
| 5164 | 201 return r; \ |
| 202 } | |
| 203 | |
| 204 #define SPARSE_A2S_OP_2(OP) \ | |
| 205 template <class T> \ | |
| 206 MSparse<T> \ | |
| 207 operator OP (const MSparse<T>& a, const T& s) \ | |
| 208 { \ | |
| 5275 | 209 octave_idx_type nr = a.rows (); \ |
| 210 octave_idx_type nc = a.cols (); \ | |
| 5681 | 211 octave_idx_type nz = a.nnz (); \ |
| 5164 | 212 \ |
| 213 MSparse<T> r (nr, nc, nz); \ | |
| 214 \ | |
| 5275 | 215 for (octave_idx_type i = 0; i < nz; i++) \ |
| 5164 | 216 { \ |
|
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
217 r.data(i) = a.data(i) OP s; \ |
|
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
218 r.ridx(i) = a.ridx(i); \ |
| 5164 | 219 } \ |
| 5275 | 220 for (octave_idx_type i = 0; i < nc + 1; i++) \ |
| 5164 | 221 r.cidx(i) = a.cidx(i); \ |
| 222 r.maybe_compress (true); \ | |
| 223 return r; \ | |
| 224 } | |
| 225 | |
| 226 | |
| 227 SPARSE_A2S_OP_1 (+) | |
| 228 SPARSE_A2S_OP_1 (-) | |
| 229 SPARSE_A2S_OP_2 (*) | |
| 230 SPARSE_A2S_OP_2 (/) | |
| 231 | |
| 232 // Element by element scalar by MSparse ops. | |
| 233 | |
| 234 #define SPARSE_SA2_OP_1(OP) \ | |
| 235 template <class T> \ | |
|
10350
12884915a8e4
merge MArray classes & improve Array interface
Jaroslav Hajek <highegg@gmail.com>
parents:
10314
diff
changeset
|
236 MArray<T> \ |
| 5164 | 237 operator OP (const T& s, const MSparse<T>& a) \ |
| 238 { \ | |
| 5275 | 239 octave_idx_type nr = a.rows (); \ |
| 240 octave_idx_type nc = a.cols (); \ | |
| 5164 | 241 \ |
|
10350
12884915a8e4
merge MArray classes & improve Array interface
Jaroslav Hajek <highegg@gmail.com>
parents:
10314
diff
changeset
|
242 MArray<T> r (nr, nc, (s OP 0.0)); \ |
| 5164 | 243 \ |
| 5275 | 244 for (octave_idx_type j = 0; j < nc; j++) \ |
|
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
245 for (octave_idx_type i = a.cidx(j); i < a.cidx(j+1); i++) \ |
|
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
246 r.elem (a.ridx (i), j) = s OP a.data (i); \ |
| 5164 | 247 return r; \ |
| 248 } | |
| 249 | |
| 250 #define SPARSE_SA2_OP_2(OP) \ | |
| 251 template <class T> \ | |
| 252 MSparse<T> \ | |
| 253 operator OP (const T& s, const MSparse<T>& a) \ | |
| 254 { \ | |
| 5275 | 255 octave_idx_type nr = a.rows (); \ |
| 256 octave_idx_type nc = a.cols (); \ | |
| 5681 | 257 octave_idx_type nz = a.nnz (); \ |
| 5164 | 258 \ |
| 259 MSparse<T> r (nr, nc, nz); \ | |
| 260 \ | |
| 5275 | 261 for (octave_idx_type i = 0; i < nz; i++) \ |
| 5164 | 262 { \ |
|
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
263 r.data(i) = s OP a.data(i); \ |
|
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
264 r.ridx(i) = a.ridx(i); \ |
| 5164 | 265 } \ |
| 5275 | 266 for (octave_idx_type i = 0; i < nc + 1; i++) \ |
| 5164 | 267 r.cidx(i) = a.cidx(i); \ |
| 268 r.maybe_compress (true); \ | |
| 269 return r; \ | |
| 270 } | |
| 271 | |
| 272 SPARSE_SA2_OP_1 (+) | |
| 273 SPARSE_SA2_OP_1 (-) | |
| 274 SPARSE_SA2_OP_2 (*) | |
| 275 SPARSE_SA2_OP_2 (/) | |
| 276 | |
| 277 // Element by element MSparse by MSparse ops. | |
| 278 | |
| 279 #define SPARSE_A2A2_OP(OP) \ | |
| 280 template <class T> \ | |
| 281 MSparse<T> \ | |
| 282 operator OP (const MSparse<T>& a, const MSparse<T>& b) \ | |
| 283 { \ | |
| 284 MSparse<T> r; \ | |
| 285 \ | |
| 5275 | 286 octave_idx_type a_nr = a.rows (); \ |
| 287 octave_idx_type a_nc = a.cols (); \ | |
| 5164 | 288 \ |
| 5275 | 289 octave_idx_type b_nr = b.rows (); \ |
| 290 octave_idx_type b_nc = b.cols (); \ | |
| 5164 | 291 \ |
| 6221 | 292 if (a_nr == 1 && a_nc == 1) \ |
| 293 { \ | |
| 294 if (a.elem(0,0) == 0.) \ | |
| 7342 | 295 r = OP MSparse<T> (b); \ |
| 6221 | 296 else \ |
| 297 { \ | |
|
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
298 r = MSparse<T> (b_nr, b_nc, a.data(0) OP 0.); \ |
| 6221 | 299 \ |
| 300 for (octave_idx_type j = 0 ; j < b_nc ; j++) \ | |
| 301 { \ | |
|
10142
829e69ec3110
make OCTAVE_QUIT a function
Jaroslav Hajek <highegg@gmail.com>
parents:
7344
diff
changeset
|
302 octave_quit (); \ |
| 6221 | 303 octave_idx_type idxj = j * b_nr; \ |
| 304 for (octave_idx_type i = b.cidx(j) ; i < b.cidx(j+1) ; i++) \ | |
| 305 { \ | |
|
10142
829e69ec3110
make OCTAVE_QUIT a function
Jaroslav Hajek <highegg@gmail.com>
parents:
7344
diff
changeset
|
306 octave_quit (); \ |
| 6221 | 307 r.data(idxj + b.ridx(i)) = a.data(0) OP b.data(i); \ |
|
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
308 } \ |
| 6221 | 309 } \ |
| 310 r.maybe_compress (); \ | |
| 311 } \ | |
| 312 } \ | |
| 313 else if (b_nr == 1 && b_nc == 1) \ | |
| 314 { \ | |
| 315 if (b.elem(0,0) == 0.) \ | |
| 316 r = MSparse<T> (a); \ | |
| 317 else \ | |
| 318 { \ | |
|
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
319 r = MSparse<T> (a_nr, a_nc, 0. OP b.data(0)); \ |
| 6221 | 320 \ |
| 321 for (octave_idx_type j = 0 ; j < a_nc ; j++) \ | |
| 322 { \ | |
|
10142
829e69ec3110
make OCTAVE_QUIT a function
Jaroslav Hajek <highegg@gmail.com>
parents:
7344
diff
changeset
|
323 octave_quit (); \ |
| 6221 | 324 octave_idx_type idxj = j * a_nr; \ |
| 325 for (octave_idx_type i = a.cidx(j) ; i < a.cidx(j+1) ; i++) \ | |
| 326 { \ | |
|
10142
829e69ec3110
make OCTAVE_QUIT a function
Jaroslav Hajek <highegg@gmail.com>
parents:
7344
diff
changeset
|
327 octave_quit (); \ |
| 6221 | 328 r.data(idxj + a.ridx(i)) = a.data(i) OP b.data(0); \ |
|
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
329 } \ |
| 6221 | 330 } \ |
| 331 r.maybe_compress (); \ | |
| 332 } \ | |
| 333 } \ | |
| 334 else if (a_nr != b_nr || a_nc != b_nc) \ | |
| 5164 | 335 gripe_nonconformant ("operator " # OP, a_nr, a_nc, b_nr, b_nc); \ |
| 336 else \ | |
| 337 { \ | |
| 5681 | 338 r = MSparse<T> (a_nr, a_nc, (a.nnz () + b.nnz ())); \ |
| 5164 | 339 \ |
| 5275 | 340 octave_idx_type jx = 0; \ |
|
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
341 r.cidx (0) = 0; \ |
| 5275 | 342 for (octave_idx_type i = 0 ; i < a_nc ; i++) \ |
| 5164 | 343 { \ |
| 5275 | 344 octave_idx_type ja = a.cidx(i); \ |
| 345 octave_idx_type ja_max = a.cidx(i+1); \ | |
| 5164 | 346 bool ja_lt_max= ja < ja_max; \ |
| 347 \ | |
| 5275 | 348 octave_idx_type jb = b.cidx(i); \ |
| 349 octave_idx_type jb_max = b.cidx(i+1); \ | |
| 5164 | 350 bool jb_lt_max = jb < jb_max; \ |
| 351 \ | |
| 352 while (ja_lt_max || jb_lt_max ) \ | |
| 353 { \ | |
|
10142
829e69ec3110
make OCTAVE_QUIT a function
Jaroslav Hajek <highegg@gmail.com>
parents:
7344
diff
changeset
|
354 octave_quit (); \ |
| 5164 | 355 if ((! jb_lt_max) || \ |
| 356 (ja_lt_max && (a.ridx(ja) < b.ridx(jb)))) \ | |
| 357 { \ | |
| 358 r.ridx(jx) = a.ridx(ja); \ | |
| 359 r.data(jx) = a.data(ja) OP 0.; \ | |
| 360 jx++; \ | |
| 361 ja++; \ | |
| 362 ja_lt_max= ja < ja_max; \ | |
| 363 } \ | |
| 364 else if (( !ja_lt_max ) || \ | |
| 365 (jb_lt_max && (b.ridx(jb) < a.ridx(ja)) ) ) \ | |
| 366 { \ | |
|
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
367 r.ridx(jx) = b.ridx(jb); \ |
|
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
368 r.data(jx) = 0. OP b.data(jb); \ |
|
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
369 jx++; \ |
| 5164 | 370 jb++; \ |
| 371 jb_lt_max= jb < jb_max; \ | |
| 372 } \ | |
| 373 else \ | |
| 374 { \ | |
|
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
375 if ((a.data(ja) OP b.data(jb)) != 0.) \ |
|
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
376 { \ |
| 5164 | 377 r.data(jx) = a.data(ja) OP b.data(jb); \ |
| 378 r.ridx(jx) = a.ridx(ja); \ | |
| 379 jx++; \ | |
| 380 } \ | |
| 381 ja++; \ | |
| 382 ja_lt_max= ja < ja_max; \ | |
| 383 jb++; \ | |
| 384 jb_lt_max= jb < jb_max; \ | |
| 385 } \ | |
| 386 } \ | |
| 387 r.cidx(i+1) = jx; \ | |
| 388 } \ | |
| 389 \ | |
|
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
390 r.maybe_compress (); \ |
| 5164 | 391 } \ |
| 392 \ | |
| 393 return r; \ | |
| 394 } | |
| 395 | |
|
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
396 #define SPARSE_A2A2_FCN_1(FCN, OP) \ |
| 5164 | 397 template <class T> \ |
| 398 MSparse<T> \ | |
| 399 FCN (const MSparse<T>& a, const MSparse<T>& b) \ | |
| 400 { \ | |
| 401 MSparse<T> r; \ | |
| 402 \ | |
| 5275 | 403 octave_idx_type a_nr = a.rows (); \ |
| 404 octave_idx_type a_nc = a.cols (); \ | |
| 5164 | 405 \ |
| 5275 | 406 octave_idx_type b_nr = b.rows (); \ |
| 407 octave_idx_type b_nc = b.cols (); \ | |
| 5164 | 408 \ |
| 6221 | 409 if (a_nr == 1 && a_nc == 1) \ |
| 410 { \ | |
| 411 if (a.elem(0,0) == 0.) \ | |
| 412 r = MSparse<T> (b_nr, b_nc); \ | |
| 413 else \ | |
| 414 { \ | |
|
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
415 r = MSparse<T> (b); \ |
| 6221 | 416 octave_idx_type b_nnz = b.nnz(); \ |
| 417 \ | |
| 418 for (octave_idx_type i = 0 ; i < b_nnz ; i++) \ | |
| 419 { \ | |
|
10142
829e69ec3110
make OCTAVE_QUIT a function
Jaroslav Hajek <highegg@gmail.com>
parents:
7344
diff
changeset
|
420 octave_quit (); \ |
| 6221 | 421 r.data (i) = a.data(0) OP r.data(i); \ |
| 422 } \ | |
| 423 r.maybe_compress (); \ | |
| 424 } \ | |
| 425 } \ | |
| 426 else if (b_nr == 1 && b_nc == 1) \ | |
| 427 { \ | |
| 428 if (b.elem(0,0) == 0.) \ | |
| 429 r = MSparse<T> (a_nr, a_nc); \ | |
| 430 else \ | |
| 431 { \ | |
|
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
432 r = MSparse<T> (a); \ |
| 6221 | 433 octave_idx_type a_nnz = a.nnz(); \ |
| 434 \ | |
| 435 for (octave_idx_type i = 0 ; i < a_nnz ; i++) \ | |
| 436 { \ | |
|
10142
829e69ec3110
make OCTAVE_QUIT a function
Jaroslav Hajek <highegg@gmail.com>
parents:
7344
diff
changeset
|
437 octave_quit (); \ |
| 6221 | 438 r.data (i) = r.data(i) OP b.data(0); \ |
| 439 } \ | |
| 440 r.maybe_compress (); \ | |
| 441 } \ | |
| 442 } \ | |
| 443 else if (a_nr != b_nr || a_nc != b_nc) \ | |
| 5164 | 444 gripe_nonconformant (#FCN, a_nr, a_nc, b_nr, b_nc); \ |
| 445 else \ | |
| 446 { \ | |
| 5681 | 447 r = MSparse<T> (a_nr, a_nc, (a.nnz () > b.nnz () ? a.nnz () : b.nnz ())); \ |
| 5164 | 448 \ |
| 5275 | 449 octave_idx_type jx = 0; \ |
|
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
450 r.cidx (0) = 0; \ |
| 5275 | 451 for (octave_idx_type i = 0 ; i < a_nc ; i++) \ |
| 5164 | 452 { \ |
| 5275 | 453 octave_idx_type ja = a.cidx(i); \ |
| 454 octave_idx_type ja_max = a.cidx(i+1); \ | |
| 5164 | 455 bool ja_lt_max= ja < ja_max; \ |
| 456 \ | |
| 5275 | 457 octave_idx_type jb = b.cidx(i); \ |
| 458 octave_idx_type jb_max = b.cidx(i+1); \ | |
| 5164 | 459 bool jb_lt_max = jb < jb_max; \ |
| 460 \ | |
| 461 while (ja_lt_max || jb_lt_max ) \ | |
| 462 { \ | |
|
10142
829e69ec3110
make OCTAVE_QUIT a function
Jaroslav Hajek <highegg@gmail.com>
parents:
7344
diff
changeset
|
463 octave_quit (); \ |
| 5164 | 464 if ((! jb_lt_max) || \ |
| 465 (ja_lt_max && (a.ridx(ja) < b.ridx(jb)))) \ | |
| 466 { \ | |
| 467 ja++; ja_lt_max= ja < ja_max; \ | |
| 468 } \ | |
| 469 else if (( !ja_lt_max ) || \ | |
| 470 (jb_lt_max && (b.ridx(jb) < a.ridx(ja)) ) ) \ | |
| 471 { \ | |
| 472 jb++; jb_lt_max= jb < jb_max; \ | |
| 473 } \ | |
| 474 else \ | |
| 475 { \ | |
|
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
476 if ((a.data(ja) OP b.data(jb)) != 0.) \ |
|
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
477 { \ |
| 5164 | 478 r.data(jx) = a.data(ja) OP b.data(jb); \ |
| 479 r.ridx(jx) = a.ridx(ja); \ | |
| 480 jx++; \ | |
| 481 } \ | |
| 482 ja++; ja_lt_max= ja < ja_max; \ | |
| 483 jb++; jb_lt_max= jb < jb_max; \ | |
| 484 } \ | |
| 485 } \ | |
| 486 r.cidx(i+1) = jx; \ | |
| 487 } \ | |
| 488 \ | |
|
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
489 r.maybe_compress (); \ |
| 5164 | 490 } \ |
| 491 \ | |
| 492 return r; \ | |
| 493 } | |
| 494 | |
|
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
495 #define SPARSE_A2A2_FCN_2(FCN, OP) \ |
| 5164 | 496 template <class T> \ |
| 497 MSparse<T> \ | |
| 498 FCN (const MSparse<T>& a, const MSparse<T>& b) \ | |
| 499 { \ | |
| 500 MSparse<T> r; \ | |
| 501 T Zero = T (); \ | |
| 502 \ | |
| 5275 | 503 octave_idx_type a_nr = a.rows (); \ |
| 504 octave_idx_type a_nc = a.cols (); \ | |
| 5164 | 505 \ |
| 5275 | 506 octave_idx_type b_nr = b.rows (); \ |
| 507 octave_idx_type b_nc = b.cols (); \ | |
| 5164 | 508 \ |
| 6221 | 509 if (a_nr == 1 && a_nc == 1) \ |
| 510 { \ | |
| 511 T val = a.elem (0,0); \ | |
| 512 T fill = val OP T(); \ | |
| 513 if (fill == T()) \ | |
| 514 { \ | |
| 515 octave_idx_type b_nnz = b.nnz(); \ | |
| 516 r = MSparse<T> (b); \ | |
| 517 for (octave_idx_type i = 0 ; i < b_nnz ; i++) \ | |
| 518 r.data (i) = val OP r.data(i); \ | |
| 519 r.maybe_compress (); \ | |
| 520 } \ | |
| 521 else \ | |
| 522 { \ | |
| 523 r = MSparse<T> (b_nr, b_nc, fill); \ | |
| 524 for (octave_idx_type j = 0 ; j < b_nc ; j++) \ | |
| 525 { \ | |
|
10142
829e69ec3110
make OCTAVE_QUIT a function
Jaroslav Hajek <highegg@gmail.com>
parents:
7344
diff
changeset
|
526 octave_quit (); \ |
| 6221 | 527 octave_idx_type idxj = j * b_nr; \ |
| 528 for (octave_idx_type i = b.cidx(j) ; i < b.cidx(j+1) ; i++) \ | |
| 529 { \ | |
|
10142
829e69ec3110
make OCTAVE_QUIT a function
Jaroslav Hajek <highegg@gmail.com>
parents:
7344
diff
changeset
|
530 octave_quit (); \ |
| 6221 | 531 r.data(idxj + b.ridx(i)) = val OP b.data(i); \ |
|
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
532 } \ |
| 6221 | 533 } \ |
| 534 r.maybe_compress (); \ | |
| 535 } \ | |
| 536 } \ | |
| 537 else if (b_nr == 1 && b_nc == 1) \ | |
| 538 { \ | |
| 539 T val = b.elem (0,0); \ | |
| 540 T fill = T() OP val; \ | |
| 541 if (fill == T()) \ | |
| 542 { \ | |
| 543 octave_idx_type a_nnz = a.nnz(); \ | |
| 544 r = MSparse<T> (a); \ | |
| 545 for (octave_idx_type i = 0 ; i < a_nnz ; i++) \ | |
| 546 r.data (i) = r.data(i) OP val; \ | |
| 547 r.maybe_compress (); \ | |
| 548 } \ | |
| 549 else \ | |
| 550 { \ | |
| 551 r = MSparse<T> (a_nr, a_nc, fill); \ | |
| 552 for (octave_idx_type j = 0 ; j < a_nc ; j++) \ | |
| 553 { \ | |
|
10142
829e69ec3110
make OCTAVE_QUIT a function
Jaroslav Hajek <highegg@gmail.com>
parents:
7344
diff
changeset
|
554 octave_quit (); \ |
| 6221 | 555 octave_idx_type idxj = j * a_nr; \ |
| 556 for (octave_idx_type i = a.cidx(j) ; i < a.cidx(j+1) ; i++) \ | |
| 557 { \ | |
|
10142
829e69ec3110
make OCTAVE_QUIT a function
Jaroslav Hajek <highegg@gmail.com>
parents:
7344
diff
changeset
|
558 octave_quit (); \ |
| 6221 | 559 r.data(idxj + a.ridx(i)) = a.data(i) OP val; \ |
|
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
560 } \ |
| 6221 | 561 } \ |
| 562 r.maybe_compress (); \ | |
| 563 } \ | |
| 564 } \ | |
| 565 else if (a_nr != b_nr || a_nc != b_nc) \ | |
| 5164 | 566 gripe_nonconformant (#FCN, a_nr, a_nc, b_nr, b_nc); \ |
| 567 else \ | |
| 568 { \ | |
| 569 r = MSparse<T>( a_nr, a_nc, (Zero OP Zero)); \ | |
| 570 \ | |
| 5275 | 571 for (octave_idx_type i = 0 ; i < a_nc ; i++) \ |
| 5164 | 572 { \ |
| 5275 | 573 octave_idx_type ja = a.cidx(i); \ |
| 574 octave_idx_type ja_max = a.cidx(i+1); \ | |
| 5164 | 575 bool ja_lt_max= ja < ja_max; \ |
| 576 \ | |
| 5275 | 577 octave_idx_type jb = b.cidx(i); \ |
| 578 octave_idx_type jb_max = b.cidx(i+1); \ | |
| 5164 | 579 bool jb_lt_max = jb < jb_max; \ |
| 580 \ | |
| 581 while (ja_lt_max || jb_lt_max ) \ | |
| 582 { \ | |
|
10142
829e69ec3110
make OCTAVE_QUIT a function
Jaroslav Hajek <highegg@gmail.com>
parents:
7344
diff
changeset
|
583 octave_quit (); \ |
| 5164 | 584 if ((! jb_lt_max) || \ |
| 585 (ja_lt_max && (a.ridx(ja) < b.ridx(jb)))) \ | |
| 586 { \ | |
|
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
587 r.elem (a.ridx(ja),i) = a.data(ja) OP Zero; \ |
| 5164 | 588 ja++; ja_lt_max= ja < ja_max; \ |
| 589 } \ | |
| 590 else if (( !ja_lt_max ) || \ | |
| 591 (jb_lt_max && (b.ridx(jb) < a.ridx(ja)) ) ) \ | |
| 592 { \ | |
|
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
593 r.elem (b.ridx(jb),i) = Zero OP b.data(jb); \ |
| 5164 | 594 jb++; jb_lt_max= jb < jb_max; \ |
| 595 } \ | |
| 596 else \ | |
| 597 { \ | |
| 598 r.elem (a.ridx(ja),i) = a.data(ja) OP b.data(jb); \ | |
| 599 ja++; ja_lt_max= ja < ja_max; \ | |
| 600 jb++; jb_lt_max= jb < jb_max; \ | |
| 601 } \ | |
| 602 } \ | |
| 603 } \ | |
| 604 \ | |
|
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
605 r.maybe_compress (true); \ |
| 5164 | 606 } \ |
| 607 \ | |
| 608 return r; \ | |
| 609 } | |
| 610 | |
| 611 SPARSE_A2A2_OP (+) | |
| 612 SPARSE_A2A2_OP (-) | |
| 613 SPARSE_A2A2_FCN_1 (product, *) | |
| 614 SPARSE_A2A2_FCN_2 (quotient, /) | |
| 615 | |
| 616 // Unary MSparse ops. | |
| 617 | |
| 618 template <class T> | |
| 619 MSparse<T> | |
| 620 operator + (const MSparse<T>& a) | |
| 621 { | |
| 622 return a; | |
| 623 } | |
| 624 | |
| 625 template <class T> | |
| 626 MSparse<T> | |
| 627 operator - (const MSparse<T>& a) | |
| 628 { | |
| 629 MSparse<T> retval (a); | |
| 5681 | 630 octave_idx_type nz = a.nnz (); |
| 5275 | 631 for (octave_idx_type i = 0; i < nz; i++) |
| 5164 | 632 retval.data(i) = - retval.data(i); |
| 633 return retval; | |
| 634 } |