|
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 |
|
|
44 int a_nr = a.rows (); |
|
|
45 int a_nc = a.cols (); |
|
|
46 |
|
|
47 int b_nr = b.rows (); |
|
|
48 int b_nc = b.cols (); |
|
|
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 |
|
|
56 int jx = 0; |
|
|
57 for (int i = 0 ; i < a_nc ; i++) |
|
|
58 { |
|
|
59 int ja = a.cidx(i); |
|
|
60 int ja_max = a.cidx(i+1); |
|
|
61 bool ja_lt_max= ja < ja_max; |
|
|
62 |
|
|
63 int jb = b.cidx(i); |
|
|
64 int jb_max = b.cidx(i+1); |
|
|
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 |
|
|
117 int a_nr = a.rows (); |
|
|
118 int a_nc = a.cols (); |
|
|
119 |
|
|
120 int b_nr = b.rows (); |
|
|
121 int b_nc = b.cols (); |
|
|
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 |
|
|
129 int jx = 0; |
|
|
130 for (int i = 0 ; i < a_nc ; i++) |
|
|
131 { |
|
|
132 int ja = a.cidx(i); |
|
|
133 int ja_max = a.cidx(i+1); |
|
|
134 bool ja_lt_max= ja < ja_max; |
|
|
135 |
|
|
136 int jb = b.cidx(i); |
|
|
137 int jb_max = b.cidx(i+1); |
|
|
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 { \ |
|
|
191 int nr = a.rows (); \ |
|
|
192 int nc = a.cols (); \ |
|
|
193 \ |
|
|
194 MArray2<T> r (nr, nc, (0.0 OP s)); \ |
|
|
195 \ |
|
|
196 for (int j = 0; j < nc; j++) \ |
|
|
197 for (int i = a.cidx(j); i < a.cidx(j+1); i++) \ |
|
|
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 { \ |
|
|
207 int nr = a.rows (); \ |
|
|
208 int nc = a.cols (); \ |
|
|
209 int nz = a.nnz (); \ |
|
|
210 \ |
|
|
211 MSparse<T> r (nr, nc, nz); \ |
|
|
212 \ |
|
|
213 for (int i = 0; i < nz; i++) \ |
|
|
214 { \ |
|
|
215 r.data(i) = a.data(i) OP s; \ |
|
|
216 r.ridx(i) = a.ridx(i); \ |
|
|
217 } \ |
|
|
218 for (int i = 0; i < nc + 1; i++) \ |
|
|
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 { \ |
|
|
237 int nr = a.rows (); \ |
|
|
238 int nc = a.cols (); \ |
|
|
239 \ |
|
|
240 MArray2<T> r (nr, nc, (s OP 0.0)); \ |
|
|
241 \ |
|
|
242 for (int j = 0; j < nc; j++) \ |
|
|
243 for (int i = a.cidx(j); i < a.cidx(j+1); i++) \ |
|
|
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 { \ |
|
|
253 int nr = a.rows (); \ |
|
|
254 int nc = a.cols (); \ |
|
|
255 int nz = a.nnz (); \ |
|
|
256 \ |
|
|
257 MSparse<T> r (nr, nc, nz); \ |
|
|
258 \ |
|
|
259 for (int i = 0; i < nz; i++) \ |
|
|
260 { \ |
|
|
261 r.data(i) = s OP a.data(i); \ |
|
|
262 r.ridx(i) = a.ridx(i); \ |
|
|
263 } \ |
|
|
264 for (int i = 0; i < nc + 1; i++) \ |
|
|
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 \ |
|
|
284 int a_nr = a.rows (); \ |
|
|
285 int a_nc = a.cols (); \ |
|
|
286 \ |
|
|
287 int b_nr = b.rows (); \ |
|
|
288 int b_nc = b.cols (); \ |
|
|
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 \ |
|
|
296 int jx = 0; \ |
|
|
297 r.cidx (0) = 0; \ |
|
|
298 for (int i = 0 ; i < a_nc ; i++) \ |
|
|
299 { \ |
|
|
300 int ja = a.cidx(i); \ |
|
|
301 int ja_max = a.cidx(i+1); \ |
|
|
302 bool ja_lt_max= ja < ja_max; \ |
|
|
303 \ |
|
|
304 int jb = b.cidx(i); \ |
|
|
305 int jb_max = b.cidx(i+1); \ |
|
|
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 \ |
|
|
359 int a_nr = a.rows (); \ |
|
|
360 int a_nc = a.cols (); \ |
|
|
361 \ |
|
|
362 int b_nr = b.rows (); \ |
|
|
363 int b_nc = b.cols (); \ |
|
|
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 \ |
|
|
371 int jx = 0; \ |
|
|
372 r.cidx (0) = 0; \ |
|
|
373 for (int i = 0 ; i < a_nc ; i++) \ |
|
|
374 { \ |
|
|
375 int ja = a.cidx(i); \ |
|
|
376 int ja_max = a.cidx(i+1); \ |
|
|
377 bool ja_lt_max= ja < ja_max; \ |
|
|
378 \ |
|
|
379 int jb = b.cidx(i); \ |
|
|
380 int jb_max = b.cidx(i+1); \ |
|
|
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 \ |
|
|
425 int a_nr = a.rows (); \ |
|
|
426 int a_nc = a.cols (); \ |
|
|
427 \ |
|
|
428 int b_nr = b.rows (); \ |
|
|
429 int b_nc = b.cols (); \ |
|
|
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 \ |
|
|
437 for (int i = 0 ; i < a_nc ; i++) \ |
|
|
438 { \ |
|
|
439 int ja = a.cidx(i); \ |
|
|
440 int ja_max = a.cidx(i+1); \ |
|
|
441 bool ja_lt_max= ja < ja_max; \ |
|
|
442 \ |
|
|
443 int jb = b.cidx(i); \ |
|
|
444 int jb_max = b.cidx(i+1); \ |
|
|
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); |
|
|
496 int nz = a.nnz (); |
|
|
497 for (int i = 0; i < nz; i++) |
|
|
498 retval.data(i) = - retval.data(i); |
|
|
499 return retval; |
|
|
500 } |
|
|
501 |
|
|
502 /* |
|
|
503 ;;; Local Variables: *** |
|
|
504 ;;; mode: C++ *** |
|
|
505 ;;; End: *** |
|
|
506 */ |
