Mercurial > jwe > octave

comparison liboctave/array/MSparse.cc @ 21120:499b851fbfae

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Replace pattern if/err_XXX/else/code with if/err_XXX/ code. * schur.cc, ov-complex.h, ov-cx-mat.cc, ov-cx-sparse.cc, ov-float.h, ov-flt-complex.h, ov-flt-cx-mat.cc, ov-flt-re-mat.cc, ov-range.cc, ov-re-mat.cc, ov-re-sparse.cc, ov-scalar.h, ov-str-mat.cc, ops.h, pt-idx.cc, Array.cc, CColVector.cc, CMatrix.cc, CRowVector.cc, MSparse.cc, PermMatrix.cc, Sparse.cc, dColVector.cc, dMatrix.cc, dRowVector.cc, dSparse.cc, fCColVector.cc, fCMatrix.cc, fCRowVector.cc, fColVector.cc, fMatrix.cc, fRowVector.cc: Replace pattern if/err_XXX/else/code with if/err_XXX/ code.
author Rik <rik@octave.org>
date Wed, 20 Jan 2016 16:58:59 -0800
parents e39e05d90788
children 538b57866b90
comparison
equal deleted inserted replaced
21119:90cd0f9442d5 21120:499b851fbfae
37 octave_idx_type b_nr = b.rows (); 37 octave_idx_type b_nr = b.rows ();
38 octave_idx_type b_nc = b.cols (); 38 octave_idx_type b_nc = b.cols ();
39 39
40 if (a_nr != b_nr || a_nc != b_nc) 40 if (a_nr != b_nr || a_nc != b_nc)
41 err_nonconformant (op_name , a_nr, a_nc, b_nr, b_nc); 41 err_nonconformant (op_name , a_nr, a_nc, b_nr, b_nc);
42
43 r = MSparse<T> (a_nr, a_nc, (a.nnz () + b.nnz ()));
44
45 octave_idx_type jx = 0;
46 for (octave_idx_type i = 0 ; i < a_nc ; i++)
47 {
48 octave_idx_type ja = a.cidx (i);
49 octave_idx_type ja_max = a.cidx (i+1);
50 bool ja_lt_max= ja < ja_max;
51
52 octave_idx_type jb = b.cidx (i);
53 octave_idx_type jb_max = b.cidx (i+1);
54 bool jb_lt_max = jb < jb_max;
55
56 while (ja_lt_max || jb_lt_max)
57 {
58 octave_quit ();
59 if ((! jb_lt_max) || (ja_lt_max && (a.ridx (ja) < b.ridx (jb))))
60 {
61 r.ridx (jx) = a.ridx (ja);
62 r.data (jx) = op (a.data (ja), 0.);
63 jx++;
64 ja++;
65 ja_lt_max= ja < ja_max;
66 }
67 else if ((! ja_lt_max)
68 || (jb_lt_max && (b.ridx (jb) < a.ridx (ja))))
69 {
70 r.ridx (jx) = b.ridx (jb);
71 r.data (jx) = op (0., b.data (jb));
72 jx++;
73 jb++;
74 jb_lt_max= jb < jb_max;
75 }
76 else
77 {
78 if (op (a.data (ja), b.data (jb)) != 0.)
79 {
80 r.data (jx) = op (a.data (ja), b.data (jb));
81 r.ridx (jx) = a.ridx (ja);
82 jx++;
83 }
84 ja++;
85 ja_lt_max= ja < ja_max;
86 jb++;
87 jb_lt_max= jb < jb_max;
88 }
89 }
90 r.cidx (i+1) = jx;
91 }
92
93 a = r.maybe_compress ();
94
95 return a;
96 }
97
98 template <typename T>
99 MSparse<T>&
100 operator += (MSparse<T>& a, const MSparse<T>& b)
101 {
102 return plus_or_minus (a, b, std::plus<T> (), "operator +=");
103 }
104
105 template <typename T>
106 MSparse<T>&
107 operator -= (MSparse<T>& a, const MSparse<T>& b)
108 {
109 return plus_or_minus (a, b, std::minus<T> (), "operator -=");
110 }
111
112
113 // Element by element MSparse by scalar ops.
114
115 template <class T, class OP>
116 MArray<T>
117 plus_or_minus (const MSparse<T>& a, const T& s, OP op)
118 {
119 octave_idx_type nr = a.rows ();
120 octave_idx_type nc = a.cols ();
121
122 MArray<T> r (dim_vector (nr, nc), op (0.0, s));
123
124 for (octave_idx_type j = 0; j < nc; j++)
125 for (octave_idx_type i = a.cidx (j); i < a.cidx (j+1); i++)
126 r.elem (a.ridx (i), j) = op (a.data (i), s);
127 return r;
128 }
129
130 template <typename T>
131 MArray<T>
132 operator + (const MSparse<T>& a, const T& s)
133 {
134 return plus_or_minus (a, s, std::plus<T> ());
135 }
136
137 template <typename T>
138 MArray<T>
139 operator - (const MSparse<T>& a, const T& s)
140 {
141 return plus_or_minus (a, s, std::minus<T> ());
142 }
143
144
145 template <class T, class OP>
146 MSparse<T>
147 times_or_divide (const MSparse<T>& a, const T& s, OP op)
148 {
149 octave_idx_type nr = a.rows ();
150 octave_idx_type nc = a.cols ();
151 octave_idx_type nz = a.nnz ();
152
153 MSparse<T> r (nr, nc, nz);
154
155 for (octave_idx_type i = 0; i < nz; i++)
156 {
157 r.data (i) = op (a.data (i), s);
158 r.ridx (i) = a.ridx (i);
159 }
160 for (octave_idx_type i = 0; i < nc + 1; i++)
161 r.cidx (i) = a.cidx (i);
162 r.maybe_compress (true);
163 return r;
164 }
165
166 template <typename T>
167 MSparse<T>
168 operator * (const MSparse<T>& a, const T& s)
169 {
170 return times_or_divide (a, s, std::multiplies<T> ());
171 }
172
173 template <typename T>
174 MSparse<T>
175 operator / (const MSparse<T>& a, const T& s)
176 {
177 return times_or_divide (a, s, std::divides<T> ());
178 }
179
180
181 // Element by element scalar by MSparse ops.
182
183 template <class T, class OP>
184 MArray<T>
185 plus_or_minus (const T& s, const MSparse<T>& a, OP op)
186 {
187 octave_idx_type nr = a.rows ();
188 octave_idx_type nc = a.cols ();
189
190 MArray<T> r (dim_vector (nr, nc), op (s, 0.0));
191
192 for (octave_idx_type j = 0; j < nc; j++)
193 for (octave_idx_type i = a.cidx (j); i < a.cidx (j+1); i++)
194 r.elem (a.ridx (i), j) = op (s, a.data (i));
195 return r;
196 }
197
198 template <typename T>
199 MArray<T>
200 operator + (const T& s, const MSparse<T>& a)
201 {
202 return plus_or_minus (s, a, std::plus<T> ());
203 }
204
205 template <typename T>
206 MArray<T>
207 operator - (const T& s, const MSparse<T>& a)
208 {
209 return plus_or_minus (s, a, std::minus<T> ());
210 }
211
212 template <class T, class OP>
213 MSparse<T>
214 times_or_divides (const T& s, const MSparse<T>& a, OP op)
215 {
216 octave_idx_type nr = a.rows ();
217 octave_idx_type nc = a.cols ();
218 octave_idx_type nz = a.nnz ();
219
220 MSparse<T> r (nr, nc, nz);
221
222 for (octave_idx_type i = 0; i < nz; i++)
223 {
224 r.data (i) = op (s, a.data (i));
225 r.ridx (i) = a.ridx (i);
226 }
227 for (octave_idx_type i = 0; i < nc + 1; i++)
228 r.cidx (i) = a.cidx (i);
229 r.maybe_compress (true);
230 return r;
231 }
232
233 template <class T>
234 MSparse<T>
235 operator * (const T& s, const MSparse<T>& a)
236 {
237 return times_or_divides (s, a, std::multiplies<T> ());
238 }
239
240 template <class T>
241 MSparse<T>
242 operator / (const T& s, const MSparse<T>& a)
243 {
244 return times_or_divides (s, a, std::divides<T> ());
245 }
246
247
248 // Element by element MSparse by MSparse ops.
249
250 template <class T, class OP>
251 MSparse<T>
252 plus_or_minus (const MSparse<T>& a, const MSparse<T>& b, OP op,
253 const char* op_name, bool negate)
254 {
255 MSparse<T> r;
256
257 octave_idx_type a_nr = a.rows ();
258 octave_idx_type a_nc = a.cols ();
259
260 octave_idx_type b_nr = b.rows ();
261 octave_idx_type b_nc = b.cols ();
262
263 if (a_nr == 1 && a_nc == 1)
264 {
265 if (a.elem (0,0) == 0.)
266 if (negate)
267 r = -MSparse<T> (b);
268 else
269 r = MSparse<T> (b);
270 else
271 {
272 r = MSparse<T> (b_nr, b_nc, op (a.data (0), 0.));
273
274 for (octave_idx_type j = 0 ; j < b_nc ; j++)
275 {
276 octave_quit ();
277 octave_idx_type idxj = j * b_nr;
278 for (octave_idx_type i = b.cidx (j) ; i < b.cidx (j+1) ; i++)
279 {
280 octave_quit ();
281 r.data (idxj + b.ridx (i)) = op (a.data (0), b.data (i));
282 }
283 }
284 r.maybe_compress ();
285 }
286 }
287 else if (b_nr == 1 && b_nc == 1)
288 {
289 if (b.elem (0,0) == 0.)
290 r = MSparse<T> (a);
291 else
292 {
293 r = MSparse<T> (a_nr, a_nc, op (0.0, b.data (0)));
294
295 for (octave_idx_type j = 0 ; j < a_nc ; j++)
296 {
297 octave_quit ();
298 octave_idx_type idxj = j * a_nr;
299 for (octave_idx_type i = a.cidx (j) ; i < a.cidx (j+1) ; i++)
300 {
301 octave_quit ();
302 r.data (idxj + a.ridx (i)) = op (a.data (i), b.data (0));
303 }
304 }
305 r.maybe_compress ();
306 }
307 }
308 else if (a_nr != b_nr || a_nc != b_nc)
309 err_nonconformant (op_name, a_nr, a_nc, b_nr, b_nc);
42 else 310 else
43 { 311 {
44 r = MSparse<T> (a_nr, a_nc, (a.nnz () + b.nnz ())); 312 r = MSparse<T> (a_nr, a_nc, (a.nnz () + b.nnz ()));
45 313
46 octave_idx_type jx = 0; 314 octave_idx_type jx = 0;
315 r.cidx (0) = 0;
47 for (octave_idx_type i = 0 ; i < a_nc ; i++) 316 for (octave_idx_type i = 0 ; i < a_nc ; i++)
48 { 317 {
49 octave_idx_type ja = a.cidx (i); 318 octave_idx_type ja = a.cidx (i);
50 octave_idx_type ja_max = a.cidx (i+1); 319 octave_idx_type ja_max = a.cidx (i+1);
51 bool ja_lt_max= ja < ja_max; 320 bool ja_lt_max= ja < ja_max;
67 } 336 }
68 else if ((! ja_lt_max) 337 else if ((! ja_lt_max)
69 || (jb_lt_max && (b.ridx (jb) < a.ridx (ja)))) 338 || (jb_lt_max && (b.ridx (jb) < a.ridx (ja))))
70 { 339 {
71 r.ridx (jx) = b.ridx (jb); 340 r.ridx (jx) = b.ridx (jb);
72 r.data (jx) = op (0., b.data (jb)); 341 r.data (jx) = op (0., b.data (jb));
73 jx++; 342 jx++;
74 jb++; 343 jb++;
75 jb_lt_max= jb < jb_max; 344 jb_lt_max= jb < jb_max;
76 } 345 }
77 else 346 else
89 } 358 }
90 } 359 }
91 r.cidx (i+1) = jx; 360 r.cidx (i+1) = jx;
92 } 361 }
93 362
94 a = r.maybe_compress ();
95 }
96
97 return a;
98 }
99
100 template <typename T>
101 MSparse<T>&
102 operator += (MSparse<T>& a, const MSparse<T>& b)
103 {
104 return plus_or_minus (a, b, std::plus<T> (), "operator +=");
105 }
106
107 template <typename T>
108 MSparse<T>&
109 operator -= (MSparse<T>& a, const MSparse<T>& b)
110 {
111 return plus_or_minus (a, b, std::minus<T> (), "operator -=");
112 }
113
114
115 // Element by element MSparse by scalar ops.
116
117 template <class T, class OP>
118 MArray<T>
119 plus_or_minus (const MSparse<T>& a, const T& s, OP op)
120 {
121 octave_idx_type nr = a.rows ();
122 octave_idx_type nc = a.cols ();
123
124 MArray<T> r (dim_vector (nr, nc), op (0.0, s));
125
126 for (octave_idx_type j = 0; j < nc; j++)
127 for (octave_idx_type i = a.cidx (j); i < a.cidx (j+1); i++)
128 r.elem (a.ridx (i), j) = op (a.data (i), s);
129 return r;
130 }
131
132 template <typename T>
133 MArray<T>
134 operator + (const MSparse<T>& a, const T& s)
135 {
136 return plus_or_minus (a, s, std::plus<T> ());
137 }
138
139 template <typename T>
140 MArray<T>
141 operator - (const MSparse<T>& a, const T& s)
142 {
143 return plus_or_minus (a, s, std::minus<T> ());
144 }
145
146
147 template <class T, class OP>
148 MSparse<T>
149 times_or_divide (const MSparse<T>& a, const T& s, OP op)
150 {
151 octave_idx_type nr = a.rows ();
152 octave_idx_type nc = a.cols ();
153 octave_idx_type nz = a.nnz ();
154
155 MSparse<T> r (nr, nc, nz);
156
157 for (octave_idx_type i = 0; i < nz; i++)
158 {
159 r.data (i) = op (a.data (i), s);
160 r.ridx (i) = a.ridx (i);
161 }
162 for (octave_idx_type i = 0; i < nc + 1; i++)
163 r.cidx (i) = a.cidx (i);
164 r.maybe_compress (true);
165 return r;
166 }
167
168 template <typename T>
169 MSparse<T>
170 operator * (const MSparse<T>& a, const T& s)
171 {
172 return times_or_divide (a, s, std::multiplies<T> ());
173 }
174
175 template <typename T>
176 MSparse<T>
177 operator / (const MSparse<T>& a, const T& s)
178 {
179 return times_or_divide (a, s, std::divides<T> ());
180 }
181
182
183 // Element by element scalar by MSparse ops.
184
185 template <class T, class OP>
186 MArray<T>
187 plus_or_minus (const T& s, const MSparse<T>& a, OP op)
188 {
189 octave_idx_type nr = a.rows ();
190 octave_idx_type nc = a.cols ();
191
192 MArray<T> r (dim_vector (nr, nc), op (s, 0.0));
193
194 for (octave_idx_type j = 0; j < nc; j++)
195 for (octave_idx_type i = a.cidx (j); i < a.cidx (j+1); i++)
196 r.elem (a.ridx (i), j) = op (s, a.data (i));
197 return r;
198 }
199
200 template <typename T>
201 MArray<T>
202 operator + (const T& s, const MSparse<T>& a)
203 {
204 return plus_or_minus (s, a, std::plus<T> ());
205 }
206
207 template <typename T>
208 MArray<T>
209 operator - (const T& s, const MSparse<T>& a)
210 {
211 return plus_or_minus (s, a, std::minus<T> ());
212 }
213
214 template <class T, class OP>
215 MSparse<T>
216 times_or_divides (const T& s, const MSparse<T>& a, OP op)
217 {
218 octave_idx_type nr = a.rows ();
219 octave_idx_type nc = a.cols ();
220 octave_idx_type nz = a.nnz ();
221
222 MSparse<T> r (nr, nc, nz);
223
224 for (octave_idx_type i = 0; i < nz; i++)
225 {
226 r.data (i) = op (s, a.data (i));
227 r.ridx (i) = a.ridx (i);
228 }
229 for (octave_idx_type i = 0; i < nc + 1; i++)
230 r.cidx (i) = a.cidx (i);
231 r.maybe_compress (true);
232 return r;
233 }
234
235 template <class T>
236 MSparse<T>
237 operator * (const T& s, const MSparse<T>& a)
238 {
239 return times_or_divides (s, a, std::multiplies<T> ());
240 }
241
242 template <class T>
243 MSparse<T>
244 operator / (const T& s, const MSparse<T>& a)
245 {
246 return times_or_divides (s, a, std::divides<T> ());
247 }
248
249
250 // Element by element MSparse by MSparse ops.
251
252 template <class T, class OP>
253 MSparse<T>
254 plus_or_minus (const MSparse<T>& a, const MSparse<T>& b, OP op,
255 const char* op_name, bool negate)
256 {
257 MSparse<T> r;
258
259 octave_idx_type a_nr = a.rows ();
260 octave_idx_type a_nc = a.cols ();
261
262 octave_idx_type b_nr = b.rows ();
263 octave_idx_type b_nc = b.cols ();
264
265 if (a_nr == 1 && a_nc == 1)
266 {
267 if (a.elem (0,0) == 0.)
268 if (negate)
269 r = -MSparse<T> (b);
270 else
271 r = MSparse<T> (b);
272 else
273 {
274 r = MSparse<T> (b_nr, b_nc, op (a.data (0), 0.));
275
276 for (octave_idx_type j = 0 ; j < b_nc ; j++)
277 {
278 octave_quit ();
279 octave_idx_type idxj = j * b_nr;
280 for (octave_idx_type i = b.cidx (j) ; i < b.cidx (j+1) ; i++)
281 {
282 octave_quit ();
283 r.data (idxj + b.ridx (i)) = op (a.data (0), b.data (i));
284 }
285 }
286 r.maybe_compress ();
287 }
288 }
289 else if (b_nr == 1 && b_nc == 1)
290 {
291 if (b.elem (0,0) == 0.)
292 r = MSparse<T> (a);
293 else
294 {
295 r = MSparse<T> (a_nr, a_nc, op (0.0, b.data (0)));
296
297 for (octave_idx_type j = 0 ; j < a_nc ; j++)
298 {
299 octave_quit ();
300 octave_idx_type idxj = j * a_nr;
301 for (octave_idx_type i = a.cidx (j) ; i < a.cidx (j+1) ; i++)
302 {
303 octave_quit ();
304 r.data (idxj + a.ridx (i)) = op (a.data (i), b.data (0));
305 }
306 }
307 r.maybe_compress ();
308 }
309 }
310 else if (a_nr != b_nr || a_nc != b_nc)
311 err_nonconformant (op_name, a_nr, a_nc, b_nr, b_nc);
312 else
313 {
314 r = MSparse<T> (a_nr, a_nc, (a.nnz () + b.nnz ()));
315
316 octave_idx_type jx = 0;
317 r.cidx (0) = 0;
318 for (octave_idx_type i = 0 ; i < a_nc ; i++)
319 {
320 octave_idx_type ja = a.cidx (i);
321 octave_idx_type ja_max = a.cidx (i+1);
322 bool ja_lt_max= ja < ja_max;
323
324 octave_idx_type jb = b.cidx (i);
325 octave_idx_type jb_max = b.cidx (i+1);
326 bool jb_lt_max = jb < jb_max;
327
328 while (ja_lt_max || jb_lt_max)
329 {
330 octave_quit ();
331 if ((! jb_lt_max) || (ja_lt_max && (a.ridx (ja) < b.ridx (jb))))
332 {
333 r.ridx (jx) = a.ridx (ja);
334 r.data (jx) = op (a.data (ja), 0.);
335 jx++;
336 ja++;
337 ja_lt_max= ja < ja_max;
338 }
339 else if ((! ja_lt_max)
340 || (jb_lt_max && (b.ridx (jb) < a.ridx (ja))))
341 {
342 r.ridx (jx) = b.ridx (jb);
343 r.data (jx) = op (0., b.data (jb));
344 jx++;
345 jb++;
346 jb_lt_max= jb < jb_max;
347 }
348 else
349 {
350 if (op (a.data (ja), b.data (jb)) != 0.)
351 {
352 r.data (jx) = op (a.data (ja), b.data (jb));
353 r.ridx (jx) = a.ridx (ja);
354 jx++;
355 }
356 ja++;
357 ja_lt_max= ja < ja_max;
358 jb++;
359 jb_lt_max= jb < jb_max;
360 }
361 }
362 r.cidx (i+1) = jx;
363 }
364
365 r.maybe_compress (); 363 r.maybe_compress ();
366 } 364 }
367 365
368 return r; 366 return r;
369 } 367 }