5164
|
1 /* |
|
2 |
|
3 Copyright (C) 2004 David Bateman |
|
4 Copyright (C) 1998-2004 Andy Adler |
|
5 |
|
6 Octave is free software; you can redistribute it and/or modify it |
|
7 under the terms of the GNU General Public License as published by the |
|
8 Free Software Foundation; either version 2, or (at your option) any |
|
9 later version. |
|
10 |
|
11 Octave is distributed in the hope that it will be useful, but WITHOUT |
|
12 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or |
|
13 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License |
|
14 for more details. |
|
15 |
|
16 You should have received a copy of the GNU General Public License |
|
17 along with this program; see the file COPYING. If not, write to the Free |
|
18 Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. |
|
19 |
|
20 */ |
|
21 |
|
22 #ifdef HAVE_CONFIG_H |
|
23 #include <config.h> |
|
24 #endif |
|
25 |
|
26 #include "quit.h" |
|
27 #include "lo-error.h" |
|
28 #include "MArray2.h" |
|
29 #include "Array-util.h" |
|
30 |
|
31 #include "MSparse.h" |
|
32 #include "MSparse-defs.h" |
|
33 |
|
34 // sparse array with math ops. |
|
35 |
|
36 // Element by element MSparse by MSparse ops. |
|
37 |
|
38 template <class T> |
|
39 MSparse<T>& |
|
40 operator += (MSparse<T>& a, const MSparse<T>& b) |
|
41 { |
|
42 MSparse<T> r; |
|
43 |
5275
|
44 octave_idx_type a_nr = a.rows (); |
|
45 octave_idx_type a_nc = a.cols (); |
5164
|
46 |
5275
|
47 octave_idx_type b_nr = b.rows (); |
|
48 octave_idx_type b_nc = b.cols (); |
5164
|
49 |
|
50 if (a_nr != b_nr || a_nc != b_nc) |
|
51 gripe_nonconformant ("operator +=" , a_nr, a_nc, b_nr, b_nc); |
|
52 else |
|
53 { |
|
54 r = MSparse<T> (a_nr, a_nc, (a.nnz () + b.nnz ())); |
|
55 |
5275
|
56 octave_idx_type jx = 0; |
|
57 for (octave_idx_type i = 0 ; i < a_nc ; i++) |
5164
|
58 { |
5275
|
59 octave_idx_type ja = a.cidx(i); |
|
60 octave_idx_type ja_max = a.cidx(i+1); |
5164
|
61 bool ja_lt_max= ja < ja_max; |
|
62 |
5275
|
63 octave_idx_type jb = b.cidx(i); |
|
64 octave_idx_type jb_max = b.cidx(i+1); |
5164
|
65 bool jb_lt_max = jb < jb_max; |
|
66 |
|
67 while (ja_lt_max || jb_lt_max ) |
|
68 { |
|
69 OCTAVE_QUIT; |
|
70 if ((! jb_lt_max) || |
|
71 (ja_lt_max && (a.ridx(ja) < b.ridx(jb)))) |
|
72 { |
|
73 r.ridx(jx) = a.ridx(ja); |
|
74 r.data(jx) = a.data(ja) + 0.; |
|
75 jx++; |
|
76 ja++; |
|
77 ja_lt_max= ja < ja_max; |
|
78 } |
|
79 else if (( !ja_lt_max ) || |
|
80 (jb_lt_max && (b.ridx(jb) < a.ridx(ja)) ) ) |
|
81 { |
|
82 r.ridx(jx) = b.ridx(jb); |
|
83 r.data(jx) = 0. + b.data(jb); |
|
84 jx++; |
|
85 jb++; |
|
86 jb_lt_max= jb < jb_max; |
|
87 } |
|
88 else |
|
89 { |
|
90 if ((a.data(ja) + b.data(jb)) != 0.) |
|
91 { |
|
92 r.data(jx) = a.data(ja) + b.data(jb); |
|
93 r.ridx(jx) = a.ridx(ja); |
|
94 jx++; |
|
95 } |
|
96 ja++; |
|
97 ja_lt_max= ja < ja_max; |
|
98 jb++; |
|
99 jb_lt_max= jb < jb_max; |
|
100 } |
|
101 } |
|
102 r.cidx(i+1) = jx; |
|
103 } |
|
104 |
|
105 a = r.maybe_compress (); |
|
106 } |
|
107 |
|
108 return a; |
|
109 } |
|
110 |
|
111 template <class T> |
|
112 MSparse<T>& |
|
113 operator -= (MSparse<T>& a, const MSparse<T>& b) |
|
114 { |
|
115 MSparse<T> r; |
|
116 |
5275
|
117 octave_idx_type a_nr = a.rows (); |
|
118 octave_idx_type a_nc = a.cols (); |
5164
|
119 |
5275
|
120 octave_idx_type b_nr = b.rows (); |
|
121 octave_idx_type b_nc = b.cols (); |
5164
|
122 |
|
123 if (a_nr != b_nr || a_nc != b_nc) |
|
124 gripe_nonconformant ("operator -=" , a_nr, a_nc, b_nr, b_nc); |
|
125 else |
|
126 { |
|
127 r = MSparse<T> (a_nr, a_nc, (a.nnz () + b.nnz ())); |
|
128 |
5275
|
129 octave_idx_type jx = 0; |
|
130 for (octave_idx_type i = 0 ; i < a_nc ; i++) |
5164
|
131 { |
5275
|
132 octave_idx_type ja = a.cidx(i); |
|
133 octave_idx_type ja_max = a.cidx(i+1); |
5164
|
134 bool ja_lt_max= ja < ja_max; |
|
135 |
5275
|
136 octave_idx_type jb = b.cidx(i); |
|
137 octave_idx_type jb_max = b.cidx(i+1); |
5164
|
138 bool jb_lt_max = jb < jb_max; |
|
139 |
|
140 while (ja_lt_max || jb_lt_max ) |
|
141 { |
|
142 OCTAVE_QUIT; |
|
143 if ((! jb_lt_max) || |
|
144 (ja_lt_max && (a.ridx(ja) < b.ridx(jb)))) |
|
145 { |
|
146 r.ridx(jx) = a.ridx(ja); |
|
147 r.data(jx) = a.data(ja) - 0.; |
|
148 jx++; |
|
149 ja++; |
|
150 ja_lt_max= ja < ja_max; |
|
151 } |
|
152 else if (( !ja_lt_max ) || |
|
153 (jb_lt_max && (b.ridx(jb) < a.ridx(ja)) ) ) |
|
154 { |
|
155 r.ridx(jx) = b.ridx(jb); |
|
156 r.data(jx) = 0. - b.data(jb); |
|
157 jx++; |
|
158 jb++; |
|
159 jb_lt_max= jb < jb_max; |
|
160 } |
|
161 else |
|
162 { |
|
163 if ((a.data(ja) - b.data(jb)) != 0.) |
|
164 { |
|
165 r.data(jx) = a.data(ja) - b.data(jb); |
|
166 r.ridx(jx) = a.ridx(ja); |
|
167 jx++; |
|
168 } |
|
169 ja++; |
|
170 ja_lt_max= ja < ja_max; |
|
171 jb++; |
|
172 jb_lt_max= jb < jb_max; |
|
173 } |
|
174 } |
|
175 r.cidx(i+1) = jx; |
|
176 } |
|
177 |
|
178 a = r.maybe_compress (); |
|
179 } |
|
180 |
|
181 return a; |
|
182 } |
|
183 |
|
184 // Element by element MSparse by scalar ops. |
|
185 |
|
186 #define SPARSE_A2S_OP_1(OP) \ |
|
187 template <class T> \ |
|
188 MArray2<T> \ |
|
189 operator OP (const MSparse<T>& a, const T& s) \ |
|
190 { \ |
5275
|
191 octave_idx_type nr = a.rows (); \ |
|
192 octave_idx_type nc = a.cols (); \ |
5164
|
193 \ |
|
194 MArray2<T> r (nr, nc, (0.0 OP s)); \ |
|
195 \ |
5275
|
196 for (octave_idx_type j = 0; j < nc; j++) \ |
|
197 for (octave_idx_type i = a.cidx(j); i < a.cidx(j+1); i++) \ |
5164
|
198 r.elem (a.ridx (i), j) = a.data (i) OP s; \ |
|
199 return r; \ |
|
200 } |
|
201 |
|
202 #define SPARSE_A2S_OP_2(OP) \ |
|
203 template <class T> \ |
|
204 MSparse<T> \ |
|
205 operator OP (const MSparse<T>& a, const T& s) \ |
|
206 { \ |
5275
|
207 octave_idx_type nr = a.rows (); \ |
|
208 octave_idx_type nc = a.cols (); \ |
|
209 octave_idx_type nz = a.nnz (); \ |
5164
|
210 \ |
|
211 MSparse<T> r (nr, nc, nz); \ |
|
212 \ |
5275
|
213 for (octave_idx_type i = 0; i < nz; i++) \ |
5164
|
214 { \ |
|
215 r.data(i) = a.data(i) OP s; \ |
|
216 r.ridx(i) = a.ridx(i); \ |
|
217 } \ |
5275
|
218 for (octave_idx_type i = 0; i < nc + 1; i++) \ |
5164
|
219 r.cidx(i) = a.cidx(i); \ |
|
220 r.maybe_compress (true); \ |
|
221 return r; \ |
|
222 } |
|
223 |
|
224 |
|
225 SPARSE_A2S_OP_1 (+) |
|
226 SPARSE_A2S_OP_1 (-) |
|
227 SPARSE_A2S_OP_2 (*) |
|
228 SPARSE_A2S_OP_2 (/) |
|
229 |
|
230 // Element by element scalar by MSparse ops. |
|
231 |
|
232 #define SPARSE_SA2_OP_1(OP) \ |
|
233 template <class T> \ |
|
234 MArray2<T> \ |
|
235 operator OP (const T& s, const MSparse<T>& a) \ |
|
236 { \ |
5275
|
237 octave_idx_type nr = a.rows (); \ |
|
238 octave_idx_type nc = a.cols (); \ |
5164
|
239 \ |
|
240 MArray2<T> r (nr, nc, (s OP 0.0)); \ |
|
241 \ |
5275
|
242 for (octave_idx_type j = 0; j < nc; j++) \ |
|
243 for (octave_idx_type i = a.cidx(j); i < a.cidx(j+1); i++) \ |
5164
|
244 r.elem (a.ridx (i), j) = s OP a.data (i); \ |
|
245 return r; \ |
|
246 } |
|
247 |
|
248 #define SPARSE_SA2_OP_2(OP) \ |
|
249 template <class T> \ |
|
250 MSparse<T> \ |
|
251 operator OP (const T& s, const MSparse<T>& a) \ |
|
252 { \ |
5275
|
253 octave_idx_type nr = a.rows (); \ |
|
254 octave_idx_type nc = a.cols (); \ |
|
255 octave_idx_type nz = a.nnz (); \ |
5164
|
256 \ |
|
257 MSparse<T> r (nr, nc, nz); \ |
|
258 \ |
5275
|
259 for (octave_idx_type i = 0; i < nz; i++) \ |
5164
|
260 { \ |
|
261 r.data(i) = s OP a.data(i); \ |
|
262 r.ridx(i) = a.ridx(i); \ |
|
263 } \ |
5275
|
264 for (octave_idx_type i = 0; i < nc + 1; i++) \ |
5164
|
265 r.cidx(i) = a.cidx(i); \ |
|
266 r.maybe_compress (true); \ |
|
267 return r; \ |
|
268 } |
|
269 |
|
270 SPARSE_SA2_OP_1 (+) |
|
271 SPARSE_SA2_OP_1 (-) |
|
272 SPARSE_SA2_OP_2 (*) |
|
273 SPARSE_SA2_OP_2 (/) |
|
274 |
|
275 // Element by element MSparse by MSparse ops. |
|
276 |
|
277 #define SPARSE_A2A2_OP(OP) \ |
|
278 template <class T> \ |
|
279 MSparse<T> \ |
|
280 operator OP (const MSparse<T>& a, const MSparse<T>& b) \ |
|
281 { \ |
|
282 MSparse<T> r; \ |
|
283 \ |
5275
|
284 octave_idx_type a_nr = a.rows (); \ |
|
285 octave_idx_type a_nc = a.cols (); \ |
5164
|
286 \ |
5275
|
287 octave_idx_type b_nr = b.rows (); \ |
|
288 octave_idx_type b_nc = b.cols (); \ |
5164
|
289 \ |
|
290 if (a_nr != b_nr || a_nc != b_nc) \ |
|
291 gripe_nonconformant ("operator " # OP, a_nr, a_nc, b_nr, b_nc); \ |
|
292 else \ |
|
293 { \ |
|
294 r = MSparse<T> (a_nr, a_nc, (a.nnz () + b.nnz ())); \ |
|
295 \ |
5275
|
296 octave_idx_type jx = 0; \ |
5164
|
297 r.cidx (0) = 0; \ |
5275
|
298 for (octave_idx_type i = 0 ; i < a_nc ; i++) \ |
5164
|
299 { \ |
5275
|
300 octave_idx_type ja = a.cidx(i); \ |
|
301 octave_idx_type ja_max = a.cidx(i+1); \ |
5164
|
302 bool ja_lt_max= ja < ja_max; \ |
|
303 \ |
5275
|
304 octave_idx_type jb = b.cidx(i); \ |
|
305 octave_idx_type jb_max = b.cidx(i+1); \ |
5164
|
306 bool jb_lt_max = jb < jb_max; \ |
|
307 \ |
|
308 while (ja_lt_max || jb_lt_max ) \ |
|
309 { \ |
|
310 OCTAVE_QUIT; \ |
|
311 if ((! jb_lt_max) || \ |
|
312 (ja_lt_max && (a.ridx(ja) < b.ridx(jb)))) \ |
|
313 { \ |
|
314 r.ridx(jx) = a.ridx(ja); \ |
|
315 r.data(jx) = a.data(ja) OP 0.; \ |
|
316 jx++; \ |
|
317 ja++; \ |
|
318 ja_lt_max= ja < ja_max; \ |
|
319 } \ |
|
320 else if (( !ja_lt_max ) || \ |
|
321 (jb_lt_max && (b.ridx(jb) < a.ridx(ja)) ) ) \ |
|
322 { \ |
|
323 r.ridx(jx) = b.ridx(jb); \ |
|
324 r.data(jx) = 0. OP b.data(jb); \ |
|
325 jx++; \ |
|
326 jb++; \ |
|
327 jb_lt_max= jb < jb_max; \ |
|
328 } \ |
|
329 else \ |
|
330 { \ |
|
331 if ((a.data(ja) OP b.data(jb)) != 0.) \ |
|
332 { \ |
|
333 r.data(jx) = a.data(ja) OP b.data(jb); \ |
|
334 r.ridx(jx) = a.ridx(ja); \ |
|
335 jx++; \ |
|
336 } \ |
|
337 ja++; \ |
|
338 ja_lt_max= ja < ja_max; \ |
|
339 jb++; \ |
|
340 jb_lt_max= jb < jb_max; \ |
|
341 } \ |
|
342 } \ |
|
343 r.cidx(i+1) = jx; \ |
|
344 } \ |
|
345 \ |
|
346 r.maybe_compress (); \ |
|
347 } \ |
|
348 \ |
|
349 return r; \ |
|
350 } |
|
351 |
|
352 #define SPARSE_A2A2_FCN_1(FCN, OP) \ |
|
353 template <class T> \ |
|
354 MSparse<T> \ |
|
355 FCN (const MSparse<T>& a, const MSparse<T>& b) \ |
|
356 { \ |
|
357 MSparse<T> r; \ |
|
358 \ |
5275
|
359 octave_idx_type a_nr = a.rows (); \ |
|
360 octave_idx_type a_nc = a.cols (); \ |
5164
|
361 \ |
5275
|
362 octave_idx_type b_nr = b.rows (); \ |
|
363 octave_idx_type b_nc = b.cols (); \ |
5164
|
364 \ |
|
365 if (a_nr != b_nr || a_nc != b_nc) \ |
|
366 gripe_nonconformant (#FCN, a_nr, a_nc, b_nr, b_nc); \ |
|
367 else \ |
|
368 { \ |
|
369 r = MSparse<T> (a_nr, a_nc, (a.nnz() > b.nnz() ? a.nnz() : b.nnz())); \ |
|
370 \ |
5275
|
371 octave_idx_type jx = 0; \ |
5164
|
372 r.cidx (0) = 0; \ |
5275
|
373 for (octave_idx_type i = 0 ; i < a_nc ; i++) \ |
5164
|
374 { \ |
5275
|
375 octave_idx_type ja = a.cidx(i); \ |
|
376 octave_idx_type ja_max = a.cidx(i+1); \ |
5164
|
377 bool ja_lt_max= ja < ja_max; \ |
|
378 \ |
5275
|
379 octave_idx_type jb = b.cidx(i); \ |
|
380 octave_idx_type jb_max = b.cidx(i+1); \ |
5164
|
381 bool jb_lt_max = jb < jb_max; \ |
|
382 \ |
|
383 while (ja_lt_max || jb_lt_max ) \ |
|
384 { \ |
|
385 OCTAVE_QUIT; \ |
|
386 if ((! jb_lt_max) || \ |
|
387 (ja_lt_max && (a.ridx(ja) < b.ridx(jb)))) \ |
|
388 { \ |
|
389 ja++; ja_lt_max= ja < ja_max; \ |
|
390 } \ |
|
391 else if (( !ja_lt_max ) || \ |
|
392 (jb_lt_max && (b.ridx(jb) < a.ridx(ja)) ) ) \ |
|
393 { \ |
|
394 jb++; jb_lt_max= jb < jb_max; \ |
|
395 } \ |
|
396 else \ |
|
397 { \ |
|
398 if ((a.data(ja) OP b.data(jb)) != 0.) \ |
|
399 { \ |
|
400 r.data(jx) = a.data(ja) OP b.data(jb); \ |
|
401 r.ridx(jx) = a.ridx(ja); \ |
|
402 jx++; \ |
|
403 } \ |
|
404 ja++; ja_lt_max= ja < ja_max; \ |
|
405 jb++; jb_lt_max= jb < jb_max; \ |
|
406 } \ |
|
407 } \ |
|
408 r.cidx(i+1) = jx; \ |
|
409 } \ |
|
410 \ |
|
411 r.maybe_compress (); \ |
|
412 } \ |
|
413 \ |
|
414 return r; \ |
|
415 } |
|
416 |
|
417 #define SPARSE_A2A2_FCN_2(FCN, OP) \ |
|
418 template <class T> \ |
|
419 MSparse<T> \ |
|
420 FCN (const MSparse<T>& a, const MSparse<T>& b) \ |
|
421 { \ |
|
422 MSparse<T> r; \ |
|
423 T Zero = T (); \ |
|
424 \ |
5275
|
425 octave_idx_type a_nr = a.rows (); \ |
|
426 octave_idx_type a_nc = a.cols (); \ |
5164
|
427 \ |
5275
|
428 octave_idx_type b_nr = b.rows (); \ |
|
429 octave_idx_type b_nc = b.cols (); \ |
5164
|
430 \ |
|
431 if (a_nr != b_nr || a_nc != b_nc) \ |
|
432 gripe_nonconformant (#FCN, a_nr, a_nc, b_nr, b_nc); \ |
|
433 else \ |
|
434 { \ |
|
435 r = MSparse<T>( a_nr, a_nc, (Zero OP Zero)); \ |
|
436 \ |
5275
|
437 for (octave_idx_type i = 0 ; i < a_nc ; i++) \ |
5164
|
438 { \ |
5275
|
439 octave_idx_type ja = a.cidx(i); \ |
|
440 octave_idx_type ja_max = a.cidx(i+1); \ |
5164
|
441 bool ja_lt_max= ja < ja_max; \ |
|
442 \ |
5275
|
443 octave_idx_type jb = b.cidx(i); \ |
|
444 octave_idx_type jb_max = b.cidx(i+1); \ |
5164
|
445 bool jb_lt_max = jb < jb_max; \ |
|
446 \ |
|
447 while (ja_lt_max || jb_lt_max ) \ |
|
448 { \ |
|
449 OCTAVE_QUIT; \ |
|
450 if ((! jb_lt_max) || \ |
|
451 (ja_lt_max && (a.ridx(ja) < b.ridx(jb)))) \ |
|
452 { \ |
|
453 r.elem (a.ridx(ja),i) = a.data(ja) OP Zero; \ |
|
454 ja++; ja_lt_max= ja < ja_max; \ |
|
455 } \ |
|
456 else if (( !ja_lt_max ) || \ |
|
457 (jb_lt_max && (b.ridx(jb) < a.ridx(ja)) ) ) \ |
|
458 { \ |
|
459 r.elem (b.ridx(jb),i) = Zero OP b.data(jb); \ |
|
460 jb++; jb_lt_max= jb < jb_max; \ |
|
461 } \ |
|
462 else \ |
|
463 { \ |
|
464 r.elem (a.ridx(ja),i) = a.data(ja) OP b.data(jb); \ |
|
465 ja++; ja_lt_max= ja < ja_max; \ |
|
466 jb++; jb_lt_max= jb < jb_max; \ |
|
467 } \ |
|
468 } \ |
|
469 } \ |
|
470 \ |
|
471 r.maybe_compress (true); \ |
|
472 } \ |
|
473 \ |
|
474 return r; \ |
|
475 } |
|
476 |
|
477 SPARSE_A2A2_OP (+) |
|
478 SPARSE_A2A2_OP (-) |
|
479 SPARSE_A2A2_FCN_1 (product, *) |
|
480 SPARSE_A2A2_FCN_2 (quotient, /) |
|
481 |
|
482 // Unary MSparse ops. |
|
483 |
|
484 template <class T> |
|
485 MSparse<T> |
|
486 operator + (const MSparse<T>& a) |
|
487 { |
|
488 return a; |
|
489 } |
|
490 |
|
491 template <class T> |
|
492 MSparse<T> |
|
493 operator - (const MSparse<T>& a) |
|
494 { |
|
495 MSparse<T> retval (a); |
5275
|
496 octave_idx_type nz = a.nnz (); |
|
497 for (octave_idx_type i = 0; i < nz; i++) |
5164
|
498 retval.data(i) = - retval.data(i); |
|
499 return retval; |
|
500 } |
|
501 |
|
502 /* |
|
503 ;;; Local Variables: *** |
|
504 ;;; mode: C++ *** |
|
505 ;;; End: *** |
|
506 */ |