|
5164
|
1 /* |
|
|
2 |
|
7017
|
3 Copyright (C) 2004, 2005, 2006, 2007 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" |
|
|
30 #include "MArray2.h" |
|
|
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 { |
|
|
71 OCTAVE_QUIT; |
|
|
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 { |
|
|
84 r.ridx(jx) = b.ridx(jb); |
|
|
85 r.data(jx) = 0. + b.data(jb); |
|
|
86 jx++; |
|
|
87 jb++; |
|
|
88 jb_lt_max= jb < jb_max; |
|
|
89 } |
|
|
90 else |
|
|
91 { |
|
|
92 if ((a.data(ja) + b.data(jb)) != 0.) |
|
|
93 { |
|
|
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 |
|
|
107 a = r.maybe_compress (); |
|
|
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 { |
|
|
144 OCTAVE_QUIT; |
|
|
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 { |
|
|
157 r.ridx(jx) = b.ridx(jb); |
|
|
158 r.data(jx) = 0. - b.data(jb); |
|
|
159 jx++; |
|
|
160 jb++; |
|
|
161 jb_lt_max= jb < jb_max; |
|
|
162 } |
|
|
163 else |
|
|
164 { |
|
|
165 if ((a.data(ja) - b.data(jb)) != 0.) |
|
|
166 { |
|
|
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 |
|
|
180 a = r.maybe_compress (); |
|
|
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> \ |
|
|
190 MArray2<T> \ |
|
|
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 \ |
|
|
196 MArray2<T> r (nr, nc, (0.0 OP s)); \ |
|
|
197 \ |
|
5275
|
198 for (octave_idx_type j = 0; j < nc; j++) \ |
|
|
199 for (octave_idx_type i = a.cidx(j); i < a.cidx(j+1); i++) \ |
|
5164
|
200 r.elem (a.ridx (i), j) = a.data (i) OP s; \ |
|
|
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 { \ |
|
|
217 r.data(i) = a.data(i) OP s; \ |
|
|
218 r.ridx(i) = a.ridx(i); \ |
|
|
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> \ |
|
|
236 MArray2<T> \ |
|
|
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 \ |
|
|
242 MArray2<T> r (nr, nc, (s OP 0.0)); \ |
|
|
243 \ |
|
5275
|
244 for (octave_idx_type j = 0; j < nc; j++) \ |
|
|
245 for (octave_idx_type i = a.cidx(j); i < a.cidx(j+1); i++) \ |
|
5164
|
246 r.elem (a.ridx (i), j) = s OP a.data (i); \ |
|
|
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 { \ |
|
|
263 r.data(i) = s OP a.data(i); \ |
|
|
264 r.ridx(i) = a.ridx(i); \ |
|
|
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 { \ |
|
|
298 r = MSparse<T> (b_nr, b_nc, a.data(0) OP 0.); \ |
|
|
299 \ |
|
|
300 for (octave_idx_type j = 0 ; j < b_nc ; j++) \ |
|
|
301 { \ |
|
|
302 OCTAVE_QUIT; \ |
|
|
303 octave_idx_type idxj = j * b_nr; \ |
|
|
304 for (octave_idx_type i = b.cidx(j) ; i < b.cidx(j+1) ; i++) \ |
|
|
305 { \ |
|
|
306 OCTAVE_QUIT; \ |
|
|
307 r.data(idxj + b.ridx(i)) = a.data(0) OP b.data(i); \ |
|
|
308 } \ |
|
|
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 { \ |
|
|
319 r = MSparse<T> (a_nr, a_nc, 0. OP b.data(0)); \ |
|
|
320 \ |
|
|
321 for (octave_idx_type j = 0 ; j < a_nc ; j++) \ |
|
|
322 { \ |
|
|
323 OCTAVE_QUIT; \ |
|
|
324 octave_idx_type idxj = j * a_nr; \ |
|
|
325 for (octave_idx_type i = a.cidx(j) ; i < a.cidx(j+1) ; i++) \ |
|
|
326 { \ |
|
|
327 OCTAVE_QUIT; \ |
|
|
328 r.data(idxj + a.ridx(i)) = a.data(i) OP b.data(0); \ |
|
|
329 } \ |
|
|
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; \ |
|
5164
|
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 { \ |
|
|
354 OCTAVE_QUIT; \ |
|
|
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 { \ |
|
|
367 r.ridx(jx) = b.ridx(jb); \ |
|
|
368 r.data(jx) = 0. OP b.data(jb); \ |
|
|
369 jx++; \ |
|
|
370 jb++; \ |
|
|
371 jb_lt_max= jb < jb_max; \ |
|
|
372 } \ |
|
|
373 else \ |
|
|
374 { \ |
|
|
375 if ((a.data(ja) OP b.data(jb)) != 0.) \ |
|
|
376 { \ |
|
|
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 \ |
|
|
390 r.maybe_compress (); \ |
|
|
391 } \ |
|
|
392 \ |
|
|
393 return r; \ |
|
|
394 } |
|
|
395 |
|
|
396 #define SPARSE_A2A2_FCN_1(FCN, OP) \ |
|
|
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 { \ |
|
|
415 r = MSparse<T> (b); \ |
|
|
416 octave_idx_type b_nnz = b.nnz(); \ |
|
|
417 \ |
|
|
418 for (octave_idx_type i = 0 ; i < b_nnz ; i++) \ |
|
|
419 { \ |
|
|
420 OCTAVE_QUIT; \ |
|
|
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 { \ |
|
|
432 r = MSparse<T> (a); \ |
|
|
433 octave_idx_type a_nnz = a.nnz(); \ |
|
|
434 \ |
|
|
435 for (octave_idx_type i = 0 ; i < a_nnz ; i++) \ |
|
|
436 { \ |
|
|
437 OCTAVE_QUIT; \ |
|
|
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; \ |
|
5164
|
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 { \ |
|
|
463 OCTAVE_QUIT; \ |
|
|
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 { \ |
|
|
476 if ((a.data(ja) OP b.data(jb)) != 0.) \ |
|
|
477 { \ |
|
|
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 \ |
|
|
489 r.maybe_compress (); \ |
|
|
490 } \ |
|
|
491 \ |
|
|
492 return r; \ |
|
|
493 } |
|
|
494 |
|
|
495 #define SPARSE_A2A2_FCN_2(FCN, OP) \ |
|
|
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 { \ |
|
|
526 OCTAVE_QUIT; \ |
|
|
527 octave_idx_type idxj = j * b_nr; \ |
|
|
528 for (octave_idx_type i = b.cidx(j) ; i < b.cidx(j+1) ; i++) \ |
|
|
529 { \ |
|
|
530 OCTAVE_QUIT; \ |
|
|
531 r.data(idxj + b.ridx(i)) = val OP b.data(i); \ |
|
|
532 } \ |
|
|
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 { \ |
|
|
554 OCTAVE_QUIT; \ |
|
|
555 octave_idx_type idxj = j * a_nr; \ |
|
|
556 for (octave_idx_type i = a.cidx(j) ; i < a.cidx(j+1) ; i++) \ |
|
|
557 { \ |
|
|
558 OCTAVE_QUIT; \ |
|
|
559 r.data(idxj + a.ridx(i)) = a.data(i) OP val; \ |
|
|
560 } \ |
|
|
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 { \ |
|
|
583 OCTAVE_QUIT; \ |
|
|
584 if ((! jb_lt_max) || \ |
|
|
585 (ja_lt_max && (a.ridx(ja) < b.ridx(jb)))) \ |
|
|
586 { \ |
|
|
587 r.elem (a.ridx(ja),i) = a.data(ja) OP Zero; \ |
|
|
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 { \ |
|
|
593 r.elem (b.ridx(jb),i) = Zero OP b.data(jb); \ |
|
|
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 \ |
|
|
605 r.maybe_compress (true); \ |
|
|
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 } |
|
|
635 |
|
|
636 /* |
|
|
637 ;;; Local Variables: *** |
|
|
638 ;;; mode: C++ *** |
|
|
639 ;;; End: *** |
|
|
640 */ |
