Mercurial > octave-nkf
annotate liboctave/Sparse.cc @ 8710:739141cde75a ss-3-1-52
fix typo in Array-f.cc
author | Jaroslav Hajek <highegg@gmail.com> |
---|---|
date | Mon, 09 Feb 2009 21:51:31 +0100 |
parents | 095ae5e0a831 |
children | 06b9903a029b |
rev | line source |
---|---|
5164 | 1 // Template sparse array class |
2 /* | |
3 | |
7017 | 4 Copyright (C) 2004, 2005, 2006, 2007 David Bateman |
7016 | 5 Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004 Andy Adler |
6 | |
7 This file is part of Octave. | |
5164 | 8 |
9 Octave is free software; you can redistribute it and/or modify it | |
10 under the terms of the GNU General Public License as published by the | |
7016 | 11 Free Software Foundation; either version 3 of the License, or (at your |
12 option) any later version. | |
5164 | 13 |
14 Octave is distributed in the hope that it will be useful, but WITHOUT | |
15 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or | |
16 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License | |
17 for more details. | |
18 | |
19 You should have received a copy of the GNU General Public License | |
7016 | 20 along with Octave; see the file COPYING. If not, see |
21 <http://www.gnu.org/licenses/>. | |
5164 | 22 |
23 */ | |
24 | |
25 #ifdef HAVE_CONFIG_H | |
26 #include <config.h> | |
27 #endif | |
28 | |
29 #include <cassert> | |
30 #include <climits> | |
31 | |
32 #include <iostream> | |
5765 | 33 #include <sstream> |
5164 | 34 #include <vector> |
35 | |
36 #include "Array.h" | |
37 #include "Array-util.h" | |
38 #include "Range.h" | |
39 #include "idx-vector.h" | |
40 #include "lo-error.h" | |
41 #include "quit.h" | |
8377
25bc2d31e1bf
improve OCTAVE_LOCAL_BUFFER
Jaroslav Hajek <highegg@gmail.com>
parents:
8290
diff
changeset
|
42 #include "oct-locbuf.h" |
5164 | 43 |
44 #include "Sparse.h" | |
45 #include "sparse-sort.h" | |
46 #include "oct-spparms.h" | |
47 | |
48 template <class T> | |
49 T& | |
5275 | 50 Sparse<T>::SparseRep::elem (octave_idx_type _r, octave_idx_type _c) |
5164 | 51 { |
5275 | 52 octave_idx_type i; |
5164 | 53 |
5604 | 54 if (nzmx > 0) |
5164 | 55 { |
56 for (i = c[_c]; i < c[_c + 1]; i++) | |
57 if (r[i] == _r) | |
58 return d[i]; | |
59 else if (r[i] > _r) | |
60 break; | |
61 | |
62 // Ok, If we've gotten here, we're in trouble.. Have to create a | |
63 // new element in the sparse array. This' gonna be slow!!! | |
5869 | 64 if (c[ncols] == nzmx) |
5164 | 65 { |
66 (*current_liboctave_error_handler) | |
5275 | 67 ("Sparse::SparseRep::elem (octave_idx_type, octave_idx_type): sparse matrix filled"); |
5164 | 68 return *d; |
69 } | |
70 | |
5275 | 71 octave_idx_type to_move = c[ncols] - i; |
5164 | 72 if (to_move != 0) |
73 { | |
5275 | 74 for (octave_idx_type j = c[ncols]; j > i; j--) |
5164 | 75 { |
76 d[j] = d[j-1]; | |
77 r[j] = r[j-1]; | |
78 } | |
79 } | |
80 | |
5275 | 81 for (octave_idx_type j = _c + 1; j < ncols + 1; j++) |
5164 | 82 c[j] = c[j] + 1; |
83 | |
84 d[i] = 0.; | |
85 r[i] = _r; | |
86 | |
87 return d[i]; | |
88 } | |
89 else | |
90 { | |
91 (*current_liboctave_error_handler) | |
5275 | 92 ("Sparse::SparseRep::elem (octave_idx_type, octave_idx_type): sparse matrix filled"); |
5164 | 93 return *d; |
94 } | |
95 } | |
96 | |
97 template <class T> | |
98 T | |
5275 | 99 Sparse<T>::SparseRep::celem (octave_idx_type _r, octave_idx_type _c) const |
5164 | 100 { |
5604 | 101 if (nzmx > 0) |
5275 | 102 for (octave_idx_type i = c[_c]; i < c[_c + 1]; i++) |
5164 | 103 if (r[i] == _r) |
104 return d[i]; | |
105 return T (); | |
106 } | |
107 | |
108 template <class T> | |
109 void | |
110 Sparse<T>::SparseRep::maybe_compress (bool remove_zeros) | |
111 { | |
5604 | 112 octave_idx_type ndel = nzmx - c[ncols]; |
5275 | 113 octave_idx_type nzero = 0; |
5164 | 114 |
115 if (remove_zeros) | |
5604 | 116 for (octave_idx_type i = 0; i < nzmx - ndel; i++) |
5164 | 117 if (d[i] == T ()) |
118 nzero++; | |
119 | |
120 if (!ndel && !nzero) | |
121 return; | |
122 | |
123 if (!nzero) | |
124 { | |
5604 | 125 octave_idx_type new_nzmx = nzmx - ndel; |
126 | |
127 T *new_data = new T [new_nzmx]; | |
128 for (octave_idx_type i = 0; i < new_nzmx; i++) | |
5164 | 129 new_data[i] = d[i]; |
130 delete [] d; | |
131 d = new_data; | |
132 | |
5604 | 133 octave_idx_type *new_ridx = new octave_idx_type [new_nzmx]; |
134 for (octave_idx_type i = 0; i < new_nzmx; i++) | |
5164 | 135 new_ridx[i] = r[i]; |
136 delete [] r; | |
137 r = new_ridx; | |
138 } | |
139 else | |
140 { | |
5604 | 141 octave_idx_type new_nzmx = nzmx - ndel - nzero; |
142 | |
143 T *new_data = new T [new_nzmx]; | |
144 octave_idx_type *new_ridx = new octave_idx_type [new_nzmx]; | |
5275 | 145 |
146 octave_idx_type ii = 0; | |
147 octave_idx_type ic = 0; | |
148 for (octave_idx_type j = 0; j < ncols; j++) | |
5164 | 149 { |
5275 | 150 for (octave_idx_type k = ic; k < c[j+1]; k++) |
5164 | 151 if (d[k] != T ()) |
152 { | |
153 new_data [ii] = d[k]; | |
154 new_ridx [ii++] = r[k]; | |
155 } | |
156 ic = c[j+1]; | |
157 c[j+1] = ii; | |
158 } | |
159 | |
160 delete [] d; | |
161 d = new_data; | |
162 | |
163 delete [] r; | |
164 r = new_ridx; | |
165 } | |
166 | |
5604 | 167 nzmx -= ndel + nzero; |
5164 | 168 } |
169 | |
170 template <class T> | |
171 void | |
5275 | 172 Sparse<T>::SparseRep::change_length (octave_idx_type nz) |
5164 | 173 { |
5604 | 174 if (nz != nzmx) |
5164 | 175 { |
5604 | 176 octave_idx_type min_nzmx = (nz < nzmx ? nz : nzmx); |
5275 | 177 |
178 octave_idx_type * new_ridx = new octave_idx_type [nz]; | |
5604 | 179 for (octave_idx_type i = 0; i < min_nzmx; i++) |
5164 | 180 new_ridx[i] = r[i]; |
181 | |
182 delete [] r; | |
183 r = new_ridx; | |
184 | |
185 T * new_data = new T [nz]; | |
5604 | 186 for (octave_idx_type i = 0; i < min_nzmx; i++) |
5164 | 187 new_data[i] = d[i]; |
188 | |
189 delete [] d; | |
190 d = new_data; | |
191 | |
5604 | 192 if (nz < nzmx) |
5275 | 193 for (octave_idx_type i = 0; i <= ncols; i++) |
5164 | 194 if (c[i] > nz) |
195 c[i] = nz; | |
196 | |
5604 | 197 nzmx = nz; |
5164 | 198 } |
199 } | |
200 | |
201 template <class T> | |
202 template <class U> | |
203 Sparse<T>::Sparse (const Sparse<U>& a) | |
204 : dimensions (a.dimensions), idx (0), idx_count (0) | |
205 { | |
5681 | 206 if (a.nnz () == 0) |
5164 | 207 rep = new typename Sparse<T>::SparseRep (rows (), cols()); |
208 else | |
209 { | |
5681 | 210 rep = new typename Sparse<T>::SparseRep (rows (), cols (), a.nnz ()); |
5164 | 211 |
5681 | 212 octave_idx_type nz = a.nnz (); |
5275 | 213 octave_idx_type nc = cols (); |
214 for (octave_idx_type i = 0; i < nz; i++) | |
5164 | 215 { |
216 xdata (i) = T (a.data (i)); | |
217 xridx (i) = a.ridx (i); | |
218 } | |
5275 | 219 for (octave_idx_type i = 0; i < nc + 1; i++) |
5164 | 220 xcidx (i) = a.cidx (i); |
221 } | |
222 } | |
223 | |
224 template <class T> | |
5275 | 225 Sparse<T>::Sparse (octave_idx_type nr, octave_idx_type nc, T val) |
5630 | 226 : dimensions (dim_vector (nr, nc)), idx (0), idx_count (0) |
5164 | 227 { |
5630 | 228 if (val != T ()) |
5164 | 229 { |
5630 | 230 rep = new typename Sparse<T>::SparseRep (nr, nc, nr*nc); |
231 | |
232 octave_idx_type ii = 0; | |
233 xcidx (0) = 0; | |
234 for (octave_idx_type j = 0; j < nc; j++) | |
5164 | 235 { |
5630 | 236 for (octave_idx_type i = 0; i < nr; i++) |
237 { | |
238 xdata (ii) = val; | |
239 xridx (ii++) = i; | |
240 } | |
241 xcidx (j+1) = ii; | |
242 } | |
243 } | |
244 else | |
245 { | |
246 rep = new typename Sparse<T>::SparseRep (nr, nc, 0); | |
247 for (octave_idx_type j = 0; j < nc+1; j++) | |
248 xcidx(j) = 0; | |
5164 | 249 } |
250 } | |
251 | |
252 template <class T> | |
253 Sparse<T>::Sparse (const dim_vector& dv) | |
254 : dimensions (dv), idx (0), idx_count (0) | |
255 { | |
256 if (dv.length() != 2) | |
257 (*current_liboctave_error_handler) | |
258 ("Sparse::Sparse (const dim_vector&): dimension mismatch"); | |
259 else | |
260 rep = new typename Sparse<T>::SparseRep (dv(0), dv(1)); | |
261 } | |
262 | |
263 template <class T> | |
264 Sparse<T>::Sparse (const Sparse<T>& a, const dim_vector& dv) | |
265 : dimensions (dv), idx (0), idx_count (0) | |
266 { | |
267 | |
268 // Work in unsigned long long to avoid overflow issues with numel | |
269 unsigned long long a_nel = static_cast<unsigned long long>(a.rows ()) * | |
270 static_cast<unsigned long long>(a.cols ()); | |
271 unsigned long long dv_nel = static_cast<unsigned long long>(dv (0)) * | |
272 static_cast<unsigned long long>(dv (1)); | |
273 | |
274 if (a_nel != dv_nel) | |
275 (*current_liboctave_error_handler) | |
276 ("Sparse::Sparse (const Sparse&, const dim_vector&): dimension mismatch"); | |
277 else | |
278 { | |
279 dim_vector old_dims = a.dims(); | |
5681 | 280 octave_idx_type new_nzmx = a.nnz (); |
5275 | 281 octave_idx_type new_nr = dv (0); |
282 octave_idx_type new_nc = dv (1); | |
283 octave_idx_type old_nr = old_dims (0); | |
284 octave_idx_type old_nc = old_dims (1); | |
5164 | 285 |
5604 | 286 rep = new typename Sparse<T>::SparseRep (new_nr, new_nc, new_nzmx); |
5164 | 287 |
5275 | 288 octave_idx_type kk = 0; |
5164 | 289 xcidx(0) = 0; |
5275 | 290 for (octave_idx_type i = 0; i < old_nc; i++) |
291 for (octave_idx_type j = a.cidx(i); j < a.cidx(i+1); j++) | |
5164 | 292 { |
5275 | 293 octave_idx_type tmp = i * old_nr + a.ridx(j); |
294 octave_idx_type ii = tmp % new_nr; | |
295 octave_idx_type jj = (tmp - ii) / new_nr; | |
296 for (octave_idx_type k = kk; k < jj; k++) | |
5164 | 297 xcidx(k+1) = j; |
298 kk = jj; | |
299 xdata(j) = a.data(j); | |
300 xridx(j) = ii; | |
301 } | |
5275 | 302 for (octave_idx_type k = kk; k < new_nc; k++) |
5604 | 303 xcidx(k+1) = new_nzmx; |
5164 | 304 } |
305 } | |
306 | |
307 template <class T> | |
5275 | 308 Sparse<T>::Sparse (const Array<T>& a, const Array<octave_idx_type>& r, |
309 const Array<octave_idx_type>& c, octave_idx_type nr, | |
310 octave_idx_type nc, bool sum_terms) | |
5164 | 311 : dimensions (dim_vector (nr, nc)), idx (0), idx_count (0) |
312 { | |
5275 | 313 octave_idx_type a_len = a.length (); |
314 octave_idx_type r_len = r.length (); | |
315 octave_idx_type c_len = c.length (); | |
5164 | 316 bool ri_scalar = (r_len == 1); |
317 bool ci_scalar = (c_len == 1); | |
318 bool cf_scalar = (a_len == 1); | |
319 | |
320 if ((a_len != r_len && !cf_scalar && !ri_scalar) || | |
321 (a_len != c_len && !cf_scalar && !ci_scalar) || | |
322 (r_len != c_len && !ri_scalar && !ci_scalar) || nr < 0 || nc < 0) | |
323 { | |
324 (*current_liboctave_error_handler) | |
5275 | 325 ("Sparse::Sparse (const Array<T>&, const Array<octave_idx_type>&, ...): dimension mismatch"); |
5164 | 326 rep = nil_rep (); |
327 dimensions = dim_vector (0, 0); | |
328 } | |
329 else | |
330 { | |
5604 | 331 octave_idx_type max_nzmx = (r_len > c_len ? r_len : c_len); |
332 | |
333 OCTAVE_LOCAL_BUFFER (octave_sparse_sort_idxl *, sidx, max_nzmx); | |
334 OCTAVE_LOCAL_BUFFER (octave_sparse_sort_idxl, sidxX, max_nzmx); | |
335 | |
336 for (octave_idx_type i = 0; i < max_nzmx; i++) | |
5164 | 337 sidx[i] = &sidxX[i]; |
338 | |
5604 | 339 octave_idx_type actual_nzmx = 0; |
5164 | 340 OCTAVE_QUIT; |
5604 | 341 for (octave_idx_type i = 0; i < max_nzmx; i++) |
5164 | 342 { |
5275 | 343 octave_idx_type rowidx = (ri_scalar ? r(0) : r(i)); |
344 octave_idx_type colidx = (ci_scalar ? c(0) : c(i)); | |
5164 | 345 if (rowidx < nr && rowidx >= 0 && |
346 colidx < nc && colidx >= 0 ) | |
347 { | |
348 if ( a (cf_scalar ? 0 : i ) != T ()) | |
349 { | |
5604 | 350 sidx[actual_nzmx]->r = rowidx; |
351 sidx[actual_nzmx]->c = colidx; | |
352 sidx[actual_nzmx]->idx = i; | |
353 actual_nzmx++; | |
5164 | 354 } |
355 } | |
356 else | |
357 { | |
358 (*current_liboctave_error_handler) | |
359 ("Sparse::Sparse : index (%d,%d) out of range", | |
360 rowidx + 1, colidx + 1); | |
361 rep = nil_rep (); | |
362 dimensions = dim_vector (0, 0); | |
363 return; | |
364 } | |
365 } | |
366 | |
5604 | 367 if (actual_nzmx == 0) |
5164 | 368 rep = new typename Sparse<T>::SparseRep (nr, nc); |
369 else | |
370 { | |
371 OCTAVE_QUIT; | |
372 octave_sort<octave_sparse_sort_idxl *> | |
7433 | 373 lsort (octave_sparse_sidxl_comp); |
374 | |
375 lsort.sort (sidx, actual_nzmx); | |
5164 | 376 OCTAVE_QUIT; |
377 | |
378 // Now count the unique non-zero values | |
5604 | 379 octave_idx_type real_nzmx = 1; |
380 for (octave_idx_type i = 1; i < actual_nzmx; i++) | |
5164 | 381 if (sidx[i-1]->r != sidx[i]->r || sidx[i-1]->c != sidx[i]->c) |
5604 | 382 real_nzmx++; |
383 | |
384 rep = new typename Sparse<T>::SparseRep (nr, nc, real_nzmx); | |
5164 | 385 |
5275 | 386 octave_idx_type cx = 0; |
387 octave_idx_type prev_rval = -1; | |
388 octave_idx_type prev_cval = -1; | |
389 octave_idx_type ii = -1; | |
5164 | 390 xcidx (0) = 0; |
5604 | 391 for (octave_idx_type i = 0; i < actual_nzmx; i++) |
5164 | 392 { |
393 OCTAVE_QUIT; | |
5275 | 394 octave_idx_type iidx = sidx[i]->idx; |
395 octave_idx_type rval = sidx[i]->r; | |
396 octave_idx_type cval = sidx[i]->c; | |
5164 | 397 |
398 if (prev_cval < cval || (prev_rval < rval && prev_cval == cval)) | |
399 { | |
5275 | 400 octave_idx_type ci = static_cast<octave_idx_type> (c (ci_scalar ? 0 : iidx)); |
5164 | 401 ii++; |
402 while (cx < ci) | |
403 xcidx (++cx) = ii; | |
404 xdata(ii) = a (cf_scalar ? 0 : iidx); | |
5275 | 405 xridx(ii) = static_cast<octave_idx_type> (r (ri_scalar ? 0 : iidx)); |
5164 | 406 } |
407 else | |
408 { | |
409 if (sum_terms) | |
410 xdata(ii) += a (cf_scalar ? 0 : iidx); | |
411 else | |
412 xdata(ii) = a (cf_scalar ? 0 : iidx); | |
413 } | |
414 prev_rval = rval; | |
415 prev_cval = cval; | |
416 } | |
417 | |
418 while (cx < nc) | |
419 xcidx (++cx) = ii + 1; | |
420 } | |
421 } | |
422 } | |
423 | |
424 template <class T> | |
425 Sparse<T>::Sparse (const Array<T>& a, const Array<double>& r, | |
5275 | 426 const Array<double>& c, octave_idx_type nr, |
427 octave_idx_type nc, bool sum_terms) | |
5164 | 428 : dimensions (dim_vector (nr, nc)), idx (0), idx_count (0) |
429 { | |
5275 | 430 octave_idx_type a_len = a.length (); |
431 octave_idx_type r_len = r.length (); | |
432 octave_idx_type c_len = c.length (); | |
5164 | 433 bool ri_scalar = (r_len == 1); |
434 bool ci_scalar = (c_len == 1); | |
435 bool cf_scalar = (a_len == 1); | |
436 | |
437 if ((a_len != r_len && !cf_scalar && !ri_scalar) || | |
438 (a_len != c_len && !cf_scalar && !ci_scalar) || | |
439 (r_len != c_len && !ri_scalar && !ci_scalar) || nr < 0 || nc < 0) | |
440 { | |
441 (*current_liboctave_error_handler) | |
442 ("Sparse::Sparse (const Array<T>&, const Array<double>&, ...): dimension mismatch"); | |
443 rep = nil_rep (); | |
444 dimensions = dim_vector (0, 0); | |
445 } | |
446 else | |
447 { | |
5604 | 448 octave_idx_type max_nzmx = (r_len > c_len ? r_len : c_len); |
5164 | 449 |
5604 | 450 OCTAVE_LOCAL_BUFFER (octave_sparse_sort_idxl *, sidx, max_nzmx); |
451 OCTAVE_LOCAL_BUFFER (octave_sparse_sort_idxl, sidxX, max_nzmx); | |
452 | |
453 for (octave_idx_type i = 0; i < max_nzmx; i++) | |
5164 | 454 sidx[i] = &sidxX[i]; |
455 | |
5604 | 456 octave_idx_type actual_nzmx = 0; |
5164 | 457 OCTAVE_QUIT; |
458 | |
5604 | 459 for (octave_idx_type i = 0; i < max_nzmx; i++) |
5164 | 460 { |
5275 | 461 octave_idx_type rowidx = static_cast<octave_idx_type> (ri_scalar ? r(0) : r(i)); |
462 octave_idx_type colidx = static_cast<octave_idx_type> (ci_scalar ? c(0) : c(i)); | |
5164 | 463 if (rowidx < nr && rowidx >= 0 && |
464 colidx < nc && colidx >= 0 ) | |
465 { | |
466 if ( a (cf_scalar ? 0 : i ) != T ()) | |
467 { | |
5604 | 468 sidx[actual_nzmx]->r = rowidx; |
469 sidx[actual_nzmx]->c = colidx; | |
470 sidx[actual_nzmx]->idx = i; | |
471 actual_nzmx++; | |
5164 | 472 } |
473 } | |
474 else | |
475 { | |
476 (*current_liboctave_error_handler) | |
477 ("Sparse::Sparse : index (%d,%d) out of range", | |
478 rowidx + 1, colidx + 1); | |
479 rep = nil_rep (); | |
480 dimensions = dim_vector (0, 0); | |
481 return; | |
482 } | |
483 } | |
484 | |
5604 | 485 if (actual_nzmx == 0) |
5164 | 486 rep = new typename Sparse<T>::SparseRep (nr, nc); |
487 else | |
488 { | |
489 OCTAVE_QUIT; | |
490 octave_sort<octave_sparse_sort_idxl *> | |
7433 | 491 lsort (octave_sparse_sidxl_comp); |
492 | |
493 lsort.sort (sidx, actual_nzmx); | |
5164 | 494 OCTAVE_QUIT; |
495 | |
496 // Now count the unique non-zero values | |
5604 | 497 octave_idx_type real_nzmx = 1; |
498 for (octave_idx_type i = 1; i < actual_nzmx; i++) | |
5164 | 499 if (sidx[i-1]->r != sidx[i]->r || sidx[i-1]->c != sidx[i]->c) |
5604 | 500 real_nzmx++; |
501 | |
502 rep = new typename Sparse<T>::SparseRep (nr, nc, real_nzmx); | |
5164 | 503 |
5275 | 504 octave_idx_type cx = 0; |
505 octave_idx_type prev_rval = -1; | |
506 octave_idx_type prev_cval = -1; | |
507 octave_idx_type ii = -1; | |
5164 | 508 xcidx (0) = 0; |
5604 | 509 for (octave_idx_type i = 0; i < actual_nzmx; i++) |
5164 | 510 { |
511 OCTAVE_QUIT; | |
5275 | 512 octave_idx_type iidx = sidx[i]->idx; |
513 octave_idx_type rval = sidx[i]->r; | |
514 octave_idx_type cval = sidx[i]->c; | |
5164 | 515 |
516 if (prev_cval < cval || (prev_rval < rval && prev_cval == cval)) | |
517 { | |
5275 | 518 octave_idx_type ci = static_cast<octave_idx_type> (c (ci_scalar ? 0 : iidx)); |
5164 | 519 ii++; |
520 | |
521 while (cx < ci) | |
522 xcidx (++cx) = ii; | |
523 xdata(ii) = a (cf_scalar ? 0 : iidx); | |
5275 | 524 xridx(ii) = static_cast<octave_idx_type> (r (ri_scalar ? 0 : iidx)); |
5164 | 525 } |
526 else | |
527 { | |
528 if (sum_terms) | |
529 xdata(ii) += a (cf_scalar ? 0 : iidx); | |
530 else | |
531 xdata(ii) = a (cf_scalar ? 0 : iidx); | |
532 } | |
533 prev_rval = rval; | |
534 prev_cval = cval; | |
535 } | |
536 | |
537 while (cx < nc) | |
538 xcidx (++cx) = ii + 1; | |
539 } | |
540 } | |
541 } | |
542 | |
543 template <class T> | |
544 Sparse<T>::Sparse (const Array2<T>& a) | |
545 : dimensions (a.dims ()), idx (0), idx_count (0) | |
546 { | |
5275 | 547 octave_idx_type nr = rows (); |
548 octave_idx_type nc = cols (); | |
549 octave_idx_type len = a.length (); | |
5604 | 550 octave_idx_type new_nzmx = 0; |
5164 | 551 |
552 // First count the number of non-zero terms | |
5275 | 553 for (octave_idx_type i = 0; i < len; i++) |
5164 | 554 if (a(i) != T ()) |
5604 | 555 new_nzmx++; |
556 | |
557 rep = new typename Sparse<T>::SparseRep (nr, nc, new_nzmx); | |
5164 | 558 |
5275 | 559 octave_idx_type ii = 0; |
5164 | 560 xcidx(0) = 0; |
5275 | 561 for (octave_idx_type j = 0; j < nc; j++) |
5164 | 562 { |
5275 | 563 for (octave_idx_type i = 0; i < nr; i++) |
5164 | 564 if (a.elem (i,j) != T ()) |
565 { | |
566 xdata(ii) = a.elem (i,j); | |
567 xridx(ii++) = i; | |
568 } | |
569 xcidx(j+1) = ii; | |
570 } | |
571 } | |
572 | |
573 template <class T> | |
574 Sparse<T>::Sparse (const Array<T>& a) | |
575 : dimensions (a.dims ()), idx (0), idx_count (0) | |
576 { | |
577 if (dimensions.length () > 2) | |
578 (*current_liboctave_error_handler) | |
579 ("Sparse::Sparse (const Array<T>&): dimension mismatch"); | |
580 else | |
581 { | |
5275 | 582 octave_idx_type nr = rows (); |
583 octave_idx_type nc = cols (); | |
584 octave_idx_type len = a.length (); | |
5604 | 585 octave_idx_type new_nzmx = 0; |
5164 | 586 |
587 // First count the number of non-zero terms | |
5275 | 588 for (octave_idx_type i = 0; i < len; i++) |
5164 | 589 if (a(i) != T ()) |
5604 | 590 new_nzmx++; |
591 | |
592 rep = new typename Sparse<T>::SparseRep (nr, nc, new_nzmx); | |
5164 | 593 |
5275 | 594 octave_idx_type ii = 0; |
5164 | 595 xcidx(0) = 0; |
5275 | 596 for (octave_idx_type j = 0; j < nc; j++) |
5164 | 597 { |
5275 | 598 for (octave_idx_type i = 0; i < nr; i++) |
5164 | 599 if (a.elem (i,j) != T ()) |
600 { | |
601 xdata(ii) = a.elem (i,j); | |
602 xridx(ii++) = i; | |
603 } | |
604 xcidx(j+1) = ii; | |
605 } | |
606 } | |
607 } | |
608 | |
609 template <class T> | |
610 Sparse<T>::~Sparse (void) | |
611 { | |
612 if (--rep->count <= 0) | |
613 delete rep; | |
614 | |
615 delete [] idx; | |
616 } | |
617 | |
618 template <class T> | |
7717
ff918ee1a983
Delete idx in Sparse<T> and Array<T> operator =
David Bateman <dbateman@free.fr>
parents:
7620
diff
changeset
|
619 Sparse<T>& |
ff918ee1a983
Delete idx in Sparse<T> and Array<T> operator =
David Bateman <dbateman@free.fr>
parents:
7620
diff
changeset
|
620 Sparse<T>::operator = (const Sparse<T>& a) |
ff918ee1a983
Delete idx in Sparse<T> and Array<T> operator =
David Bateman <dbateman@free.fr>
parents:
7620
diff
changeset
|
621 { |
ff918ee1a983
Delete idx in Sparse<T> and Array<T> operator =
David Bateman <dbateman@free.fr>
parents:
7620
diff
changeset
|
622 if (this != &a) |
ff918ee1a983
Delete idx in Sparse<T> and Array<T> operator =
David Bateman <dbateman@free.fr>
parents:
7620
diff
changeset
|
623 { |
ff918ee1a983
Delete idx in Sparse<T> and Array<T> operator =
David Bateman <dbateman@free.fr>
parents:
7620
diff
changeset
|
624 if (--rep->count <= 0) |
ff918ee1a983
Delete idx in Sparse<T> and Array<T> operator =
David Bateman <dbateman@free.fr>
parents:
7620
diff
changeset
|
625 delete rep; |
ff918ee1a983
Delete idx in Sparse<T> and Array<T> operator =
David Bateman <dbateman@free.fr>
parents:
7620
diff
changeset
|
626 |
ff918ee1a983
Delete idx in Sparse<T> and Array<T> operator =
David Bateman <dbateman@free.fr>
parents:
7620
diff
changeset
|
627 rep = a.rep; |
ff918ee1a983
Delete idx in Sparse<T> and Array<T> operator =
David Bateman <dbateman@free.fr>
parents:
7620
diff
changeset
|
628 rep->count++; |
ff918ee1a983
Delete idx in Sparse<T> and Array<T> operator =
David Bateman <dbateman@free.fr>
parents:
7620
diff
changeset
|
629 |
ff918ee1a983
Delete idx in Sparse<T> and Array<T> operator =
David Bateman <dbateman@free.fr>
parents:
7620
diff
changeset
|
630 dimensions = a.dimensions; |
ff918ee1a983
Delete idx in Sparse<T> and Array<T> operator =
David Bateman <dbateman@free.fr>
parents:
7620
diff
changeset
|
631 |
ff918ee1a983
Delete idx in Sparse<T> and Array<T> operator =
David Bateman <dbateman@free.fr>
parents:
7620
diff
changeset
|
632 delete [] idx; |
ff918ee1a983
Delete idx in Sparse<T> and Array<T> operator =
David Bateman <dbateman@free.fr>
parents:
7620
diff
changeset
|
633 idx_count = 0; |
ff918ee1a983
Delete idx in Sparse<T> and Array<T> operator =
David Bateman <dbateman@free.fr>
parents:
7620
diff
changeset
|
634 idx = 0; |
ff918ee1a983
Delete idx in Sparse<T> and Array<T> operator =
David Bateman <dbateman@free.fr>
parents:
7620
diff
changeset
|
635 } |
ff918ee1a983
Delete idx in Sparse<T> and Array<T> operator =
David Bateman <dbateman@free.fr>
parents:
7620
diff
changeset
|
636 |
ff918ee1a983
Delete idx in Sparse<T> and Array<T> operator =
David Bateman <dbateman@free.fr>
parents:
7620
diff
changeset
|
637 return *this; |
ff918ee1a983
Delete idx in Sparse<T> and Array<T> operator =
David Bateman <dbateman@free.fr>
parents:
7620
diff
changeset
|
638 } |
ff918ee1a983
Delete idx in Sparse<T> and Array<T> operator =
David Bateman <dbateman@free.fr>
parents:
7620
diff
changeset
|
639 |
ff918ee1a983
Delete idx in Sparse<T> and Array<T> operator =
David Bateman <dbateman@free.fr>
parents:
7620
diff
changeset
|
640 template <class T> |
5275 | 641 octave_idx_type |
642 Sparse<T>::compute_index (const Array<octave_idx_type>& ra_idx) const | |
5164 | 643 { |
5275 | 644 octave_idx_type retval = -1; |
645 | |
646 octave_idx_type n = dimensions.length (); | |
5164 | 647 |
648 if (n > 0 && n == ra_idx.length ()) | |
649 { | |
650 retval = ra_idx(--n); | |
651 | |
652 while (--n >= 0) | |
653 { | |
654 retval *= dimensions(n); | |
655 retval += ra_idx(n); | |
656 } | |
657 } | |
658 else | |
659 (*current_liboctave_error_handler) | |
660 ("Sparse<T>::compute_index: invalid ra_idxing operation"); | |
661 | |
662 return retval; | |
663 } | |
664 | |
665 template <class T> | |
666 T | |
5275 | 667 Sparse<T>::range_error (const char *fcn, octave_idx_type n) const |
5164 | 668 { |
669 (*current_liboctave_error_handler) ("%s (%d): range error", fcn, n); | |
670 return T (); | |
671 } | |
672 | |
673 template <class T> | |
674 T& | |
5275 | 675 Sparse<T>::range_error (const char *fcn, octave_idx_type n) |
5164 | 676 { |
677 (*current_liboctave_error_handler) ("%s (%d): range error", fcn, n); | |
678 static T foo; | |
679 return foo; | |
680 } | |
681 | |
682 template <class T> | |
683 T | |
5275 | 684 Sparse<T>::range_error (const char *fcn, octave_idx_type i, octave_idx_type j) const |
5164 | 685 { |
686 (*current_liboctave_error_handler) | |
687 ("%s (%d, %d): range error", fcn, i, j); | |
688 return T (); | |
689 } | |
690 | |
691 template <class T> | |
692 T& | |
5275 | 693 Sparse<T>::range_error (const char *fcn, octave_idx_type i, octave_idx_type j) |
5164 | 694 { |
695 (*current_liboctave_error_handler) | |
696 ("%s (%d, %d): range error", fcn, i, j); | |
697 static T foo; | |
698 return foo; | |
699 } | |
700 | |
701 template <class T> | |
702 T | |
5275 | 703 Sparse<T>::range_error (const char *fcn, const Array<octave_idx_type>& ra_idx) const |
5164 | 704 { |
5765 | 705 std::ostringstream buf; |
5164 | 706 |
707 buf << fcn << " ("; | |
708 | |
5275 | 709 octave_idx_type n = ra_idx.length (); |
5164 | 710 |
711 if (n > 0) | |
712 buf << ra_idx(0); | |
713 | |
5275 | 714 for (octave_idx_type i = 1; i < n; i++) |
5164 | 715 buf << ", " << ra_idx(i); |
716 | |
717 buf << "): range error"; | |
5765 | 718 |
719 std::string buf_str = buf.str (); | |
720 | |
721 (*current_liboctave_error_handler) (buf_str.c_str ()); | |
5164 | 722 |
723 return T (); | |
724 } | |
725 | |
726 template <class T> | |
727 T& | |
5275 | 728 Sparse<T>::range_error (const char *fcn, const Array<octave_idx_type>& ra_idx) |
5164 | 729 { |
5765 | 730 std::ostringstream buf; |
5164 | 731 |
732 buf << fcn << " ("; | |
733 | |
5275 | 734 octave_idx_type n = ra_idx.length (); |
5164 | 735 |
736 if (n > 0) | |
737 buf << ra_idx(0); | |
738 | |
5275 | 739 for (octave_idx_type i = 1; i < n; i++) |
5164 | 740 buf << ", " << ra_idx(i); |
741 | |
742 buf << "): range error"; | |
743 | |
5765 | 744 std::string buf_str = buf.str (); |
745 | |
746 (*current_liboctave_error_handler) (buf_str.c_str ()); | |
5164 | 747 |
748 static T foo; | |
749 return foo; | |
750 } | |
751 | |
752 template <class T> | |
753 Sparse<T> | |
754 Sparse<T>::reshape (const dim_vector& new_dims) const | |
755 { | |
756 Sparse<T> retval; | |
6689 | 757 dim_vector dims2 = new_dims; |
758 | |
759 if (dims2.length () > 2) | |
5164 | 760 { |
6814 | 761 (*current_liboctave_warning_handler) |
762 ("reshape: sparse reshape to N-d array smashes dims"); | |
763 | |
6689 | 764 for (octave_idx_type i = 2; i < dims2.length(); i++) |
6814 | 765 dims2(1) *= dims2(i); |
766 | |
6689 | 767 dims2.resize (2); |
768 } | |
769 | |
770 if (dimensions != dims2) | |
771 { | |
772 if (dimensions.numel () == dims2.numel ()) | |
5164 | 773 { |
5681 | 774 octave_idx_type new_nnz = nnz (); |
6689 | 775 octave_idx_type new_nr = dims2 (0); |
776 octave_idx_type new_nc = dims2 (1); | |
5275 | 777 octave_idx_type old_nr = rows (); |
778 octave_idx_type old_nc = cols (); | |
5681 | 779 retval = Sparse<T> (new_nr, new_nc, new_nnz); |
5164 | 780 |
5275 | 781 octave_idx_type kk = 0; |
5164 | 782 retval.xcidx(0) = 0; |
5275 | 783 for (octave_idx_type i = 0; i < old_nc; i++) |
784 for (octave_idx_type j = cidx(i); j < cidx(i+1); j++) | |
5164 | 785 { |
5275 | 786 octave_idx_type tmp = i * old_nr + ridx(j); |
787 octave_idx_type ii = tmp % new_nr; | |
788 octave_idx_type jj = (tmp - ii) / new_nr; | |
789 for (octave_idx_type k = kk; k < jj; k++) | |
5164 | 790 retval.xcidx(k+1) = j; |
791 kk = jj; | |
792 retval.xdata(j) = data(j); | |
793 retval.xridx(j) = ii; | |
794 } | |
5275 | 795 for (octave_idx_type k = kk; k < new_nc; k++) |
5681 | 796 retval.xcidx(k+1) = new_nnz; |
5164 | 797 } |
798 else | |
8527
6b074f37e8d7
reshape: improve error message
John W. Eaton <jwe@octave.org>
parents:
8526
diff
changeset
|
799 { |
6b074f37e8d7
reshape: improve error message
John W. Eaton <jwe@octave.org>
parents:
8526
diff
changeset
|
800 std::string dimensions_str = dimensions.str (); |
6b074f37e8d7
reshape: improve error message
John W. Eaton <jwe@octave.org>
parents:
8526
diff
changeset
|
801 std::string new_dims_str = new_dims.str (); |
6b074f37e8d7
reshape: improve error message
John W. Eaton <jwe@octave.org>
parents:
8526
diff
changeset
|
802 |
6b074f37e8d7
reshape: improve error message
John W. Eaton <jwe@octave.org>
parents:
8526
diff
changeset
|
803 (*current_liboctave_error_handler) |
6b074f37e8d7
reshape: improve error message
John W. Eaton <jwe@octave.org>
parents:
8526
diff
changeset
|
804 ("reshape: can't reshape %s array to %s array", |
6b074f37e8d7
reshape: improve error message
John W. Eaton <jwe@octave.org>
parents:
8526
diff
changeset
|
805 dimensions_str.c_str (), new_dims_str.c_str ()); |
6b074f37e8d7
reshape: improve error message
John W. Eaton <jwe@octave.org>
parents:
8526
diff
changeset
|
806 } |
5164 | 807 } |
808 else | |
809 retval = *this; | |
810 | |
811 return retval; | |
812 } | |
813 | |
814 template <class T> | |
815 Sparse<T> | |
5275 | 816 Sparse<T>::permute (const Array<octave_idx_type>& perm_vec, bool) const |
5164 | 817 { |
6813 | 818 // The only valid permutations of a sparse array are [1, 2] and [2, 1]. |
819 | |
820 bool fail = false; | |
6817 | 821 bool trans = false; |
6813 | 822 |
823 if (perm_vec.length () == 2) | |
5164 | 824 { |
6813 | 825 if (perm_vec(0) == 0 && perm_vec(1) == 1) |
826 /* do nothing */; | |
827 else if (perm_vec(0) == 1 && perm_vec(1) == 0) | |
6817 | 828 trans = true; |
5164 | 829 else |
6813 | 830 fail = true; |
5164 | 831 } |
832 else | |
6813 | 833 fail = true; |
834 | |
835 if (fail) | |
836 (*current_liboctave_error_handler) | |
837 ("permutation vector contains an invalid element"); | |
838 | |
6817 | 839 return trans ? this->transpose () : *this; |
5164 | 840 } |
841 | |
842 template <class T> | |
843 void | |
844 Sparse<T>::resize_no_fill (const dim_vector& dv) | |
845 { | |
5275 | 846 octave_idx_type n = dv.length (); |
5164 | 847 |
848 if (n != 2) | |
849 { | |
850 (*current_liboctave_error_handler) ("sparse array must be 2-D"); | |
851 return; | |
852 } | |
853 | |
854 resize_no_fill (dv(0), dv(1)); | |
855 } | |
856 | |
857 template <class T> | |
858 void | |
5275 | 859 Sparse<T>::resize_no_fill (octave_idx_type r, octave_idx_type c) |
5164 | 860 { |
861 if (r < 0 || c < 0) | |
862 { | |
863 (*current_liboctave_error_handler) | |
864 ("can't resize to negative dimension"); | |
865 return; | |
866 } | |
867 | |
868 if (ndims () == 0) | |
869 dimensions = dim_vector (0, 0); | |
870 | |
871 if (r == dim1 () && c == dim2 ()) | |
872 return; | |
873 | |
5731 | 874 typename Sparse<T>::SparseRep *old_rep = rep; |
875 | |
5275 | 876 octave_idx_type nc = cols (); |
877 octave_idx_type nr = rows (); | |
5164 | 878 |
5681 | 879 if (nnz () == 0 || r == 0 || c == 0) |
5164 | 880 // Special case of redimensioning to/from a sparse matrix with |
881 // no elements | |
882 rep = new typename Sparse<T>::SparseRep (r, c); | |
883 else | |
884 { | |
5275 | 885 octave_idx_type n = 0; |
5164 | 886 Sparse<T> tmpval; |
887 if (r >= nr) | |
888 { | |
889 if (c > nc) | |
5731 | 890 n = xcidx(nc); |
5164 | 891 else |
5731 | 892 n = xcidx(c); |
5164 | 893 |
894 tmpval = Sparse<T> (r, c, n); | |
895 | |
896 if (c > nc) | |
897 { | |
6101 | 898 for (octave_idx_type i = 0; i < nc + 1; i++) |
5731 | 899 tmpval.cidx(i) = xcidx(i); |
6101 | 900 for (octave_idx_type i = nc + 1; i < c + 1; i++) |
5164 | 901 tmpval.cidx(i) = tmpval.cidx(i-1); |
902 } | |
903 else if (c <= nc) | |
6101 | 904 for (octave_idx_type i = 0; i < c + 1; i++) |
5731 | 905 tmpval.cidx(i) = xcidx(i); |
5164 | 906 |
5275 | 907 for (octave_idx_type i = 0; i < n; i++) |
5164 | 908 { |
5731 | 909 tmpval.data(i) = xdata(i); |
910 tmpval.ridx(i) = xridx(i); | |
5164 | 911 } |
912 } | |
913 else | |
914 { | |
915 // Count how many non zero terms before we do anything | |
6101 | 916 octave_idx_type min_nc = (c < nc ? c : nc); |
917 for (octave_idx_type i = 0; i < min_nc; i++) | |
5731 | 918 for (octave_idx_type j = xcidx(i); j < xcidx(i+1); j++) |
919 if (xridx(j) < r) | |
5164 | 920 n++; |
921 | |
922 if (n) | |
923 { | |
924 // Now that we know the size we can do something | |
925 tmpval = Sparse<T> (r, c, n); | |
926 | |
927 tmpval.cidx(0); | |
6101 | 928 for (octave_idx_type i = 0, ii = 0; i < min_nc; i++) |
5164 | 929 { |
5731 | 930 for (octave_idx_type j = xcidx(i); j < xcidx(i+1); j++) |
931 if (xridx(j) < r) | |
5164 | 932 { |
5731 | 933 tmpval.data(ii) = xdata(j); |
934 tmpval.ridx(ii++) = xridx(j); | |
5164 | 935 } |
936 tmpval.cidx(i+1) = ii; | |
937 } | |
6101 | 938 if (c > min_nc) |
939 for (octave_idx_type i = nc; i < c; i++) | |
940 tmpval.cidx(i+1) = tmpval.cidx(i); | |
5164 | 941 } |
942 else | |
943 tmpval = Sparse<T> (r, c); | |
944 } | |
945 | |
946 rep = tmpval.rep; | |
947 rep->count++; | |
948 } | |
949 | |
950 dimensions = dim_vector (r, c); | |
951 | |
952 if (--old_rep->count <= 0) | |
953 delete old_rep; | |
954 } | |
955 | |
956 template <class T> | |
957 Sparse<T>& | |
5275 | 958 Sparse<T>::insert (const Sparse<T>& a, octave_idx_type r, octave_idx_type c) |
5164 | 959 { |
5275 | 960 octave_idx_type a_rows = a.rows (); |
961 octave_idx_type a_cols = a.cols (); | |
962 octave_idx_type nr = rows (); | |
963 octave_idx_type nc = cols (); | |
5164 | 964 |
965 if (r < 0 || r + a_rows > rows () || c < 0 || c + a_cols > cols ()) | |
966 { | |
967 (*current_liboctave_error_handler) ("range error for insert"); | |
968 return *this; | |
969 } | |
970 | |
971 // First count the number of elements in the final array | |
5681 | 972 octave_idx_type nel = cidx(c) + a.nnz (); |
5164 | 973 |
974 if (c + a_cols < nc) | |
975 nel += cidx(nc) - cidx(c + a_cols); | |
976 | |
5275 | 977 for (octave_idx_type i = c; i < c + a_cols; i++) |
978 for (octave_idx_type j = cidx(i); j < cidx(i+1); j++) | |
5164 | 979 if (ridx(j) < r || ridx(j) >= r + a_rows) |
980 nel++; | |
981 | |
982 Sparse<T> tmp (*this); | |
983 --rep->count; | |
984 rep = new typename Sparse<T>::SparseRep (nr, nc, nel); | |
985 | |
5275 | 986 for (octave_idx_type i = 0; i < tmp.cidx(c); i++) |
5164 | 987 { |
988 data(i) = tmp.data(i); | |
989 ridx(i) = tmp.ridx(i); | |
990 } | |
5275 | 991 for (octave_idx_type i = 0; i < c + 1; i++) |
5164 | 992 cidx(i) = tmp.cidx(i); |
993 | |
5275 | 994 octave_idx_type ii = cidx(c); |
995 | |
996 for (octave_idx_type i = c; i < c + a_cols; i++) | |
5164 | 997 { |
998 OCTAVE_QUIT; | |
999 | |
5275 | 1000 for (octave_idx_type j = tmp.cidx(i); j < tmp.cidx(i+1); j++) |
5164 | 1001 if (tmp.ridx(j) < r) |
1002 { | |
1003 data(ii) = tmp.data(j); | |
1004 ridx(ii++) = tmp.ridx(j); | |
1005 } | |
1006 | |
1007 OCTAVE_QUIT; | |
1008 | |
5275 | 1009 for (octave_idx_type j = a.cidx(i-c); j < a.cidx(i-c+1); j++) |
5164 | 1010 { |
1011 data(ii) = a.data(j); | |
1012 ridx(ii++) = r + a.ridx(j); | |
1013 } | |
1014 | |
1015 OCTAVE_QUIT; | |
1016 | |
5275 | 1017 for (octave_idx_type j = tmp.cidx(i); j < tmp.cidx(i+1); j++) |
5164 | 1018 if (tmp.ridx(j) >= r + a_rows) |
1019 { | |
1020 data(ii) = tmp.data(j); | |
1021 ridx(ii++) = tmp.ridx(j); | |
1022 } | |
1023 | |
1024 cidx(i+1) = ii; | |
1025 } | |
1026 | |
5275 | 1027 for (octave_idx_type i = c + a_cols; i < nc; i++) |
5164 | 1028 { |
5275 | 1029 for (octave_idx_type j = tmp.cidx(i); j < tmp.cidx(i+1); j++) |
5164 | 1030 { |
1031 data(ii) = tmp.data(j); | |
1032 ridx(ii++) = tmp.ridx(j); | |
1033 } | |
1034 cidx(i+1) = ii; | |
1035 } | |
1036 | |
1037 return *this; | |
1038 } | |
1039 | |
1040 template <class T> | |
1041 Sparse<T>& | |
5275 | 1042 Sparse<T>::insert (const Sparse<T>& a, const Array<octave_idx_type>& ra_idx) |
5164 | 1043 { |
1044 | |
1045 if (ra_idx.length () != 2) | |
1046 { | |
1047 (*current_liboctave_error_handler) ("range error for insert"); | |
1048 return *this; | |
1049 } | |
1050 | |
1051 return insert (a, ra_idx (0), ra_idx (1)); | |
1052 } | |
1053 | |
1054 template <class T> | |
1055 Sparse<T> | |
1056 Sparse<T>::transpose (void) const | |
1057 { | |
1058 assert (ndims () == 2); | |
1059 | |
5275 | 1060 octave_idx_type nr = rows (); |
1061 octave_idx_type nc = cols (); | |
5648 | 1062 octave_idx_type nz = nnz (); |
5164 | 1063 Sparse<T> retval (nc, nr, nz); |
1064 | |
5648 | 1065 OCTAVE_LOCAL_BUFFER (octave_idx_type, w, nr + 1); |
1066 for (octave_idx_type i = 0; i < nr; i++) | |
1067 w[i] = 0; | |
1068 for (octave_idx_type i = 0; i < nz; i++) | |
1069 w[ridx(i)]++; | |
1070 nz = 0; | |
1071 for (octave_idx_type i = 0; i < nr; i++) | |
5164 | 1072 { |
5648 | 1073 retval.xcidx(i) = nz; |
1074 nz += w[i]; | |
1075 w[i] = retval.xcidx(i); | |
5164 | 1076 } |
5648 | 1077 retval.xcidx(nr) = nz; |
1078 w[nr] = nz; | |
1079 | |
1080 for (octave_idx_type j = 0; j < nc; j++) | |
1081 for (octave_idx_type k = cidx(j); k < cidx(j+1); k++) | |
1082 { | |
1083 octave_idx_type q = w [ridx(k)]++; | |
1084 retval.xridx (q) = j; | |
1085 retval.xdata (q) = data (k); | |
1086 } | |
5164 | 1087 |
1088 return retval; | |
1089 } | |
1090 | |
1091 template <class T> | |
1092 void | |
1093 Sparse<T>::clear_index (void) | |
1094 { | |
1095 delete [] idx; | |
1096 idx = 0; | |
1097 idx_count = 0; | |
1098 } | |
1099 | |
1100 template <class T> | |
1101 void | |
1102 Sparse<T>::set_index (const idx_vector& idx_arg) | |
1103 { | |
5275 | 1104 octave_idx_type nd = ndims (); |
5164 | 1105 |
1106 if (! idx && nd > 0) | |
1107 idx = new idx_vector [nd]; | |
1108 | |
1109 if (idx_count < nd) | |
1110 { | |
1111 idx[idx_count++] = idx_arg; | |
1112 } | |
1113 else | |
1114 { | |
1115 idx_vector *new_idx = new idx_vector [idx_count+1]; | |
1116 | |
5275 | 1117 for (octave_idx_type i = 0; i < idx_count; i++) |
5164 | 1118 new_idx[i] = idx[i]; |
1119 | |
1120 new_idx[idx_count++] = idx_arg; | |
1121 | |
1122 delete [] idx; | |
1123 | |
1124 idx = new_idx; | |
1125 } | |
1126 } | |
1127 | |
1128 template <class T> | |
1129 void | |
1130 Sparse<T>::maybe_delete_elements (idx_vector& idx_arg) | |
1131 { | |
5275 | 1132 octave_idx_type nr = dim1 (); |
1133 octave_idx_type nc = dim2 (); | |
5164 | 1134 |
1135 if (nr == 0 && nc == 0) | |
1136 return; | |
1137 | |
5275 | 1138 octave_idx_type n; |
5164 | 1139 if (nr == 1) |
1140 n = nc; | |
1141 else if (nc == 1) | |
1142 n = nr; | |
1143 else | |
1144 { | |
1145 // Reshape to row vector for Matlab compatibility. | |
1146 | |
1147 n = nr * nc; | |
1148 nr = 1; | |
1149 nc = n; | |
1150 } | |
1151 | |
1152 if (idx_arg.is_colon_equiv (n, 1)) | |
1153 { | |
1154 // Either A(:) = [] or A(idx) = [] with idx enumerating all | |
1155 // elements, so we delete all elements and return [](0x0). To | |
1156 // preserve the orientation of the vector, you have to use | |
1157 // A(idx,:) = [] (delete rows) or A(:,idx) (delete columns). | |
1158 | |
1159 resize_no_fill (0, 0); | |
1160 return; | |
1161 } | |
1162 | |
1163 idx_arg.sort (true); | |
1164 | |
5275 | 1165 octave_idx_type num_to_delete = idx_arg.length (n); |
5164 | 1166 |
1167 if (num_to_delete != 0) | |
1168 { | |
5275 | 1169 octave_idx_type new_n = n; |
5681 | 1170 octave_idx_type new_nnz = nnz (); |
5275 | 1171 |
1172 octave_idx_type iidx = 0; | |
5164 | 1173 |
1174 const Sparse<T> tmp (*this); | |
1175 | |
5275 | 1176 for (octave_idx_type i = 0; i < n; i++) |
5164 | 1177 { |
1178 OCTAVE_QUIT; | |
1179 | |
1180 if (i == idx_arg.elem (iidx)) | |
1181 { | |
1182 iidx++; | |
1183 new_n--; | |
1184 | |
1185 if (tmp.elem (i) != T ()) | |
5681 | 1186 new_nnz--; |
5164 | 1187 |
1188 if (iidx == num_to_delete) | |
1189 break; | |
1190 } | |
1191 } | |
1192 | |
1193 if (new_n > 0) | |
1194 { | |
1195 rep->count--; | |
1196 | |
1197 if (nr == 1) | |
5681 | 1198 rep = new typename Sparse<T>::SparseRep (1, new_n, new_nnz); |
5164 | 1199 else |
5681 | 1200 rep = new typename Sparse<T>::SparseRep (new_n, 1, new_nnz); |
5164 | 1201 |
5275 | 1202 octave_idx_type ii = 0; |
1203 octave_idx_type jj = 0; | |
5164 | 1204 iidx = 0; |
5275 | 1205 for (octave_idx_type i = 0; i < n; i++) |
5164 | 1206 { |
1207 OCTAVE_QUIT; | |
1208 | |
1209 if (iidx < num_to_delete && i == idx_arg.elem (iidx)) | |
1210 iidx++; | |
1211 else | |
1212 { | |
1213 T el = tmp.elem (i); | |
1214 if (el != T ()) | |
1215 { | |
1216 data(ii) = el; | |
1217 ridx(ii++) = jj; | |
1218 } | |
1219 jj++; | |
1220 } | |
1221 } | |
1222 | |
1223 dimensions.resize (2); | |
1224 | |
1225 if (nr == 1) | |
1226 { | |
1227 ii = 0; | |
1228 cidx(0) = 0; | |
5275 | 1229 for (octave_idx_type i = 0; i < new_n; i++) |
5164 | 1230 { |
1231 OCTAVE_QUIT; | |
1232 if (ridx(ii) == i) | |
1233 ridx(ii++) = 0; | |
1234 cidx(i+1) = ii; | |
1235 } | |
1236 | |
1237 dimensions(0) = 1; | |
1238 dimensions(1) = new_n; | |
1239 } | |
1240 else | |
1241 { | |
1242 cidx(0) = 0; | |
5681 | 1243 cidx(1) = new_nnz; |
5164 | 1244 dimensions(0) = new_n; |
1245 dimensions(1) = 1; | |
1246 } | |
1247 } | |
1248 else | |
1249 (*current_liboctave_error_handler) | |
1250 ("A(idx) = []: index out of range"); | |
1251 } | |
1252 } | |
1253 | |
1254 template <class T> | |
1255 void | |
1256 Sparse<T>::maybe_delete_elements (idx_vector& idx_i, idx_vector& idx_j) | |
1257 { | |
1258 assert (ndims () == 2); | |
1259 | |
5275 | 1260 octave_idx_type nr = dim1 (); |
1261 octave_idx_type nc = dim2 (); | |
5164 | 1262 |
1263 if (nr == 0 && nc == 0) | |
1264 return; | |
1265 | |
1266 if (idx_i.is_colon ()) | |
1267 { | |
1268 if (idx_j.is_colon ()) | |
1269 { | |
1270 // A(:,:) -- We are deleting columns and rows, so the result | |
1271 // is [](0x0). | |
1272 | |
1273 resize_no_fill (0, 0); | |
1274 return; | |
1275 } | |
1276 | |
1277 if (idx_j.is_colon_equiv (nc, 1)) | |
1278 { | |
1279 // A(:,j) -- We are deleting columns by enumerating them, | |
1280 // If we enumerate all of them, we should have zero columns | |
1281 // with the same number of rows that we started with. | |
1282 | |
1283 resize_no_fill (nr, 0); | |
1284 return; | |
1285 } | |
1286 } | |
1287 | |
1288 if (idx_j.is_colon () && idx_i.is_colon_equiv (nr, 1)) | |
1289 { | |
1290 // A(i,:) -- We are deleting rows by enumerating them. If we | |
1291 // enumerate all of them, we should have zero rows with the | |
1292 // same number of columns that we started with. | |
1293 | |
1294 resize_no_fill (0, nc); | |
1295 return; | |
1296 } | |
1297 | |
1298 if (idx_i.is_colon_equiv (nr, 1)) | |
1299 { | |
1300 if (idx_j.is_colon_equiv (nc, 1)) | |
1301 resize_no_fill (0, 0); | |
1302 else | |
1303 { | |
1304 idx_j.sort (true); | |
1305 | |
5275 | 1306 octave_idx_type num_to_delete = idx_j.length (nc); |
5164 | 1307 |
1308 if (num_to_delete != 0) | |
1309 { | |
1310 if (nr == 1 && num_to_delete == nc) | |
1311 resize_no_fill (0, 0); | |
1312 else | |
1313 { | |
5275 | 1314 octave_idx_type new_nc = nc; |
5681 | 1315 octave_idx_type new_nnz = nnz (); |
5275 | 1316 |
1317 octave_idx_type iidx = 0; | |
1318 | |
1319 for (octave_idx_type j = 0; j < nc; j++) | |
5164 | 1320 { |
1321 OCTAVE_QUIT; | |
1322 | |
1323 if (j == idx_j.elem (iidx)) | |
1324 { | |
1325 iidx++; | |
1326 new_nc--; | |
1327 | |
5681 | 1328 new_nnz -= cidx(j+1) - cidx(j); |
5164 | 1329 |
1330 if (iidx == num_to_delete) | |
1331 break; | |
1332 } | |
1333 } | |
1334 | |
1335 if (new_nc > 0) | |
1336 { | |
1337 const Sparse<T> tmp (*this); | |
1338 --rep->count; | |
1339 rep = new typename Sparse<T>::SparseRep (nr, new_nc, | |
5681 | 1340 new_nnz); |
5275 | 1341 octave_idx_type ii = 0; |
1342 octave_idx_type jj = 0; | |
5164 | 1343 iidx = 0; |
1344 cidx(0) = 0; | |
5275 | 1345 for (octave_idx_type j = 0; j < nc; j++) |
5164 | 1346 { |
1347 OCTAVE_QUIT; | |
1348 | |
1349 if (iidx < num_to_delete && j == idx_j.elem (iidx)) | |
1350 iidx++; | |
1351 else | |
1352 { | |
5275 | 1353 for (octave_idx_type i = tmp.cidx(j); |
5164 | 1354 i < tmp.cidx(j+1); i++) |
1355 { | |
1356 data(jj) = tmp.data(i); | |
1357 ridx(jj++) = tmp.ridx(i); | |
1358 } | |
1359 cidx(++ii) = jj; | |
1360 } | |
1361 } | |
1362 | |
1363 dimensions.resize (2); | |
1364 dimensions(1) = new_nc; | |
1365 } | |
1366 else | |
1367 (*current_liboctave_error_handler) | |
1368 ("A(idx) = []: index out of range"); | |
1369 } | |
1370 } | |
1371 } | |
1372 } | |
1373 else if (idx_j.is_colon_equiv (nc, 1)) | |
1374 { | |
1375 if (idx_i.is_colon_equiv (nr, 1)) | |
1376 resize_no_fill (0, 0); | |
1377 else | |
1378 { | |
1379 idx_i.sort (true); | |
1380 | |
5275 | 1381 octave_idx_type num_to_delete = idx_i.length (nr); |
5164 | 1382 |
1383 if (num_to_delete != 0) | |
1384 { | |
1385 if (nc == 1 && num_to_delete == nr) | |
1386 resize_no_fill (0, 0); | |
1387 else | |
1388 { | |
5275 | 1389 octave_idx_type new_nr = nr; |
5681 | 1390 octave_idx_type new_nnz = nnz (); |
5275 | 1391 |
1392 octave_idx_type iidx = 0; | |
1393 | |
1394 for (octave_idx_type i = 0; i < nr; i++) | |
5164 | 1395 { |
1396 OCTAVE_QUIT; | |
1397 | |
1398 if (i == idx_i.elem (iidx)) | |
1399 { | |
1400 iidx++; | |
1401 new_nr--; | |
1402 | |
5681 | 1403 for (octave_idx_type j = 0; j < nnz (); j++) |
5164 | 1404 if (ridx(j) == i) |
5681 | 1405 new_nnz--; |
5164 | 1406 |
1407 if (iidx == num_to_delete) | |
1408 break; | |
1409 } | |
1410 } | |
1411 | |
1412 if (new_nr > 0) | |
1413 { | |
1414 const Sparse<T> tmp (*this); | |
1415 --rep->count; | |
1416 rep = new typename Sparse<T>::SparseRep (new_nr, nc, | |
5681 | 1417 new_nnz); |
5164 | 1418 |
5275 | 1419 octave_idx_type jj = 0; |
5164 | 1420 cidx(0) = 0; |
5275 | 1421 for (octave_idx_type i = 0; i < nc; i++) |
5164 | 1422 { |
1423 iidx = 0; | |
5275 | 1424 for (octave_idx_type j = tmp.cidx(i); j < tmp.cidx(i+1); j++) |
5164 | 1425 { |
1426 OCTAVE_QUIT; | |
1427 | |
5275 | 1428 octave_idx_type ri = tmp.ridx(j); |
5164 | 1429 |
1430 while (iidx < num_to_delete && | |
1431 ri > idx_i.elem (iidx)) | |
1432 { | |
1433 iidx++; | |
1434 } | |
1435 | |
1436 if (iidx == num_to_delete || | |
1437 ri != idx_i.elem(iidx)) | |
1438 { | |
1439 data(jj) = tmp.data(j); | |
1440 ridx(jj++) = ri - iidx; | |
1441 } | |
1442 } | |
1443 cidx(i+1) = jj; | |
1444 } | |
1445 | |
1446 dimensions.resize (2); | |
1447 dimensions(0) = new_nr; | |
1448 } | |
1449 else | |
1450 (*current_liboctave_error_handler) | |
1451 ("A(idx) = []: index out of range"); | |
1452 } | |
1453 } | |
1454 } | |
1455 } | |
1456 } | |
1457 | |
1458 template <class T> | |
1459 void | |
1460 Sparse<T>::maybe_delete_elements (Array<idx_vector>& ra_idx) | |
1461 { | |
1462 if (ra_idx.length () == 1) | |
1463 maybe_delete_elements (ra_idx(0)); | |
1464 else if (ra_idx.length () == 2) | |
1465 maybe_delete_elements (ra_idx(0), ra_idx(1)); | |
1466 else | |
1467 (*current_liboctave_error_handler) | |
1468 ("range error for maybe_delete_elements"); | |
1469 } | |
1470 | |
1471 template <class T> | |
1472 Sparse<T> | |
1473 Sparse<T>::value (void) | |
1474 { | |
1475 Sparse<T> retval; | |
1476 | |
1477 int n_idx = index_count (); | |
1478 | |
1479 if (n_idx == 2) | |
1480 { | |
1481 idx_vector *tmp = get_idx (); | |
1482 | |
1483 idx_vector idx_i = tmp[0]; | |
1484 idx_vector idx_j = tmp[1]; | |
1485 | |
1486 retval = index (idx_i, idx_j); | |
1487 } | |
1488 else if (n_idx == 1) | |
1489 { | |
1490 retval = index (idx[0]); | |
1491 } | |
1492 else | |
1493 (*current_liboctave_error_handler) | |
1494 ("Sparse<T>::value: invalid number of indices specified"); | |
1495 | |
1496 clear_index (); | |
1497 | |
1498 return retval; | |
1499 } | |
1500 | |
1501 template <class T> | |
1502 Sparse<T> | |
1503 Sparse<T>::index (idx_vector& idx_arg, int resize_ok) const | |
1504 { | |
1505 Sparse<T> retval; | |
1506 | |
1507 assert (ndims () == 2); | |
1508 | |
5275 | 1509 octave_idx_type nr = dim1 (); |
1510 octave_idx_type nc = dim2 (); | |
5681 | 1511 octave_idx_type nz = nnz (); |
5275 | 1512 |
1513 octave_idx_type orig_len = nr * nc; | |
5164 | 1514 |
1515 dim_vector idx_orig_dims = idx_arg.orig_dimensions (); | |
1516 | |
5275 | 1517 octave_idx_type idx_orig_rows = idx_arg.orig_rows (); |
1518 octave_idx_type idx_orig_columns = idx_arg.orig_columns (); | |
5164 | 1519 |
1520 if (idx_orig_dims.length () > 2) | |
1521 (*current_liboctave_error_handler) | |
1522 ("Sparse<T>::index: Can not index Sparse<T> with an N-D Array"); | |
1523 else if (idx_arg.is_colon ()) | |
1524 { | |
1525 // Fast magic colon processing. | |
1526 retval = Sparse<T> (nr * nc, 1, nz); | |
1527 | |
5275 | 1528 for (octave_idx_type i = 0; i < nc; i++) |
1529 for (octave_idx_type j = cidx(i); j < cidx(i+1); j++) | |
5164 | 1530 { |
1531 OCTAVE_QUIT; | |
1532 retval.xdata(j) = data(j); | |
1533 retval.xridx(j) = ridx(j) + i * nr; | |
1534 } | |
1535 retval.xcidx(0) = 0; | |
1536 retval.xcidx(1) = nz; | |
1537 } | |
1538 else if (nr == 1 && nc == 1) | |
1539 { | |
1540 // You have to be pretty sick to get to this bit of code, | |
1541 // since you have a scalar stored as a sparse matrix, and | |
1542 // then want to make a dense matrix with sparse | |
1543 // representation. Ok, we'll do it, but you deserve what | |
1544 // you get!! | |
5275 | 1545 octave_idx_type n = idx_arg.freeze (length (), "sparse vector", resize_ok); |
5164 | 1546 if (n == 0) |
7573
755bf7ecc29b
eliminate one_zero stuff from idx_vector
John W. Eaton <jwe@octave.org>
parents:
7546
diff
changeset
|
1547 |
7299 | 1548 retval = Sparse<T> (idx_orig_dims); |
5164 | 1549 else if (nz < 1) |
1550 if (n >= idx_orig_dims.numel ()) | |
1551 retval = Sparse<T> (idx_orig_dims); | |
1552 else | |
1553 retval = Sparse<T> (dim_vector (n, 1)); | |
1554 else if (n >= idx_orig_dims.numel ()) | |
1555 { | |
1556 T el = elem (0); | |
5275 | 1557 octave_idx_type new_nr = idx_orig_rows; |
1558 octave_idx_type new_nc = idx_orig_columns; | |
1559 for (octave_idx_type i = 2; i < idx_orig_dims.length (); i++) | |
5164 | 1560 new_nc *= idx_orig_dims (i); |
1561 | |
1562 retval = Sparse<T> (new_nr, new_nc, idx_arg.ones_count ()); | |
1563 | |
5275 | 1564 octave_idx_type ic = 0; |
1565 for (octave_idx_type i = 0; i < n; i++) | |
5164 | 1566 { |
1567 if (i % new_nr == 0) | |
7322 | 1568 retval.xcidx(i / new_nr) = ic; |
5164 | 1569 |
5275 | 1570 octave_idx_type ii = idx_arg.elem (i); |
5164 | 1571 if (ii == 0) |
1572 { | |
1573 OCTAVE_QUIT; | |
1574 retval.xdata(ic) = el; | |
1575 retval.xridx(ic++) = i % new_nr; | |
1576 } | |
1577 } | |
1578 retval.xcidx (new_nc) = ic; | |
1579 } | |
1580 else | |
1581 { | |
1582 T el = elem (0); | |
1583 retval = Sparse<T> (n, 1, nz); | |
1584 | |
5275 | 1585 for (octave_idx_type i = 0; i < nz; i++) |
5164 | 1586 { |
1587 OCTAVE_QUIT; | |
1588 retval.xdata(i) = el; | |
1589 retval.xridx(i) = i; | |
1590 } | |
1591 retval.xcidx(0) = 0; | |
1592 retval.xcidx(1) = n; | |
1593 } | |
1594 } | |
1595 else if (nr == 1 || nc == 1) | |
1596 { | |
1597 // If indexing a vector with a matrix, return value has same | |
1598 // shape as the index. Otherwise, it has same orientation as | |
1599 // indexed object. | |
5275 | 1600 octave_idx_type len = length (); |
1601 octave_idx_type n = idx_arg.freeze (len, "sparse vector", resize_ok); | |
5164 | 1602 |
1603 if (n == 0) | |
1604 if (nr == 1) | |
1605 retval = Sparse<T> (dim_vector (1, 0)); | |
1606 else | |
1607 retval = Sparse<T> (dim_vector (0, 1)); | |
1608 else if (nz < 1) | |
7573
755bf7ecc29b
eliminate one_zero stuff from idx_vector
John W. Eaton <jwe@octave.org>
parents:
7546
diff
changeset
|
1609 if (idx_orig_rows == 1 || idx_orig_columns == 1) |
5164 | 1610 retval = Sparse<T> ((nr == 1 ? 1 : n), (nr == 1 ? n : 1)); |
1611 else | |
1612 retval = Sparse<T> (idx_orig_dims); | |
1613 else | |
1614 { | |
1615 | |
5604 | 1616 octave_idx_type new_nzmx = 0; |
5164 | 1617 if (nr == 1) |
5275 | 1618 for (octave_idx_type i = 0; i < n; i++) |
5164 | 1619 { |
1620 OCTAVE_QUIT; | |
1621 | |
5275 | 1622 octave_idx_type ii = idx_arg.elem (i); |
5164 | 1623 if (ii < len) |
1624 if (cidx(ii) != cidx(ii+1)) | |
5604 | 1625 new_nzmx++; |
5164 | 1626 } |
1627 else | |
5275 | 1628 for (octave_idx_type i = 0; i < n; i++) |
5164 | 1629 { |
5275 | 1630 octave_idx_type ii = idx_arg.elem (i); |
5164 | 1631 if (ii < len) |
5275 | 1632 for (octave_idx_type j = 0; j < nz; j++) |
5164 | 1633 { |
1634 OCTAVE_QUIT; | |
1635 | |
1636 if (ridx(j) == ii) | |
5604 | 1637 new_nzmx++; |
5164 | 1638 if (ridx(j) >= ii) |
1639 break; | |
1640 } | |
1641 } | |
1642 | |
7573
755bf7ecc29b
eliminate one_zero stuff from idx_vector
John W. Eaton <jwe@octave.org>
parents:
7546
diff
changeset
|
1643 if (idx_orig_rows == 1 || idx_orig_columns == 1) |
5164 | 1644 { |
1645 if (nr == 1) | |
1646 { | |
5604 | 1647 retval = Sparse<T> (1, n, new_nzmx); |
5275 | 1648 octave_idx_type jj = 0; |
5164 | 1649 retval.xcidx(0) = 0; |
5275 | 1650 for (octave_idx_type i = 0; i < n; i++) |
5164 | 1651 { |
1652 OCTAVE_QUIT; | |
1653 | |
5275 | 1654 octave_idx_type ii = idx_arg.elem (i); |
5164 | 1655 if (ii < len) |
1656 if (cidx(ii) != cidx(ii+1)) | |
1657 { | |
1658 retval.xdata(jj) = data(cidx(ii)); | |
1659 retval.xridx(jj++) = 0; | |
1660 } | |
1661 retval.xcidx(i+1) = jj; | |
1662 } | |
1663 } | |
1664 else | |
1665 { | |
5604 | 1666 retval = Sparse<T> (n, 1, new_nzmx); |
5164 | 1667 retval.xcidx(0) = 0; |
5604 | 1668 retval.xcidx(1) = new_nzmx; |
5275 | 1669 octave_idx_type jj = 0; |
1670 for (octave_idx_type i = 0; i < n; i++) | |
5164 | 1671 { |
5275 | 1672 octave_idx_type ii = idx_arg.elem (i); |
5164 | 1673 if (ii < len) |
5275 | 1674 for (octave_idx_type j = 0; j < nz; j++) |
5164 | 1675 { |
1676 OCTAVE_QUIT; | |
1677 | |
1678 if (ridx(j) == ii) | |
1679 { | |
1680 retval.xdata(jj) = data(j); | |
1681 retval.xridx(jj++) = i; | |
1682 } | |
1683 if (ridx(j) >= ii) | |
1684 break; | |
1685 } | |
1686 } | |
1687 } | |
1688 } | |
1689 else | |
1690 { | |
5275 | 1691 octave_idx_type new_nr; |
1692 octave_idx_type new_nc; | |
5164 | 1693 if (n >= idx_orig_dims.numel ()) |
1694 { | |
1695 new_nr = idx_orig_rows; | |
1696 new_nc = idx_orig_columns; | |
1697 } | |
1698 else | |
1699 { | |
1700 new_nr = n; | |
1701 new_nc = 1; | |
1702 } | |
1703 | |
5604 | 1704 retval = Sparse<T> (new_nr, new_nc, new_nzmx); |
5164 | 1705 |
1706 if (nr == 1) | |
1707 { | |
5275 | 1708 octave_idx_type jj = 0; |
5164 | 1709 retval.xcidx(0) = 0; |
5275 | 1710 for (octave_idx_type i = 0; i < n; i++) |
5164 | 1711 { |
1712 OCTAVE_QUIT; | |
1713 | |
5275 | 1714 octave_idx_type ii = idx_arg.elem (i); |
5164 | 1715 if (ii < len) |
1716 if (cidx(ii) != cidx(ii+1)) | |
1717 { | |
1718 retval.xdata(jj) = data(cidx(ii)); | |
1719 retval.xridx(jj++) = 0; | |
1720 } | |
1721 retval.xcidx(i/new_nr+1) = jj; | |
1722 } | |
1723 } | |
1724 else | |
1725 { | |
5275 | 1726 octave_idx_type jj = 0; |
5164 | 1727 retval.xcidx(0) = 0; |
5275 | 1728 for (octave_idx_type i = 0; i < n; i++) |
5164 | 1729 { |
5275 | 1730 octave_idx_type ii = idx_arg.elem (i); |
5164 | 1731 if (ii < len) |
5275 | 1732 for (octave_idx_type j = 0; j < nz; j++) |
5164 | 1733 { |
1734 OCTAVE_QUIT; | |
1735 | |
1736 if (ridx(j) == ii) | |
1737 { | |
1738 retval.xdata(jj) = data(j); | |
1739 retval.xridx(jj++) = i; | |
1740 } | |
1741 if (ridx(j) >= ii) | |
1742 break; | |
1743 } | |
1744 retval.xcidx(i/new_nr+1) = jj; | |
1745 } | |
1746 } | |
1747 } | |
1748 } | |
1749 } | |
1750 else | |
1751 { | |
7573
755bf7ecc29b
eliminate one_zero stuff from idx_vector
John W. Eaton <jwe@octave.org>
parents:
7546
diff
changeset
|
1752 (*current_liboctave_warning_with_id_handler) |
755bf7ecc29b
eliminate one_zero stuff from idx_vector
John W. Eaton <jwe@octave.org>
parents:
7546
diff
changeset
|
1753 ("Octave:fortran-indexing", "single index used for sparse matrix"); |
5164 | 1754 |
1755 // This code is only for indexing matrices. The vector | |
1756 // cases are handled above. | |
1757 | |
1758 idx_arg.freeze (nr * nc, "matrix", resize_ok); | |
1759 | |
1760 if (idx_arg) | |
1761 { | |
5275 | 1762 octave_idx_type result_nr = idx_orig_rows; |
1763 octave_idx_type result_nc = idx_orig_columns; | |
5164 | 1764 |
1765 if (nz < 1) | |
7573
755bf7ecc29b
eliminate one_zero stuff from idx_vector
John W. Eaton <jwe@octave.org>
parents:
7546
diff
changeset
|
1766 retval = Sparse<T> (result_nr, result_nc); |
5164 | 1767 else |
1768 { | |
1769 // Count number of non-zero elements | |
5604 | 1770 octave_idx_type new_nzmx = 0; |
5275 | 1771 octave_idx_type kk = 0; |
1772 for (octave_idx_type j = 0; j < result_nc; j++) | |
5164 | 1773 { |
5275 | 1774 for (octave_idx_type i = 0; i < result_nr; i++) |
5164 | 1775 { |
1776 OCTAVE_QUIT; | |
1777 | |
5275 | 1778 octave_idx_type ii = idx_arg.elem (kk++); |
5164 | 1779 if (ii < orig_len) |
1780 { | |
5275 | 1781 octave_idx_type fr = ii % nr; |
1782 octave_idx_type fc = (ii - fr) / nr; | |
1783 for (octave_idx_type k = cidx(fc); k < cidx(fc+1); k++) | |
5164 | 1784 { |
1785 if (ridx(k) == fr) | |
5604 | 1786 new_nzmx++; |
5164 | 1787 if (ridx(k) >= fr) |
1788 break; | |
1789 } | |
1790 } | |
1791 } | |
1792 } | |
1793 | |
5604 | 1794 retval = Sparse<T> (result_nr, result_nc, new_nzmx); |
5164 | 1795 |
1796 kk = 0; | |
5275 | 1797 octave_idx_type jj = 0; |
5164 | 1798 retval.xcidx(0) = 0; |
5275 | 1799 for (octave_idx_type j = 0; j < result_nc; j++) |
5164 | 1800 { |
5275 | 1801 for (octave_idx_type i = 0; i < result_nr; i++) |
5164 | 1802 { |
1803 OCTAVE_QUIT; | |
1804 | |
5275 | 1805 octave_idx_type ii = idx_arg.elem (kk++); |
5164 | 1806 if (ii < orig_len) |
1807 { | |
5275 | 1808 octave_idx_type fr = ii % nr; |
1809 octave_idx_type fc = (ii - fr) / nr; | |
1810 for (octave_idx_type k = cidx(fc); k < cidx(fc+1); k++) | |
5164 | 1811 { |
1812 if (ridx(k) == fr) | |
1813 { | |
1814 retval.xdata(jj) = data(k); | |
1815 retval.xridx(jj++) = i; | |
1816 } | |
1817 if (ridx(k) >= fr) | |
1818 break; | |
1819 } | |
1820 } | |
1821 } | |
1822 retval.xcidx(j+1) = jj; | |
1823 } | |
1824 } | |
1825 // idx_vector::freeze() printed an error message for us. | |
1826 } | |
1827 } | |
1828 | |
1829 return retval; | |
1830 } | |
1831 | |
6988 | 1832 struct |
1833 idx_node | |
1834 { | |
1835 octave_idx_type i; | |
1836 struct idx_node *next; | |
1837 }; | |
1838 | |
5164 | 1839 template <class T> |
1840 Sparse<T> | |
1841 Sparse<T>::index (idx_vector& idx_i, idx_vector& idx_j, int resize_ok) const | |
1842 { | |
1843 Sparse<T> retval; | |
1844 | |
1845 assert (ndims () == 2); | |
1846 | |
5275 | 1847 octave_idx_type nr = dim1 (); |
1848 octave_idx_type nc = dim2 (); | |
1849 | |
1850 octave_idx_type n = idx_i.freeze (nr, "row", resize_ok); | |
1851 octave_idx_type m = idx_j.freeze (nc, "column", resize_ok); | |
5164 | 1852 |
1853 if (idx_i && idx_j) | |
1854 { | |
1855 if (idx_i.orig_empty () || idx_j.orig_empty () || n == 0 || m == 0) | |
1856 { | |
1857 retval.resize_no_fill (n, m); | |
1858 } | |
5681 | 1859 else |
5164 | 1860 { |
5681 | 1861 int idx_i_colon = idx_i.is_colon_equiv (nr); |
1862 int idx_j_colon = idx_j.is_colon_equiv (nc); | |
1863 | |
1864 if (idx_i_colon && idx_j_colon) | |
1865 { | |
1866 retval = *this; | |
1867 } | |
1868 else | |
5164 | 1869 { |
5681 | 1870 // Identify if the indices have any repeated values |
1871 bool permutation = true; | |
1872 | |
1873 OCTAVE_LOCAL_BUFFER (octave_idx_type, itmp, | |
1874 (nr > nc ? nr : nc)); | |
7433 | 1875 octave_sort<octave_idx_type> lsort; |
5681 | 1876 |
1877 if (n > nr || m > nc) | |
1878 permutation = false; | |
1879 | |
1880 if (permutation && ! idx_i_colon) | |
1881 { | |
1882 // Can't use something like | |
1883 // idx_vector tmp_idx = idx_i; | |
1884 // tmp_idx.sort (true); | |
1885 // if (tmp_idx.length(nr) != n) | |
1886 // permutation = false; | |
1887 // here as there is no make_unique function | |
1888 // for idx_vector type. | |
1889 for (octave_idx_type i = 0; i < n; i++) | |
1890 itmp [i] = idx_i.elem (i); | |
7433 | 1891 lsort.sort (itmp, n); |
5681 | 1892 for (octave_idx_type i = 1; i < n; i++) |
1893 if (itmp[i-1] == itmp[i]) | |
1894 { | |
1895 permutation = false; | |
1896 break; | |
1897 } | |
1898 } | |
1899 if (permutation && ! idx_j_colon) | |
1900 { | |
1901 for (octave_idx_type i = 0; i < m; i++) | |
1902 itmp [i] = idx_j.elem (i); | |
7433 | 1903 lsort.sort (itmp, m); |
5681 | 1904 for (octave_idx_type i = 1; i < m; i++) |
1905 if (itmp[i-1] == itmp[i]) | |
1906 { | |
1907 permutation = false; | |
1908 break; | |
1909 } | |
1910 } | |
1911 | |
1912 if (permutation) | |
5164 | 1913 { |
5681 | 1914 // Special case permutation like indexing for speed |
1915 retval = Sparse<T> (n, m, nnz ()); | |
1916 octave_idx_type *ri = retval.xridx (); | |
1917 | |
5766 | 1918 std::vector<T> X (n); |
5681 | 1919 for (octave_idx_type i = 0; i < nr; i++) |
1920 itmp [i] = -1; | |
1921 for (octave_idx_type i = 0; i < n; i++) | |
1922 itmp[idx_i.elem(i)] = i; | |
1923 | |
1924 octave_idx_type kk = 0; | |
1925 retval.xcidx(0) = 0; | |
1926 for (octave_idx_type j = 0; j < m; j++) | |
1927 { | |
1928 octave_idx_type jj = idx_j.elem (j); | |
1929 for (octave_idx_type i = cidx(jj); i < cidx(jj+1); i++) | |
1930 { | |
6988 | 1931 OCTAVE_QUIT; |
1932 | |
5681 | 1933 octave_idx_type ii = itmp [ridx(i)]; |
1934 if (ii >= 0) | |
1935 { | |
1936 X [ii] = data (i); | |
1937 retval.xridx (kk++) = ii; | |
1938 } | |
1939 } | |
7433 | 1940 lsort.sort (ri + retval.xcidx (j), kk - retval.xcidx (j)); |
5681 | 1941 for (octave_idx_type p = retval.xcidx (j); p < kk; p++) |
1942 retval.xdata (p) = X [retval.xridx (p)]; | |
1943 retval.xcidx(j+1) = kk; | |
1944 } | |
1945 retval.maybe_compress (); | |
1946 } | |
1947 else | |
1948 { | |
6988 | 1949 OCTAVE_LOCAL_BUFFER (struct idx_node, nodes, n); |
1950 OCTAVE_LOCAL_BUFFER (octave_idx_type, start_nodes, nr); | |
1951 | |
1952 for (octave_idx_type i = 0; i < nr; i++) | |
1953 start_nodes[i] = -1; | |
1954 | |
1955 for (octave_idx_type i = 0; i < n; i++) | |
1956 { | |
1957 octave_idx_type ii = idx_i.elem (i); | |
1958 nodes[i].i = i; | |
1959 nodes[i].next = 0; | |
1960 | |
1961 octave_idx_type node = start_nodes[ii]; | |
1962 if (node == -1) | |
1963 start_nodes[ii] = i; | |
1964 else | |
1965 { | |
7322 | 1966 while (nodes[node].next) |
1967 node = nodes[node].next->i; | |
1968 nodes[node].next = nodes + i; | |
6988 | 1969 } |
1970 } | |
1971 | |
5681 | 1972 // First count the number of non-zero elements |
1973 octave_idx_type new_nzmx = 0; | |
1974 for (octave_idx_type j = 0; j < m; j++) | |
5164 | 1975 { |
5681 | 1976 octave_idx_type jj = idx_j.elem (j); |
6988 | 1977 |
1978 if (jj < nc) | |
5164 | 1979 { |
6988 | 1980 for (octave_idx_type i = cidx(jj); |
1981 i < cidx(jj+1); i++) | |
5681 | 1982 { |
6988 | 1983 OCTAVE_QUIT; |
1984 | |
1985 octave_idx_type ii = start_nodes [ridx(i)]; | |
1986 | |
1987 if (ii >= 0) | |
5681 | 1988 { |
6988 | 1989 struct idx_node inode = nodes[ii]; |
1990 | |
1991 while (true) | |
1992 { | |
7326 | 1993 if (idx_i.elem (inode.i) < nr) |
6988 | 1994 new_nzmx ++; |
1995 if (inode.next == 0) | |
1996 break; | |
1997 else | |
1998 inode = *inode.next; | |
1999 } | |
5681 | 2000 } |
2001 } | |
5164 | 2002 } |
2003 } | |
5681 | 2004 |
6988 | 2005 std::vector<T> X (n); |
5681 | 2006 retval = Sparse<T> (n, m, new_nzmx); |
6988 | 2007 octave_idx_type *ri = retval.xridx (); |
5681 | 2008 |
2009 octave_idx_type kk = 0; | |
2010 retval.xcidx(0) = 0; | |
2011 for (octave_idx_type j = 0; j < m; j++) | |
2012 { | |
2013 octave_idx_type jj = idx_j.elem (j); | |
6988 | 2014 if (jj < nc) |
5681 | 2015 { |
6988 | 2016 for (octave_idx_type i = cidx(jj); |
2017 i < cidx(jj+1); i++) | |
5681 | 2018 { |
6988 | 2019 OCTAVE_QUIT; |
2020 | |
2021 octave_idx_type ii = start_nodes [ridx(i)]; | |
2022 | |
2023 if (ii >= 0) | |
5681 | 2024 { |
6988 | 2025 struct idx_node inode = nodes[ii]; |
2026 | |
2027 while (true) | |
5681 | 2028 { |
7326 | 2029 if (idx_i.elem (inode.i) < nr) |
6988 | 2030 { |
2031 X [inode.i] = data (i); | |
2032 retval.xridx (kk++) = inode.i; | |
2033 } | |
2034 | |
2035 if (inode.next == 0) | |
2036 break; | |
2037 else | |
2038 inode = *inode.next; | |
5681 | 2039 } |
2040 } | |
2041 } | |
7433 | 2042 lsort.sort (ri + retval.xcidx (j), |
6988 | 2043 kk - retval.xcidx (j)); |
2044 for (octave_idx_type p = retval.xcidx (j); | |
2045 p < kk; p++) | |
2046 retval.xdata (p) = X [retval.xridx (p)]; | |
2047 retval.xcidx(j+1) = kk; | |
5681 | 2048 } |
2049 } | |
5164 | 2050 } |
2051 } | |
2052 } | |
2053 } | |
2054 // idx_vector::freeze() printed an error message for us. | |
2055 | |
2056 return retval; | |
2057 } | |
2058 | |
2059 template <class T> | |
2060 Sparse<T> | |
2061 Sparse<T>::index (Array<idx_vector>& ra_idx, int resize_ok) const | |
2062 { | |
2063 | |
2064 if (ra_idx.length () != 2) | |
2065 { | |
2066 (*current_liboctave_error_handler) ("range error for index"); | |
2067 return *this; | |
2068 } | |
2069 | |
2070 return index (ra_idx (0), ra_idx (1), resize_ok); | |
2071 } | |
2072 | |
7433 | 2073 // Can't use versions of these in Array.cc due to duplication of the |
2074 // instantiations for Array<double and Sparse<double>, etc | |
2075 template <class T> | |
2076 bool | |
2077 sparse_ascending_compare (T a, T b) | |
2078 { | |
2079 return (a < b); | |
2080 } | |
2081 | |
2082 template <class T> | |
2083 bool | |
2084 sparse_descending_compare (T a, T b) | |
2085 { | |
2086 return (a > b); | |
2087 } | |
2088 | |
2089 template <class T> | |
2090 bool | |
2091 sparse_ascending_compare (vec_index<T> *a, vec_index<T> *b) | |
2092 { | |
2093 return (a->vec < b->vec); | |
2094 } | |
2095 | |
2096 template <class T> | |
2097 bool | |
2098 sparse_descending_compare (vec_index<T> *a, vec_index<T> *b) | |
2099 { | |
2100 return (a->vec > b->vec); | |
2101 } | |
2102 | |
2103 template <class T> | |
2104 Sparse<T> | |
2105 Sparse<T>::sort (octave_idx_type dim, sortmode mode) const | |
2106 { | |
2107 Sparse<T> m = *this; | |
2108 | |
2109 octave_idx_type nr = m.rows (); | |
2110 octave_idx_type nc = m.columns (); | |
2111 | |
2112 if (m.length () < 1) | |
2113 return m; | |
2114 | |
2115 if (dim > 0) | |
2116 { | |
2117 m = m.transpose (); | |
2118 nr = m.rows (); | |
2119 nc = m.columns (); | |
2120 } | |
2121 | |
2122 octave_sort<T> lsort; | |
2123 | |
2124 if (mode == ASCENDING) | |
2125 lsort.set_compare (sparse_ascending_compare); | |
2126 else if (mode == DESCENDING) | |
2127 lsort.set_compare (sparse_descending_compare); | |
7463
2467639bd8c0
eliminate UNDEFINED sort mode
John W. Eaton <jwe@octave.org>
parents:
7433
diff
changeset
|
2128 else |
2467639bd8c0
eliminate UNDEFINED sort mode
John W. Eaton <jwe@octave.org>
parents:
7433
diff
changeset
|
2129 abort (); |
7433 | 2130 |
2131 T *v = m.data (); | |
2132 octave_idx_type *mcidx = m.cidx (); | |
2133 octave_idx_type *mridx = m.ridx (); | |
2134 | |
2135 for (octave_idx_type j = 0; j < nc; j++) | |
2136 { | |
2137 octave_idx_type ns = mcidx [j + 1] - mcidx [j]; | |
2138 lsort.sort (v, ns); | |
2139 | |
2140 octave_idx_type i; | |
2141 if (mode == ASCENDING) | |
2142 { | |
2143 for (i = 0; i < ns; i++) | |
2144 if (sparse_ascending_compare (static_cast<T> (0), v [i])) | |
2145 break; | |
2146 } | |
2147 else | |
2148 { | |
2149 for (i = 0; i < ns; i++) | |
2150 if (sparse_descending_compare (static_cast<T> (0), v [i])) | |
2151 break; | |
2152 } | |
2153 for (octave_idx_type k = 0; k < i; k++) | |
2154 mridx [k] = k; | |
2155 for (octave_idx_type k = i; k < ns; k++) | |
2156 mridx [k] = k - ns + nr; | |
2157 | |
2158 v += ns; | |
2159 mridx += ns; | |
2160 } | |
2161 | |
2162 if (dim > 0) | |
2163 m = m.transpose (); | |
2164 | |
2165 return m; | |
2166 } | |
2167 | |
2168 template <class T> | |
2169 Sparse<T> | |
2170 Sparse<T>::sort (Array<octave_idx_type> &sidx, octave_idx_type dim, | |
2171 sortmode mode) const | |
2172 { | |
2173 Sparse<T> m = *this; | |
2174 | |
2175 octave_idx_type nr = m.rows (); | |
2176 octave_idx_type nc = m.columns (); | |
2177 | |
2178 if (m.length () < 1) | |
2179 { | |
2180 sidx = Array<octave_idx_type> (dim_vector (nr, nc)); | |
2181 return m; | |
2182 } | |
2183 | |
2184 if (dim > 0) | |
2185 { | |
2186 m = m.transpose (); | |
2187 nr = m.rows (); | |
2188 nc = m.columns (); | |
2189 } | |
2190 | |
2191 octave_sort<vec_index<T> *> indexed_sort; | |
2192 | |
2193 if (mode == ASCENDING) | |
2194 indexed_sort.set_compare (sparse_ascending_compare); | |
2195 else if (mode == DESCENDING) | |
2196 indexed_sort.set_compare (sparse_descending_compare); | |
7463
2467639bd8c0
eliminate UNDEFINED sort mode
John W. Eaton <jwe@octave.org>
parents:
7433
diff
changeset
|
2197 else |
2467639bd8c0
eliminate UNDEFINED sort mode
John W. Eaton <jwe@octave.org>
parents:
7433
diff
changeset
|
2198 abort (); |
7433 | 2199 |
2200 T *v = m.data (); | |
2201 octave_idx_type *mcidx = m.cidx (); | |
2202 octave_idx_type *mridx = m.ridx (); | |
2203 | |
2204 OCTAVE_LOCAL_BUFFER (vec_index<T> *, vi, nr); | |
2205 OCTAVE_LOCAL_BUFFER (vec_index<T>, vix, nr); | |
2206 | |
2207 for (octave_idx_type i = 0; i < nr; i++) | |
2208 vi[i] = &vix[i]; | |
2209 | |
2210 sidx = Array<octave_idx_type> (dim_vector (nr, nc)); | |
2211 | |
2212 for (octave_idx_type j = 0; j < nc; j++) | |
2213 { | |
2214 octave_idx_type ns = mcidx [j + 1] - mcidx [j]; | |
2215 octave_idx_type offset = j * nr; | |
2216 | |
2217 if (ns == 0) | |
2218 { | |
2219 for (octave_idx_type k = 0; k < nr; k++) | |
2220 sidx (offset + k) = k; | |
2221 } | |
2222 else | |
2223 { | |
2224 for (octave_idx_type i = 0; i < ns; i++) | |
2225 { | |
2226 vi[i]->vec = v[i]; | |
2227 vi[i]->indx = mridx[i]; | |
2228 } | |
2229 | |
2230 indexed_sort.sort (vi, ns); | |
2231 | |
2232 octave_idx_type i; | |
2233 if (mode == ASCENDING) | |
2234 { | |
2235 for (i = 0; i < ns; i++) | |
2236 if (sparse_ascending_compare (static_cast<T> (0), | |
2237 vi [i] -> vec)) | |
2238 break; | |
2239 } | |
2240 else | |
2241 { | |
2242 for (i = 0; i < ns; i++) | |
2243 if (sparse_descending_compare (static_cast<T> (0), | |
2244 vi [i] -> vec)) | |
2245 break; | |
2246 } | |
2247 | |
2248 octave_idx_type ii = 0; | |
2249 octave_idx_type jj = i; | |
2250 for (octave_idx_type k = 0; k < nr; k++) | |
2251 { | |
2252 if (ii < ns && mridx[ii] == k) | |
2253 ii++; | |
2254 else | |
2255 sidx (offset + jj++) = k; | |
2256 } | |
2257 | |
2258 for (octave_idx_type k = 0; k < i; k++) | |
2259 { | |
2260 v [k] = vi [k] -> vec; | |
2261 sidx (k + offset) = vi [k] -> indx; | |
2262 mridx [k] = k; | |
2263 } | |
2264 | |
2265 for (octave_idx_type k = i; k < ns; k++) | |
2266 { | |
2267 v [k] = vi [k] -> vec; | |
2268 sidx (k - ns + nr + offset) = vi [k] -> indx; | |
2269 mridx [k] = k - ns + nr; | |
2270 } | |
2271 | |
2272 v += ns; | |
2273 mridx += ns; | |
2274 } | |
2275 } | |
2276 | |
2277 if (dim > 0) | |
2278 { | |
2279 m = m.transpose (); | |
2280 sidx = sidx.transpose (); | |
2281 } | |
2282 | |
2283 return m; | |
2284 } | |
2285 | |
7620
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2286 template <class T> |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2287 Sparse<T> |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2288 Sparse<T>::diag (octave_idx_type k) const |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2289 { |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2290 octave_idx_type nnr = rows (); |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2291 octave_idx_type nnc = cols (); |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2292 Sparse<T> d; |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2293 |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2294 if (nnr == 0 || nnc == 0) |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2295 ; // do nothing |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2296 else if (nnr != 1 && nnc != 1) |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2297 { |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2298 if (k > 0) |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2299 nnc -= k; |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2300 else if (k < 0) |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2301 nnr += k; |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2302 |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2303 if (nnr > 0 && nnc > 0) |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2304 { |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2305 octave_idx_type ndiag = (nnr < nnc) ? nnr : nnc; |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2306 |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2307 // Count the number of non-zero elements |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2308 octave_idx_type nel = 0; |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2309 if (k > 0) |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2310 { |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2311 for (octave_idx_type i = 0; i < ndiag; i++) |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2312 if (elem (i, i+k) != 0.) |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2313 nel++; |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2314 } |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2315 else if ( k < 0) |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2316 { |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2317 for (octave_idx_type i = 0; i < ndiag; i++) |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2318 if (elem (i-k, i) != 0.) |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2319 nel++; |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2320 } |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2321 else |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2322 { |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2323 for (octave_idx_type i = 0; i < ndiag; i++) |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2324 if (elem (i, i) != 0.) |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2325 nel++; |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2326 } |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2327 |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2328 d = Sparse<T> (ndiag, 1, nel); |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2329 d.xcidx (0) = 0; |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2330 d.xcidx (1) = nel; |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2331 |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2332 octave_idx_type ii = 0; |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2333 if (k > 0) |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2334 { |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2335 for (octave_idx_type i = 0; i < ndiag; i++) |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2336 { |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2337 T tmp = elem (i, i+k); |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2338 if (tmp != 0.) |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2339 { |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2340 d.xdata (ii) = tmp; |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2341 d.xridx (ii++) = i; |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2342 } |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2343 } |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2344 } |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2345 else if ( k < 0) |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2346 { |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2347 for (octave_idx_type i = 0; i < ndiag; i++) |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2348 { |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2349 T tmp = elem (i-k, i); |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2350 if (tmp != 0.) |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2351 { |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2352 d.xdata (ii) = tmp; |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2353 d.xridx (ii++) = i; |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2354 } |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2355 } |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2356 } |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2357 else |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2358 { |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2359 for (octave_idx_type i = 0; i < ndiag; i++) |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2360 { |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2361 T tmp = elem (i, i); |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2362 if (tmp != 0.) |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2363 { |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2364 d.xdata (ii) = tmp; |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2365 d.xridx (ii++) = i; |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2366 } |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2367 } |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2368 } |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2369 } |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2370 else |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2371 (*current_liboctave_error_handler) |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2372 ("diag: requested diagonal out of range"); |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2373 } |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2374 else if (nnr != 0 && nnc != 0) |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2375 { |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2376 octave_idx_type roff = 0; |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2377 octave_idx_type coff = 0; |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2378 if (k > 0) |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2379 { |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2380 roff = 0; |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2381 coff = k; |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2382 } |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2383 else if (k < 0) |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2384 { |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2385 roff = -k; |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2386 coff = 0; |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2387 } |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2388 |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2389 if (nnr == 1) |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2390 { |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2391 octave_idx_type n = nnc + std::abs (k); |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2392 octave_idx_type nz = nzmax (); |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2393 d = Sparse<T> (n, n, nz); |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2394 for (octave_idx_type i = 0; i < coff+1; i++) |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2395 d.xcidx (i) = 0; |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2396 for (octave_idx_type j = 0; j < nnc; j++) |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2397 { |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2398 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++) |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2399 { |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2400 d.xdata (i) = data (i); |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2401 d.xridx (i) = j + roff; |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2402 } |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2403 d.xcidx (j + coff + 1) = cidx(j+1); |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2404 } |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2405 for (octave_idx_type i = nnc + coff + 1; i < n + 1; i++) |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2406 d.xcidx (i) = nz; |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2407 } |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2408 else |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2409 { |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2410 octave_idx_type n = nnr + std::abs (k); |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2411 octave_idx_type nz = nzmax (); |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2412 octave_idx_type ii = 0; |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2413 octave_idx_type ir = ridx(0); |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2414 d = Sparse<T> (n, n, nz); |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2415 for (octave_idx_type i = 0; i < coff+1; i++) |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2416 d.xcidx (i) = 0; |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2417 for (octave_idx_type i = 0; i < nnr; i++) |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2418 { |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2419 if (ir == i) |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2420 { |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2421 d.xdata (ii) = data (ii); |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2422 d.xridx (ii++) = ir + roff; |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2423 if (ii != nz) |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2424 ir = ridx (ii); |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2425 } |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2426 d.xcidx (i + coff + 1) = ii; |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2427 } |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2428 for (octave_idx_type i = nnr + coff + 1; i < n+1; i++) |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2429 d.xcidx (i) = nz; |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2430 } |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2431 } |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2432 |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2433 return d; |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2434 } |
36594d5bbe13
Move diag function into the octave_value class
David Bateman <dbateman@free.fr>
parents:
7573
diff
changeset
|
2435 |
5775 | 2436 // FIXME |
5164 | 2437 // Unfortunately numel can overflow for very large but very sparse matrices. |
2438 // For now just flag an error when this happens. | |
2439 template <class LT, class RT> | |
2440 int | |
2441 assign1 (Sparse<LT>& lhs, const Sparse<RT>& rhs) | |
2442 { | |
2443 int retval = 1; | |
2444 | |
2445 idx_vector *idx_tmp = lhs.get_idx (); | |
2446 | |
2447 idx_vector lhs_idx = idx_tmp[0]; | |
2448 | |
5275 | 2449 octave_idx_type lhs_len = lhs.numel (); |
2450 octave_idx_type rhs_len = rhs.numel (); | |
5164 | 2451 |
5828 | 2452 uint64_t long_lhs_len = |
2453 static_cast<uint64_t> (lhs.rows ()) * | |
2454 static_cast<uint64_t> (lhs.cols ()); | |
2455 | |
2456 uint64_t long_rhs_len = | |
2457 static_cast<uint64_t> (rhs.rows ()) * | |
2458 static_cast<uint64_t> (rhs.cols ()); | |
2459 | |
2460 if (long_rhs_len != static_cast<uint64_t>(rhs_len) || | |
2461 long_lhs_len != static_cast<uint64_t>(lhs_len)) | |
5164 | 2462 { |
2463 (*current_liboctave_error_handler) | |
2464 ("A(I) = X: Matrix dimensions too large to ensure correct\n", | |
2465 "operation. This is an limitation that should be removed\n", | |
2466 "in the future."); | |
2467 | |
2468 lhs.clear_index (); | |
2469 return 0; | |
2470 } | |
2471 | |
5275 | 2472 octave_idx_type nr = lhs.rows (); |
2473 octave_idx_type nc = lhs.cols (); | |
5681 | 2474 octave_idx_type nz = lhs.nnz (); |
5275 | 2475 |
5781 | 2476 octave_idx_type n = lhs_idx.freeze (lhs_len, "vector", true); |
5164 | 2477 |
2478 if (n != 0) | |
2479 { | |
5275 | 2480 octave_idx_type max_idx = lhs_idx.max () + 1; |
5164 | 2481 max_idx = max_idx < lhs_len ? lhs_len : max_idx; |
2482 | |
2483 // Take a constant copy of lhs. This means that elem won't | |
2484 // create missing elements. | |
2485 const Sparse<LT> c_lhs (lhs); | |
2486 | |
2487 if (rhs_len == n) | |
2488 { | |
5681 | 2489 octave_idx_type new_nzmx = lhs.nnz (); |
5164 | 2490 |
5603 | 2491 OCTAVE_LOCAL_BUFFER (octave_idx_type, rhs_idx, n); |
2492 if (! lhs_idx.is_colon ()) | |
2493 { | |
2494 // Ok here we have to be careful with the indexing, | |
2495 // to treat cases like "a([3,2,1]) = b", and still | |
2496 // handle the need for strict sorting of the sparse | |
2497 // elements. | |
2498 OCTAVE_LOCAL_BUFFER (octave_idx_vector_sort *, sidx, n); | |
2499 OCTAVE_LOCAL_BUFFER (octave_idx_vector_sort, sidxX, n); | |
2500 | |
2501 for (octave_idx_type i = 0; i < n; i++) | |
2502 { | |
2503 sidx[i] = &sidxX[i]; | |
2504 sidx[i]->i = lhs_idx.elem(i); | |
2505 sidx[i]->idx = i; | |
2506 } | |
2507 | |
2508 OCTAVE_QUIT; | |
2509 octave_sort<octave_idx_vector_sort *> | |
2510 sort (octave_idx_vector_comp); | |
2511 | |
2512 sort.sort (sidx, n); | |
2513 | |
2514 intNDArray<octave_idx_type> new_idx (dim_vector (n,1)); | |
2515 | |
2516 for (octave_idx_type i = 0; i < n; i++) | |
2517 { | |
8290
7cbe01c21986
improve dense array indexing
Jaroslav Hajek <highegg@gmail.com>
parents:
8150
diff
changeset
|
2518 new_idx.xelem(i) = sidx[i]->i; |
5603 | 2519 rhs_idx[i] = sidx[i]->idx; |
2520 } | |
2521 | |
2522 lhs_idx = idx_vector (new_idx); | |
2523 } | |
2524 else | |
2525 for (octave_idx_type i = 0; i < n; i++) | |
2526 rhs_idx[i] = i; | |
2527 | |
5164 | 2528 // First count the number of non-zero elements |
5275 | 2529 for (octave_idx_type i = 0; i < n; i++) |
5164 | 2530 { |
2531 OCTAVE_QUIT; | |
2532 | |
5275 | 2533 octave_idx_type ii = lhs_idx.elem (i); |
7546
4249c6fb6e09
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7463
diff
changeset
|
2534 if (i < n - 1 && lhs_idx.elem (i + 1) == ii) |
4249c6fb6e09
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7463
diff
changeset
|
2535 continue; |
5164 | 2536 if (ii < lhs_len && c_lhs.elem(ii) != LT ()) |
5604 | 2537 new_nzmx--; |
5603 | 2538 if (rhs.elem(rhs_idx[i]) != RT ()) |
5604 | 2539 new_nzmx++; |
5164 | 2540 } |
2541 | |
2542 if (nr > 1) | |
2543 { | |
7246 | 2544 Sparse<LT> tmp ((max_idx > nr ? max_idx : nr), 1, new_nzmx); |
5164 | 2545 tmp.cidx(0) = 0; |
5681 | 2546 tmp.cidx(1) = new_nzmx; |
5164 | 2547 |
5275 | 2548 octave_idx_type i = 0; |
2549 octave_idx_type ii = 0; | |
5164 | 2550 if (i < nz) |
2551 ii = c_lhs.ridx(i); | |
2552 | |
5275 | 2553 octave_idx_type j = 0; |
2554 octave_idx_type jj = lhs_idx.elem(j); | |
2555 | |
2556 octave_idx_type kk = 0; | |
5164 | 2557 |
2558 while (j < n || i < nz) | |
2559 { | |
7546
4249c6fb6e09
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7463
diff
changeset
|
2560 if (j < n - 1 && lhs_idx.elem (j + 1) == jj) |
4249c6fb6e09
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7463
diff
changeset
|
2561 { |
4249c6fb6e09
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7463
diff
changeset
|
2562 j++; |
4249c6fb6e09
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7463
diff
changeset
|
2563 jj = lhs_idx.elem (j); |
4249c6fb6e09
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7463
diff
changeset
|
2564 continue; |
4249c6fb6e09
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7463
diff
changeset
|
2565 } |
5164 | 2566 if (j == n || (i < nz && ii < jj)) |
2567 { | |
2568 tmp.xdata (kk) = c_lhs.data (i); | |
2569 tmp.xridx (kk++) = ii; | |
2570 if (++i < nz) | |
2571 ii = c_lhs.ridx(i); | |
2572 } | |
2573 else | |
2574 { | |
5603 | 2575 RT rtmp = rhs.elem (rhs_idx[j]); |
5164 | 2576 if (rtmp != RT ()) |
2577 { | |
2578 tmp.xdata (kk) = rtmp; | |
2579 tmp.xridx (kk++) = jj; | |
2580 } | |
2581 | |
2582 if (ii == jj && i < nz) | |
2583 if (++i < nz) | |
2584 ii = c_lhs.ridx(i); | |
2585 if (++j < n) | |
2586 jj = lhs_idx.elem(j); | |
2587 } | |
2588 } | |
2589 | |
2590 lhs = tmp; | |
2591 } | |
2592 else | |
2593 { | |
7246 | 2594 Sparse<LT> tmp (1, (max_idx > nc ? max_idx : nc), new_nzmx); |
5164 | 2595 |
5275 | 2596 octave_idx_type i = 0; |
2597 octave_idx_type ii = 0; | |
5164 | 2598 while (ii < nc && c_lhs.cidx(ii+1) <= i) |
2599 ii++; | |
2600 | |
5275 | 2601 octave_idx_type j = 0; |
2602 octave_idx_type jj = lhs_idx.elem(j); | |
2603 | |
2604 octave_idx_type kk = 0; | |
2605 octave_idx_type ic = 0; | |
5164 | 2606 |
2607 while (j < n || i < nz) | |
2608 { | |
7546
4249c6fb6e09
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7463
diff
changeset
|
2609 if (j < n - 1 && lhs_idx.elem (j + 1) == jj) |
4249c6fb6e09
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7463
diff
changeset
|
2610 { |
4249c6fb6e09
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7463
diff
changeset
|
2611 j++; |
4249c6fb6e09
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7463
diff
changeset
|
2612 jj = lhs_idx.elem (j); |
4249c6fb6e09
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7463
diff
changeset
|
2613 continue; |
4249c6fb6e09
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7463
diff
changeset
|
2614 } |
5164 | 2615 if (j == n || (i < nz && ii < jj)) |
2616 { | |
2617 while (ic <= ii) | |
2618 tmp.xcidx (ic++) = kk; | |
2619 tmp.xdata (kk) = c_lhs.data (i); | |
5603 | 2620 tmp.xridx (kk++) = 0; |
5164 | 2621 i++; |
2622 while (ii < nc && c_lhs.cidx(ii+1) <= i) | |
2623 ii++; | |
2624 } | |
2625 else | |
2626 { | |
2627 while (ic <= jj) | |
2628 tmp.xcidx (ic++) = kk; | |
2629 | |
5603 | 2630 RT rtmp = rhs.elem (rhs_idx[j]); |
5164 | 2631 if (rtmp != RT ()) |
5603 | 2632 { |
2633 tmp.xdata (kk) = rtmp; | |
2634 tmp.xridx (kk++) = 0; | |
2635 } | |
5164 | 2636 if (ii == jj) |
2637 { | |
2638 i++; | |
2639 while (ii < nc && c_lhs.cidx(ii+1) <= i) | |
2640 ii++; | |
2641 } | |
2642 j++; | |
2643 if (j < n) | |
2644 jj = lhs_idx.elem(j); | |
2645 } | |
2646 } | |
2647 | |
5275 | 2648 for (octave_idx_type iidx = ic; iidx < max_idx+1; iidx++) |
5164 | 2649 tmp.xcidx(iidx) = kk; |
2650 | |
2651 lhs = tmp; | |
2652 } | |
2653 } | |
2654 else if (rhs_len == 1) | |
2655 { | |
5681 | 2656 octave_idx_type new_nzmx = lhs.nnz (); |
5164 | 2657 RT scalar = rhs.elem (0); |
2658 bool scalar_non_zero = (scalar != RT ()); | |
5603 | 2659 lhs_idx.sort (true); |
7546
4249c6fb6e09
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7463
diff
changeset
|
2660 n = lhs_idx.length (n); |
5164 | 2661 |
2662 // First count the number of non-zero elements | |
2663 if (scalar != RT ()) | |
5604 | 2664 new_nzmx += n; |
5275 | 2665 for (octave_idx_type i = 0; i < n; i++) |
5164 | 2666 { |
2667 OCTAVE_QUIT; | |
2668 | |
5275 | 2669 octave_idx_type ii = lhs_idx.elem (i); |
5164 | 2670 if (ii < lhs_len && c_lhs.elem(ii) != LT ()) |
5604 | 2671 new_nzmx--; |
5164 | 2672 } |
2673 | |
2674 if (nr > 1) | |
2675 { | |
7246 | 2676 Sparse<LT> tmp ((max_idx > nr ? max_idx : nr), 1, new_nzmx); |
5164 | 2677 tmp.cidx(0) = 0; |
5681 | 2678 tmp.cidx(1) = new_nzmx; |
5164 | 2679 |
5275 | 2680 octave_idx_type i = 0; |
2681 octave_idx_type ii = 0; | |
5164 | 2682 if (i < nz) |
2683 ii = c_lhs.ridx(i); | |
2684 | |
5275 | 2685 octave_idx_type j = 0; |
2686 octave_idx_type jj = lhs_idx.elem(j); | |
2687 | |
2688 octave_idx_type kk = 0; | |
5164 | 2689 |
2690 while (j < n || i < nz) | |
2691 { | |
2692 if (j == n || (i < nz && ii < jj)) | |
2693 { | |
2694 tmp.xdata (kk) = c_lhs.data (i); | |
2695 tmp.xridx (kk++) = ii; | |
2696 if (++i < nz) | |
2697 ii = c_lhs.ridx(i); | |
2698 } | |
2699 else | |
2700 { | |
2701 if (scalar_non_zero) | |
2702 { | |
2703 tmp.xdata (kk) = scalar; | |
2704 tmp.xridx (kk++) = jj; | |
2705 } | |
2706 | |
2707 if (ii == jj && i < nz) | |
2708 if (++i < nz) | |
2709 ii = c_lhs.ridx(i); | |
2710 if (++j < n) | |
2711 jj = lhs_idx.elem(j); | |
2712 } | |
2713 } | |
2714 | |
2715 lhs = tmp; | |
2716 } | |
2717 else | |
2718 { | |
7246 | 2719 Sparse<LT> tmp (1, (max_idx > nc ? max_idx : nc), new_nzmx); |
5164 | 2720 |
5275 | 2721 octave_idx_type i = 0; |
2722 octave_idx_type ii = 0; | |
5164 | 2723 while (ii < nc && c_lhs.cidx(ii+1) <= i) |
2724 ii++; | |
2725 | |
5275 | 2726 octave_idx_type j = 0; |
2727 octave_idx_type jj = lhs_idx.elem(j); | |
2728 | |
2729 octave_idx_type kk = 0; | |
2730 octave_idx_type ic = 0; | |
5164 | 2731 |
2732 while (j < n || i < nz) | |
2733 { | |
2734 if (j == n || (i < nz && ii < jj)) | |
2735 { | |
2736 while (ic <= ii) | |
2737 tmp.xcidx (ic++) = kk; | |
2738 tmp.xdata (kk) = c_lhs.data (i); | |
2739 i++; | |
2740 while (ii < nc && c_lhs.cidx(ii+1) <= i) | |
2741 ii++; | |
2742 } | |
2743 else | |
2744 { | |
2745 while (ic <= jj) | |
2746 tmp.xcidx (ic++) = kk; | |
2747 if (scalar_non_zero) | |
2748 tmp.xdata (kk) = scalar; | |
2749 if (ii == jj) | |
2750 { | |
2751 i++; | |
2752 while (ii < nc && c_lhs.cidx(ii+1) <= i) | |
2753 ii++; | |
2754 } | |
2755 j++; | |
2756 if (j < n) | |
2757 jj = lhs_idx.elem(j); | |
2758 } | |
2759 tmp.xridx (kk++) = 0; | |
2760 } | |
2761 | |
5275 | 2762 for (octave_idx_type iidx = ic; iidx < max_idx+1; iidx++) |
5164 | 2763 tmp.xcidx(iidx) = kk; |
2764 | |
2765 lhs = tmp; | |
2766 } | |
2767 } | |
2768 else | |
2769 { | |
2770 (*current_liboctave_error_handler) | |
2771 ("A(I) = X: X must be a scalar or a vector with same length as I"); | |
2772 | |
2773 retval = 0; | |
2774 } | |
2775 } | |
2776 else if (lhs_idx.is_colon ()) | |
2777 { | |
2778 if (lhs_len == 0) | |
2779 { | |
2780 | |
5681 | 2781 octave_idx_type new_nzmx = rhs.nnz (); |
5604 | 2782 Sparse<LT> tmp (1, rhs_len, new_nzmx); |
5164 | 2783 |
5275 | 2784 octave_idx_type ii = 0; |
2785 octave_idx_type jj = 0; | |
2786 for (octave_idx_type i = 0; i < rhs.cols(); i++) | |
2787 for (octave_idx_type j = rhs.cidx(i); j < rhs.cidx(i+1); j++) | |
5164 | 2788 { |
2789 OCTAVE_QUIT; | |
5275 | 2790 for (octave_idx_type k = jj; k <= i * rhs.rows() + rhs.ridx(j); k++) |
5164 | 2791 tmp.cidx(jj++) = ii; |
2792 | |
2793 tmp.data(ii) = rhs.data(j); | |
2794 tmp.ridx(ii++) = 0; | |
2795 } | |
2796 | |
5275 | 2797 for (octave_idx_type i = jj; i < rhs_len + 1; i++) |
5164 | 2798 tmp.cidx(i) = ii; |
2799 | |
2800 lhs = tmp; | |
2801 } | |
2802 else | |
2803 (*current_liboctave_error_handler) | |
2804 ("A(:) = X: A must be the same size as X"); | |
2805 } | |
2806 else if (! (rhs_len == 1 || rhs_len == 0)) | |
2807 { | |
2808 (*current_liboctave_error_handler) | |
2809 ("A([]) = X: X must also be an empty matrix or a scalar"); | |
2810 | |
2811 retval = 0; | |
2812 } | |
2813 | |
2814 lhs.clear_index (); | |
2815 | |
2816 return retval; | |
2817 } | |
2818 | |
2819 template <class LT, class RT> | |
2820 int | |
2821 assign (Sparse<LT>& lhs, const Sparse<RT>& rhs) | |
2822 { | |
2823 int retval = 1; | |
2824 | |
2825 int n_idx = lhs.index_count (); | |
2826 | |
5275 | 2827 octave_idx_type lhs_nr = lhs.rows (); |
2828 octave_idx_type lhs_nc = lhs.cols (); | |
5681 | 2829 octave_idx_type lhs_nz = lhs.nnz (); |
5275 | 2830 |
2831 octave_idx_type rhs_nr = rhs.rows (); | |
2832 octave_idx_type rhs_nc = rhs.cols (); | |
5164 | 2833 |
2834 idx_vector *tmp = lhs.get_idx (); | |
2835 | |
2836 idx_vector idx_i; | |
2837 idx_vector idx_j; | |
2838 | |
2839 if (n_idx > 2) | |
2840 { | |
2841 (*current_liboctave_error_handler) | |
2842 ("A(I, J) = X: can only have 1 or 2 indexes for sparse matrices"); | |
6092 | 2843 |
2844 lhs.clear_index (); | |
5164 | 2845 return 0; |
2846 } | |
2847 | |
2848 if (n_idx > 1) | |
2849 idx_j = tmp[1]; | |
2850 | |
2851 if (n_idx > 0) | |
2852 idx_i = tmp[0]; | |
2853 | |
6988 | 2854 // Take a constant copy of lhs. This means that ridx and family won't |
2855 // call make_unique. | |
2856 const Sparse<LT> c_lhs (lhs); | |
2857 | |
5164 | 2858 if (n_idx == 2) |
2859 { | |
5781 | 2860 octave_idx_type n = idx_i.freeze (lhs_nr, "row", true); |
2861 octave_idx_type m = idx_j.freeze (lhs_nc, "column", true); | |
5164 | 2862 |
2863 int idx_i_is_colon = idx_i.is_colon (); | |
2864 int idx_j_is_colon = idx_j.is_colon (); | |
2865 | |
7238 | 2866 if (lhs_nr == 0 && lhs_nc == 0) |
2867 { | |
2868 if (idx_i_is_colon) | |
2869 n = rhs_nr; | |
2870 | |
2871 if (idx_j_is_colon) | |
2872 m = rhs_nc; | |
2873 } | |
5164 | 2874 |
2875 if (idx_i && idx_j) | |
8150
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2876 { |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2877 if (rhs_nr == 1 && rhs_nc == 1 && n >= 0 && m >= 0) |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2878 { |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2879 if (n > 0 && m > 0) |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2880 { |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2881 idx_i.sort (true); |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2882 n = idx_i.length (n); |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2883 idx_j.sort (true); |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2884 m = idx_j.length (m); |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2885 |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2886 octave_idx_type max_row_idx = idx_i_is_colon ? rhs_nr : |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2887 idx_i.max () + 1; |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2888 octave_idx_type max_col_idx = idx_j_is_colon ? rhs_nc : |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2889 idx_j.max () + 1; |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2890 octave_idx_type new_nr = max_row_idx > lhs_nr ? |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2891 max_row_idx : lhs_nr; |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2892 octave_idx_type new_nc = max_col_idx > lhs_nc ? |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2893 max_col_idx : lhs_nc; |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2894 RT scalar = rhs.elem (0, 0); |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2895 |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2896 // Count the number of non-zero terms |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2897 octave_idx_type new_nzmx = lhs.nnz (); |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2898 for (octave_idx_type j = 0; j < m; j++) |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2899 { |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2900 octave_idx_type jj = idx_j.elem (j); |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2901 if (jj < lhs_nc) |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2902 { |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2903 for (octave_idx_type i = 0; i < n; i++) |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2904 { |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2905 OCTAVE_QUIT; |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2906 |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2907 octave_idx_type ii = idx_i.elem (i); |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2908 |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2909 if (ii < lhs_nr) |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2910 { |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2911 for (octave_idx_type k = c_lhs.cidx(jj); |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2912 k < c_lhs.cidx(jj+1); k++) |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2913 { |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2914 if (c_lhs.ridx(k) == ii) |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2915 new_nzmx--; |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2916 if (c_lhs.ridx(k) >= ii) |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2917 break; |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2918 } |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2919 } |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2920 } |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2921 } |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2922 } |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2923 |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2924 if (scalar != RT()) |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2925 new_nzmx += m * n; |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2926 |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2927 Sparse<LT> stmp (new_nr, new_nc, new_nzmx); |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2928 |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2929 octave_idx_type jji = 0; |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2930 octave_idx_type jj = idx_j.elem (jji); |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2931 octave_idx_type kk = 0; |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2932 stmp.cidx(0) = 0; |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2933 for (octave_idx_type j = 0; j < new_nc; j++) |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2934 { |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2935 if (jji < m && jj == j) |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2936 { |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2937 octave_idx_type iii = 0; |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2938 octave_idx_type ii = idx_i.elem (iii); |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2939 octave_idx_type ppp = 0; |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2940 octave_idx_type ppi = (j >= lhs_nc ? 0 : |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2941 c_lhs.cidx(j+1) - |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2942 c_lhs.cidx(j)); |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2943 octave_idx_type pp = (ppp < ppi ? |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2944 c_lhs.ridx(c_lhs.cidx(j)+ppp) : |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2945 new_nr); |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2946 while (ppp < ppi || iii < n) |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2947 { |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2948 if (iii < n && ii <= pp) |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2949 { |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2950 if (scalar != RT ()) |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2951 { |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2952 stmp.data(kk) = scalar; |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2953 stmp.ridx(kk++) = ii; |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2954 } |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2955 if (ii == pp) |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2956 pp = (++ppp < ppi ? c_lhs.ridx(c_lhs.cidx(j)+ppp) : new_nr); |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2957 if (++iii < n) |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2958 ii = idx_i.elem(iii); |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2959 } |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2960 else |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2961 { |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2962 stmp.data(kk) = |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2963 c_lhs.data(c_lhs.cidx(j)+ppp); |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2964 stmp.ridx(kk++) = pp; |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2965 pp = (++ppp < ppi ? c_lhs.ridx(c_lhs.cidx(j)+ppp) : new_nr); |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2966 } |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2967 } |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2968 if (++jji < m) |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2969 jj = idx_j.elem(jji); |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2970 } |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2971 else if (j < lhs_nc) |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2972 { |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2973 for (octave_idx_type i = c_lhs.cidx(j); |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2974 i < c_lhs.cidx(j+1); i++) |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2975 { |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2976 stmp.data(kk) = c_lhs.data(i); |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2977 stmp.ridx(kk++) = c_lhs.ridx(i); |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2978 } |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2979 } |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2980 stmp.cidx(j+1) = kk; |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2981 } |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2982 |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2983 lhs = stmp; |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2984 } |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2985 else |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2986 { |
7253 | 2987 #if 0 |
8150
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2988 // FIXME -- the following code will make this |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2989 // function behave the same as the full matrix |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2990 // case for things like |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2991 // |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2992 // x = sparse (ones (2)); |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2993 // x([],3) = 2; |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2994 // |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2995 // x = |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2996 // |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2997 // Compressed Column Sparse (rows = 2, cols = 3, nnz = 4) |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2998 // |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
2999 // (1, 1) -> 1 |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3000 // (2, 1) -> 1 |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3001 // (1, 2) -> 1 |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3002 // (2, 2) -> 1 |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3003 // |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3004 // However, Matlab doesn't resize in this case |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3005 // even though it does in the full matrix case. |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3006 |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3007 if (n > 0) |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3008 { |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3009 octave_idx_type max_row_idx = idx_i_is_colon ? |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3010 rhs_nr : idx_i.max () + 1; |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3011 octave_idx_type new_nr = max_row_idx > lhs_nr ? |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3012 max_row_idx : lhs_nr; |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3013 octave_idx_type new_nc = lhs_nc; |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3014 |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3015 lhs.resize (new_nr, new_nc); |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3016 } |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3017 else if (m > 0) |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3018 { |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3019 octave_idx_type max_col_idx = idx_j_is_colon ? |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3020 rhs_nc : idx_j.max () + 1; |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3021 octave_idx_type new_nr = lhs_nr; |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3022 octave_idx_type new_nc = max_col_idx > lhs_nc ? |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3023 max_col_idx : lhs_nc; |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3024 |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3025 lhs.resize (new_nr, new_nc); |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3026 } |
7253 | 3027 #endif |
8150
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3028 } |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3029 } |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3030 else if (n == rhs_nr && m == rhs_nc) |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3031 { |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3032 if (n > 0 && m > 0) |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3033 { |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3034 octave_idx_type max_row_idx = idx_i_is_colon ? rhs_nr : |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3035 idx_i.max () + 1; |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3036 octave_idx_type max_col_idx = idx_j_is_colon ? rhs_nc : |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3037 idx_j.max () + 1; |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3038 octave_idx_type new_nr = max_row_idx > lhs_nr ? |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3039 max_row_idx : lhs_nr; |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3040 octave_idx_type new_nc = max_col_idx > lhs_nc ? |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3041 max_col_idx : lhs_nc; |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3042 |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3043 OCTAVE_LOCAL_BUFFER (octave_idx_type, rhs_idx_i, n); |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3044 if (! idx_i.is_colon ()) |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3045 { |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3046 // Ok here we have to be careful with the indexing, |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3047 // to treat cases like "a([3,2,1],:) = b", and still |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3048 // handle the need for strict sorting of the sparse |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3049 // elements. |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3050 OCTAVE_LOCAL_BUFFER (octave_idx_vector_sort *, |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3051 sidx, n); |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3052 OCTAVE_LOCAL_BUFFER (octave_idx_vector_sort, |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3053 sidxX, n); |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3054 |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3055 for (octave_idx_type i = 0; i < n; i++) |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3056 { |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3057 sidx[i] = &sidxX[i]; |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3058 sidx[i]->i = idx_i.elem(i); |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3059 sidx[i]->idx = i; |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3060 } |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3061 |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3062 OCTAVE_QUIT; |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3063 octave_sort<octave_idx_vector_sort *> |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3064 sort (octave_idx_vector_comp); |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3065 |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3066 sort.sort (sidx, n); |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3067 |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3068 intNDArray<octave_idx_type> new_idx (dim_vector (n,1)); |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3069 |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3070 for (octave_idx_type i = 0; i < n; i++) |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3071 { |
8290
7cbe01c21986
improve dense array indexing
Jaroslav Hajek <highegg@gmail.com>
parents:
8150
diff
changeset
|
3072 new_idx.xelem(i) = sidx[i]->i; |
8150
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3073 rhs_idx_i[i] = sidx[i]->idx; |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3074 } |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3075 |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3076 idx_i = idx_vector (new_idx); |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3077 } |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3078 else |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3079 for (octave_idx_type i = 0; i < n; i++) |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3080 rhs_idx_i[i] = i; |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3081 |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3082 OCTAVE_LOCAL_BUFFER (octave_idx_type, rhs_idx_j, m); |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3083 if (! idx_j.is_colon ()) |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3084 { |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3085 // Ok here we have to be careful with the indexing, |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3086 // to treat cases like "a([3,2,1],:) = b", and still |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3087 // handle the need for strict sorting of the sparse |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3088 // elements. |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3089 OCTAVE_LOCAL_BUFFER (octave_idx_vector_sort *, |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3090 sidx, m); |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3091 OCTAVE_LOCAL_BUFFER (octave_idx_vector_sort, |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3092 sidxX, m); |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3093 |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3094 for (octave_idx_type i = 0; i < m; i++) |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3095 { |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3096 sidx[i] = &sidxX[i]; |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3097 sidx[i]->i = idx_j.elem(i); |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3098 sidx[i]->idx = i; |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3099 } |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3100 |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3101 OCTAVE_QUIT; |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3102 octave_sort<octave_idx_vector_sort *> |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3103 sort (octave_idx_vector_comp); |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3104 |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3105 sort.sort (sidx, m); |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3106 |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3107 intNDArray<octave_idx_type> new_idx (dim_vector (m,1)); |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3108 |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3109 for (octave_idx_type i = 0; i < m; i++) |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3110 { |
8290
7cbe01c21986
improve dense array indexing
Jaroslav Hajek <highegg@gmail.com>
parents:
8150
diff
changeset
|
3111 new_idx.xelem(i) = sidx[i]->i; |
8150
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3112 rhs_idx_j[i] = sidx[i]->idx; |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3113 } |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3114 |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3115 idx_j = idx_vector (new_idx); |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3116 } |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3117 else |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3118 for (octave_idx_type i = 0; i < m; i++) |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3119 rhs_idx_j[i] = i; |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3120 |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3121 // Maximum number of non-zero elements |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3122 octave_idx_type new_nzmx = lhs.nnz() + rhs.nnz(); |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3123 |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3124 Sparse<LT> stmp (new_nr, new_nc, new_nzmx); |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3125 |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3126 octave_idx_type jji = 0; |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3127 octave_idx_type jj = idx_j.elem (jji); |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3128 octave_idx_type kk = 0; |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3129 stmp.cidx(0) = 0; |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3130 for (octave_idx_type j = 0; j < new_nc; j++) |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3131 { |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3132 if (jji < m && jj == j) |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3133 { |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3134 octave_idx_type iii = 0; |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3135 octave_idx_type ii = idx_i.elem (iii); |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3136 octave_idx_type ppp = 0; |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3137 octave_idx_type ppi = (j >= lhs_nc ? 0 : |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3138 c_lhs.cidx(j+1) - |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3139 c_lhs.cidx(j)); |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3140 octave_idx_type pp = (ppp < ppi ? |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3141 c_lhs.ridx(c_lhs.cidx(j)+ppp) : |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3142 new_nr); |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3143 while (ppp < ppi || iii < n) |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3144 { |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3145 if (iii < n && ii <= pp) |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3146 { |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3147 if (iii < n - 1 && |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3148 idx_i.elem (iii + 1) == ii) |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3149 { |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3150 iii++; |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3151 ii = idx_i.elem(iii); |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3152 continue; |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3153 } |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3154 |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3155 RT rtmp = rhs.elem (rhs_idx_i[iii], |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3156 rhs_idx_j[jji]); |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3157 if (rtmp != RT ()) |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3158 { |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3159 stmp.data(kk) = rtmp; |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3160 stmp.ridx(kk++) = ii; |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3161 } |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3162 if (ii == pp) |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3163 pp = (++ppp < ppi ? c_lhs.ridx(c_lhs.cidx(j)+ppp) : new_nr); |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3164 if (++iii < n) |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3165 ii = idx_i.elem(iii); |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3166 } |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3167 else |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3168 { |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3169 stmp.data(kk) = |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3170 c_lhs.data(c_lhs.cidx(j)+ppp); |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3171 stmp.ridx(kk++) = pp; |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3172 pp = (++ppp < ppi ? c_lhs.ridx(c_lhs.cidx(j)+ppp) : new_nr); |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3173 } |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3174 } |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3175 if (++jji < m) |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3176 jj = idx_j.elem(jji); |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3177 } |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3178 else if (j < lhs_nc) |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3179 { |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3180 for (octave_idx_type i = c_lhs.cidx(j); |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3181 i < c_lhs.cidx(j+1); i++) |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3182 { |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3183 stmp.data(kk) = c_lhs.data(i); |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3184 stmp.ridx(kk++) = c_lhs.ridx(i); |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3185 } |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3186 } |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3187 stmp.cidx(j+1) = kk; |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3188 } |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3189 |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3190 stmp.maybe_compress(); |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3191 lhs = stmp; |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3192 } |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3193 } |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3194 else if (n == 0 && m == 0) |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3195 { |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3196 if (! ((rhs_nr == 1 && rhs_nc == 1) |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3197 || (rhs_nr == 0 || rhs_nc == 0))) |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3198 { |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3199 (*current_liboctave_error_handler) |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3200 ("A([], []) = X: X must be an empty matrix or a scalar"); |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3201 |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3202 retval = 0; |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3203 } |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3204 } |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3205 else |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3206 { |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3207 (*current_liboctave_error_handler) |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3208 ("A(I, J) = X: X must be a scalar or the number of elements in I must"); |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3209 (*current_liboctave_error_handler) |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3210 ("match the number of rows in X and the number of elements in J must"); |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3211 (*current_liboctave_error_handler) |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3212 ("match the number of columns in X"); |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3213 |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3214 retval = 0; |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3215 } |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3216 } |
5164 | 3217 // idx_vector::freeze() printed an error message for us. |
3218 } | |
3219 else if (n_idx == 1) | |
3220 { | |
3221 int lhs_is_empty = lhs_nr == 0 || lhs_nc == 0; | |
3222 | |
3223 if (lhs_is_empty || (lhs_nr == 1 && lhs_nc == 1)) | |
3224 { | |
5275 | 3225 octave_idx_type lhs_len = lhs.length (); |
3226 | |
8677
095ae5e0a831
eliminte some compiler warnings
John W. Eaton <jwe@octave.org>
parents:
8527
diff
changeset
|
3227 // Called for side-effects on idx_i. |
095ae5e0a831
eliminte some compiler warnings
John W. Eaton <jwe@octave.org>
parents:
8527
diff
changeset
|
3228 idx_i.freeze (lhs_len, 0, true); |
5164 | 3229 |
3230 if (idx_i) | |
8150
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3231 { |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3232 if (lhs_is_empty |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3233 && idx_i.is_colon () |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3234 && ! (rhs_nr == 1 || rhs_nc == 1)) |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3235 { |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3236 (*current_liboctave_warning_with_id_handler) |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3237 ("Octave:fortran-indexing", |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3238 "A(:) = X: X is not a vector or scalar"); |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3239 } |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3240 else |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3241 { |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3242 octave_idx_type idx_nr = idx_i.orig_rows (); |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3243 octave_idx_type idx_nc = idx_i.orig_columns (); |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3244 |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3245 if (! (rhs_nr == idx_nr && rhs_nc == idx_nc)) |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3246 (*current_liboctave_warning_with_id_handler) |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3247 ("Octave:fortran-indexing", |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3248 "A(I) = X: X does not have same shape as I"); |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3249 } |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3250 |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3251 if (! assign1 (lhs, rhs)) |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3252 retval = 0; |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3253 } |
5164 | 3254 // idx_vector::freeze() printed an error message for us. |
3255 } | |
3256 else if (lhs_nr == 1) | |
3257 { | |
5781 | 3258 idx_i.freeze (lhs_nc, "vector", true); |
5164 | 3259 |
3260 if (idx_i) | |
8150
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3261 { |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3262 if (! assign1 (lhs, rhs)) |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3263 retval = 0; |
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3264 } |
5164 | 3265 // idx_vector::freeze() printed an error message for us. |
3266 } | |
3267 else if (lhs_nc == 1) | |
3268 { | |
5781 | 3269 idx_i.freeze (lhs_nr, "vector", true); |
5164 | 3270 |
3271 if (idx_i) | |
3272 { | |
8150
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3273 if (! assign1 (lhs, rhs)) |
5164 | 3274 retval = 0; |
3275 } | |
3276 // idx_vector::freeze() printed an error message for us. | |
3277 } | |
3278 else | |
3279 { | |
7573
755bf7ecc29b
eliminate one_zero stuff from idx_vector
John W. Eaton <jwe@octave.org>
parents:
7546
diff
changeset
|
3280 if (! idx_i.is_colon ()) |
5781 | 3281 (*current_liboctave_warning_with_id_handler) |
3282 ("Octave:fortran-indexing", "single index used for matrix"); | |
5164 | 3283 |
5275 | 3284 octave_idx_type lhs_len = lhs.length (); |
3285 | |
3286 octave_idx_type len = idx_i.freeze (lhs_nr * lhs_nc, "matrix"); | |
5164 | 3287 |
3288 if (idx_i) | |
3289 { | |
8150
283989f2da9b
make null assignment matlab compatible
Jaroslav Hajek <highegg@gmail.com>
parents:
7717
diff
changeset
|
3290 if (len == 0) |
5164 | 3291 { |
3292 if (! ((rhs_nr == 1 && rhs_nc == 1) | |
3293 || (rhs_nr == 0 || rhs_nc == 0))) | |
3294 (*current_liboctave_error_handler) | |
3295 ("A([]) = X: X must be an empty matrix or scalar"); | |
3296 } | |
3297 else if (len == rhs_nr * rhs_nc) | |
3298 { | |
5604 | 3299 octave_idx_type new_nzmx = lhs_nz; |
5603 | 3300 OCTAVE_LOCAL_BUFFER (octave_idx_type, rhs_idx, len); |
3301 | |
3302 if (! idx_i.is_colon ()) | |
3303 { | |
3304 // Ok here we have to be careful with the indexing, to | |
3305 // treat cases like "a([3,2,1]) = b", and still handle | |
3306 // the need for strict sorting of the sparse elements. | |
3307 | |
3308 OCTAVE_LOCAL_BUFFER (octave_idx_vector_sort *, sidx, | |
3309 len); | |
3310 OCTAVE_LOCAL_BUFFER (octave_idx_vector_sort, sidxX, | |
3311 len); | |
3312 | |
3313 for (octave_idx_type i = 0; i < len; i++) | |
3314 { | |
3315 sidx[i] = &sidxX[i]; | |
3316 sidx[i]->i = idx_i.elem(i); | |
3317 sidx[i]->idx = i; | |
3318 } | |
3319 | |
3320 OCTAVE_QUIT; | |
3321 octave_sort<octave_idx_vector_sort *> | |
3322 sort (octave_idx_vector_comp); | |
3323 | |
3324 sort.sort (sidx, len); | |
3325 | |
3326 intNDArray<octave_idx_type> new_idx (dim_vector (len,1)); | |
3327 | |
3328 for (octave_idx_type i = 0; i < len; i++) | |
3329 { | |
8290
7cbe01c21986
improve dense array indexing
Jaroslav Hajek <highegg@gmail.com>
parents:
8150
diff
changeset
|
3330 new_idx.xelem(i) = sidx[i]->i; |
5603 | 3331 rhs_idx[i] = sidx[i]->idx; |
3332 } | |
3333 | |
3334 idx_i = idx_vector (new_idx); | |
3335 } | |
3336 else | |
3337 for (octave_idx_type i = 0; i < len; i++) | |
3338 rhs_idx[i] = i; | |
5164 | 3339 |
3340 // First count the number of non-zero elements | |
5275 | 3341 for (octave_idx_type i = 0; i < len; i++) |
5164 | 3342 { |
3343 OCTAVE_QUIT; | |
3344 | |
5275 | 3345 octave_idx_type ii = idx_i.elem (i); |
7546
4249c6fb6e09
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7463
diff
changeset
|
3346 if (i < len - 1 && idx_i.elem (i + 1) == ii) |
4249c6fb6e09
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7463
diff
changeset
|
3347 continue; |
5164 | 3348 if (ii < lhs_len && c_lhs.elem(ii) != LT ()) |
5604 | 3349 new_nzmx--; |
5603 | 3350 if (rhs.elem(rhs_idx[i]) != RT ()) |
5604 | 3351 new_nzmx++; |
5164 | 3352 } |
3353 | |
5604 | 3354 Sparse<LT> stmp (lhs_nr, lhs_nc, new_nzmx); |
5164 | 3355 |
5275 | 3356 octave_idx_type i = 0; |
3357 octave_idx_type ii = 0; | |
3358 octave_idx_type ic = 0; | |
5164 | 3359 if (i < lhs_nz) |
3360 { | |
3361 while (ic < lhs_nc && i >= c_lhs.cidx(ic+1)) | |
3362 ic++; | |
3363 ii = ic * lhs_nr + c_lhs.ridx(i); | |
3364 } | |
3365 | |
5275 | 3366 octave_idx_type j = 0; |
3367 octave_idx_type jj = idx_i.elem (j); | |
3368 octave_idx_type jr = jj % lhs_nr; | |
3369 octave_idx_type jc = (jj - jr) / lhs_nr; | |
3370 | |
3371 octave_idx_type kk = 0; | |
3372 octave_idx_type kc = 0; | |
5164 | 3373 |
3374 while (j < len || i < lhs_nz) | |
3375 { | |
7546
4249c6fb6e09
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7463
diff
changeset
|
3376 if (j < len - 1 && idx_i.elem (j + 1) == jj) |
4249c6fb6e09
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7463
diff
changeset
|
3377 { |
4249c6fb6e09
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7463
diff
changeset
|
3378 j++; |
4249c6fb6e09
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7463
diff
changeset
|
3379 jj = idx_i.elem (j); |
4249c6fb6e09
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7463
diff
changeset
|
3380 jr = jj % lhs_nr; |
4249c6fb6e09
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7463
diff
changeset
|
3381 jc = (jj - jr) / lhs_nr; |
4249c6fb6e09
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7463
diff
changeset
|
3382 continue; |
4249c6fb6e09
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7463
diff
changeset
|
3383 } |
4249c6fb6e09
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7463
diff
changeset
|
3384 |
5164 | 3385 if (j == len || (i < lhs_nz && ii < jj)) |
3386 { | |
3387 while (kc <= ic) | |
3388 stmp.xcidx (kc++) = kk; | |
3389 stmp.xdata (kk) = c_lhs.data (i); | |
3390 stmp.xridx (kk++) = c_lhs.ridx (i); | |
3391 i++; | |
3392 while (ic < lhs_nc && i >= c_lhs.cidx(ic+1)) | |
3393 ic++; | |
3394 if (i < lhs_nz) | |
3395 ii = ic * lhs_nr + c_lhs.ridx(i); | |
3396 } | |
3397 else | |
3398 { | |
3399 while (kc <= jc) | |
3400 stmp.xcidx (kc++) = kk; | |
5603 | 3401 RT rtmp = rhs.elem (rhs_idx[j]); |
5164 | 3402 if (rtmp != RT ()) |
3403 { | |
3404 stmp.xdata (kk) = rtmp; | |
3405 stmp.xridx (kk++) = jr; | |
3406 } | |
3407 if (ii == jj) | |
3408 { | |
3409 i++; | |
3410 while (ic < lhs_nc && i >= c_lhs.cidx(ic+1)) | |
3411 ic++; | |
3412 if (i < lhs_nz) | |
3413 ii = ic * lhs_nr + c_lhs.ridx(i); | |
3414 } | |
3415 j++; | |
3416 if (j < len) | |
3417 { | |
3418 jj = idx_i.elem (j); | |
3419 jr = jj % lhs_nr; | |
3420 jc = (jj - jr) / lhs_nr; | |
3421 } | |
3422 } | |
3423 } | |
3424 | |
5275 | 3425 for (octave_idx_type iidx = kc; iidx < lhs_nc+1; iidx++) |
5603 | 3426 stmp.xcidx(iidx) = kk; |
5164 | 3427 |
3428 lhs = stmp; | |
3429 } | |
3430 else if (rhs_nr == 1 && rhs_nc == 1) | |
3431 { | |
3432 RT scalar = rhs.elem (0, 0); | |
5604 | 3433 octave_idx_type new_nzmx = lhs_nz; |
5603 | 3434 idx_i.sort (true); |
7546
4249c6fb6e09
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7463
diff
changeset
|
3435 len = idx_i.length (len); |
5164 | 3436 |
3437 // First count the number of non-zero elements | |
3438 if (scalar != RT ()) | |
5604 | 3439 new_nzmx += len; |
5275 | 3440 for (octave_idx_type i = 0; i < len; i++) |
5164 | 3441 { |
3442 OCTAVE_QUIT; | |
5275 | 3443 octave_idx_type ii = idx_i.elem (i); |
5164 | 3444 if (ii < lhs_len && c_lhs.elem(ii) != LT ()) |
5604 | 3445 new_nzmx--; |
5164 | 3446 } |
3447 | |
5604 | 3448 Sparse<LT> stmp (lhs_nr, lhs_nc, new_nzmx); |
5164 | 3449 |
5275 | 3450 octave_idx_type i = 0; |
3451 octave_idx_type ii = 0; | |
3452 octave_idx_type ic = 0; | |
5164 | 3453 if (i < lhs_nz) |
3454 { | |
3455 while (ic < lhs_nc && i >= c_lhs.cidx(ic+1)) | |
3456 ic++; | |
3457 ii = ic * lhs_nr + c_lhs.ridx(i); | |
3458 } | |
3459 | |
5275 | 3460 octave_idx_type j = 0; |
3461 octave_idx_type jj = idx_i.elem (j); | |
3462 octave_idx_type jr = jj % lhs_nr; | |
3463 octave_idx_type jc = (jj - jr) / lhs_nr; | |
3464 | |
3465 octave_idx_type kk = 0; | |
3466 octave_idx_type kc = 0; | |
5164 | 3467 |
3468 while (j < len || i < lhs_nz) | |
3469 { | |
3470 if (j == len || (i < lhs_nz && ii < jj)) | |
3471 { | |
3472 while (kc <= ic) | |
3473 stmp.xcidx (kc++) = kk; | |
3474 stmp.xdata (kk) = c_lhs.data (i); | |
3475 stmp.xridx (kk++) = c_lhs.ridx (i); | |
3476 i++; | |
3477 while (ic < lhs_nc && i >= c_lhs.cidx(ic+1)) | |
3478 ic++; | |
3479 if (i < lhs_nz) | |
3480 ii = ic * lhs_nr + c_lhs.ridx(i); | |
3481 } | |
3482 else | |
3483 { | |
3484 while (kc <= jc) | |
3485 stmp.xcidx (kc++) = kk; | |
3486 if (scalar != RT ()) | |
3487 { | |
3488 stmp.xdata (kk) = scalar; | |
3489 stmp.xridx (kk++) = jr; | |
3490 } | |
3491 if (ii == jj) | |
3492 { | |
3493 i++; | |
3494 while (ic < lhs_nc && i >= c_lhs.cidx(ic+1)) | |
3495 ic++; | |
3496 if (i < lhs_nz) | |
3497 ii = ic * lhs_nr + c_lhs.ridx(i); | |
3498 } | |
3499 j++; | |
3500 if (j < len) | |
3501 { | |
3502 jj = idx_i.elem (j); | |
3503 jr = jj % lhs_nr; | |
3504 jc = (jj - jr) / lhs_nr; | |
3505 } | |
3506 } | |
3507 } | |
3508 | |
5275 | 3509 for (octave_idx_type iidx = kc; iidx < lhs_nc+1; iidx++) |
5164 | 3510 stmp.xcidx(iidx) = kk; |
3511 | |
3512 lhs = stmp; | |
3513 } | |
3514 else | |
3515 { | |
3516 (*current_liboctave_error_handler) | |
3517 ("A(I) = X: X must be a scalar or a matrix with the same size as I"); | |
3518 | |
3519 retval = 0; | |
3520 } | |
3521 } | |
3522 // idx_vector::freeze() printed an error message for us. | |
3523 } | |
3524 } | |
3525 else | |
3526 { | |
3527 (*current_liboctave_error_handler) | |
3528 ("invalid number of indices for matrix expression"); | |
3529 | |
3530 retval = 0; | |
3531 } | |
3532 | |
3533 lhs.clear_index (); | |
3534 | |
3535 return retval; | |
3536 } | |
3537 | |
7356 | 3538 /* |
3539 * Tests | |
3540 * | |
3541 | |
7546
4249c6fb6e09
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7463
diff
changeset
|
3542 %!function x = set_slice(x, dim, slice, arg) |
7356 | 3543 %! switch dim |
3544 %! case 11 | |
3545 %! x(slice) = 2; | |
3546 %! case 21 | |
3547 %! x(slice, :) = 2; | |
3548 %! case 22 | |
3549 %! x(:, slice) = 2; | |
3550 %! otherwise | |
3551 %! error("invalid dim, '%d'", dim); | |
3552 %! endswitch | |
3553 %! endfunction | |
3554 | |
7546
4249c6fb6e09
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7463
diff
changeset
|
3555 %!function x = set_slice2(x, dim, slice) |
4249c6fb6e09
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7463
diff
changeset
|
3556 %! switch dim |
4249c6fb6e09
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7463
diff
changeset
|
3557 %! case 11 |
4249c6fb6e09
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7463
diff
changeset
|
3558 %! x(slice) = 2 * ones (size(slice)); |
4249c6fb6e09
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7463
diff
changeset
|
3559 %! case 21 |
4249c6fb6e09
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7463
diff
changeset
|
3560 %! x(slice, :) = 2 * ones (length(slice), columns (x)); |
4249c6fb6e09
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7463
diff
changeset
|
3561 %! case 22 |
4249c6fb6e09
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7463
diff
changeset
|
3562 %! x(:, slice) = 2 * ones (rows (x), length(slice)); |
4249c6fb6e09
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7463
diff
changeset
|
3563 %! otherwise |
4249c6fb6e09
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7463
diff
changeset
|
3564 %! error("invalid dim, '%d'", dim); |
4249c6fb6e09
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7463
diff
changeset
|
3565 %! endswitch |
4249c6fb6e09
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7463
diff
changeset
|
3566 %! endfunction |
4249c6fb6e09
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7463
diff
changeset
|
3567 |
7356 | 3568 %!function test_sparse_slice(size, dim, slice) |
3569 %! x = ones(size); | |
3570 %! s = set_slice(sparse(x), dim, slice); | |
3571 %! f = set_slice(x, dim, slice); | |
7546
4249c6fb6e09
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7463
diff
changeset
|
3572 %! assert (nnz(s), nnz(f)); |
4249c6fb6e09
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7463
diff
changeset
|
3573 %! assert(full(s), f); |
4249c6fb6e09
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7463
diff
changeset
|
3574 %! s = set_slice2(sparse(x), dim, slice); |
4249c6fb6e09
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7463
diff
changeset
|
3575 %! f = set_slice2(x, dim, slice); |
4249c6fb6e09
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7463
diff
changeset
|
3576 %! assert (nnz(s), nnz(f)); |
7356 | 3577 %! assert(full(s), f); |
3578 %! endfunction | |
3579 | |
3580 #### 1d indexing | |
3581 | |
3582 ## size = [2 0] | |
3583 %!test test_sparse_slice([2 0], 11, []); | |
3584 %!assert(set_slice(sparse(ones([2 0])), 11, 1), sparse([2 0]')); # sparse different from full | |
3585 %!assert(set_slice(sparse(ones([2 0])), 11, 2), sparse([0 2]')); # sparse different from full | |
3586 %!assert(set_slice(sparse(ones([2 0])), 11, 3), sparse([0 0 2]')); # sparse different from full | |
3587 %!assert(set_slice(sparse(ones([2 0])), 11, 4), sparse([0 0 0 2]')); # sparse different from full | |
3588 | |
3589 ## size = [0 2] | |
3590 %!test test_sparse_slice([0 2], 11, []); | |
3591 %!assert(set_slice(sparse(ones([0 2])), 11, 1), sparse(1,2)); # sparse different from full | |
3592 %!test test_sparse_slice([0 2], 11, 2); | |
3593 %!test test_sparse_slice([0 2], 11, 3); | |
3594 %!test test_sparse_slice([0 2], 11, 4); | |
7546
4249c6fb6e09
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7463
diff
changeset
|
3595 %!test test_sparse_slice([0 2], 11, [4, 4]); |
7356 | 3596 |
3597 ## size = [2 1] | |
3598 %!test test_sparse_slice([2 1], 11, []); | |
3599 %!test test_sparse_slice([2 1], 11, 1); | |
3600 %!test test_sparse_slice([2 1], 11, 2); | |
3601 %!test test_sparse_slice([2 1], 11, 3); | |
3602 %!test test_sparse_slice([2 1], 11, 4); | |
7546
4249c6fb6e09
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7463
diff
changeset
|
3603 %!test test_sparse_slice([2 1], 11, [4, 4]); |
7356 | 3604 |
3605 ## size = [1 2] | |
3606 %!test test_sparse_slice([1 2], 11, []); | |
3607 %!test test_sparse_slice([1 2], 11, 1); | |
3608 %!test test_sparse_slice([1 2], 11, 2); | |
3609 %!test test_sparse_slice([1 2], 11, 3); | |
3610 %!test test_sparse_slice([1 2], 11, 4); | |
7546
4249c6fb6e09
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7463
diff
changeset
|
3611 %!test test_sparse_slice([1 2], 11, [4, 4]); |
7356 | 3612 |
3613 ## size = [2 2] | |
3614 %!test test_sparse_slice([2 2], 11, []); | |
3615 %!test test_sparse_slice([2 2], 11, 1); | |
3616 %!test test_sparse_slice([2 2], 11, 2); | |
3617 %!test test_sparse_slice([2 2], 11, 3); | |
3618 %!test test_sparse_slice([2 2], 11, 4); | |
7546
4249c6fb6e09
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7463
diff
changeset
|
3619 %!test test_sparse_slice([2 2], 11, [4, 4]); |
7356 | 3620 # These 2 errors are the same as in the full case |
3621 %!error <invalid matrix index = 5> set_slice(sparse(ones([2 2])), 11, 5); | |
3622 %!error <invalid matrix index = 6> set_slice(sparse(ones([2 2])), 11, 6); | |
3623 | |
3624 | |
3625 #### 2d indexing | |
3626 | |
3627 ## size = [2 0] | |
3628 %!test test_sparse_slice([2 0], 21, []); | |
3629 %!test test_sparse_slice([2 0], 21, 1); | |
3630 %!test test_sparse_slice([2 0], 21, 2); | |
7546
4249c6fb6e09
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7463
diff
changeset
|
3631 %!test test_sparse_slice([2 0], 21, [2,2]); |
7356 | 3632 %!assert(set_slice(sparse(ones([2 0])), 21, 3), sparse(2,0)); # sparse different from full |
3633 %!assert(set_slice(sparse(ones([2 0])), 21, 4), sparse(2,0)); # sparse different from full | |
3634 %!test test_sparse_slice([2 0], 22, []); | |
3635 %!test test_sparse_slice([2 0], 22, 1); | |
3636 %!test test_sparse_slice([2 0], 22, 2); | |
7546
4249c6fb6e09
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7463
diff
changeset
|
3637 %!test test_sparse_slice([2 0], 22, [2,2]); |
7356 | 3638 %!assert(set_slice(sparse(ones([2 0])), 22, 3), sparse([0 0 2;0 0 2])); # sparse different from full |
3639 %!assert(set_slice(sparse(ones([2 0])), 22, 4), sparse([0 0 0 2;0 0 0 2])); # sparse different from full | |
3640 | |
3641 ## size = [0 2] | |
3642 %!test test_sparse_slice([0 2], 21, []); | |
3643 %!test test_sparse_slice([0 2], 21, 1); | |
3644 %!test test_sparse_slice([0 2], 21, 2); | |
7546
4249c6fb6e09
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7463
diff
changeset
|
3645 %!test test_sparse_slice([0 2], 21, [2,2]); |
7356 | 3646 %!assert(set_slice(sparse(ones([0 2])), 21, 3), sparse([0 0;0 0;2 2])); # sparse different from full |
3647 %!assert(set_slice(sparse(ones([0 2])), 21, 4), sparse([0 0;0 0;0 0;2 2])); # sparse different from full | |
3648 %!test test_sparse_slice([0 2], 22, []); | |
3649 %!test test_sparse_slice([0 2], 22, 1); | |
3650 %!test test_sparse_slice([0 2], 22, 2); | |
7546
4249c6fb6e09
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7463
diff
changeset
|
3651 %!test test_sparse_slice([0 2], 22, [2,2]); |
7356 | 3652 %!assert(set_slice(sparse(ones([0 2])), 22, 3), sparse(0,2)); # sparse different from full |
3653 %!assert(set_slice(sparse(ones([0 2])), 22, 4), sparse(0,2)); # sparse different from full | |
3654 | |
3655 ## size = [2 1] | |
3656 %!test test_sparse_slice([2 1], 21, []); | |
3657 %!test test_sparse_slice([2 1], 21, 1); | |
3658 %!test test_sparse_slice([2 1], 21, 2); | |
7546
4249c6fb6e09
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7463
diff
changeset
|
3659 %!test test_sparse_slice([2 1], 21, [2,2]); |
7356 | 3660 %!test test_sparse_slice([2 1], 21, 3); |
3661 %!test test_sparse_slice([2 1], 21, 4); | |
3662 %!test test_sparse_slice([2 1], 22, []); | |
3663 %!test test_sparse_slice([2 1], 22, 1); | |
3664 %!test test_sparse_slice([2 1], 22, 2); | |
7546
4249c6fb6e09
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7463
diff
changeset
|
3665 %!test test_sparse_slice([2 1], 22, [2,2]); |
7356 | 3666 %!test test_sparse_slice([2 1], 22, 3); |
3667 %!test test_sparse_slice([2 1], 22, 4); | |
3668 | |
3669 ## size = [1 2] | |
3670 %!test test_sparse_slice([1 2], 21, []); | |
3671 %!test test_sparse_slice([1 2], 21, 1); | |
3672 %!test test_sparse_slice([1 2], 21, 2); | |
7546
4249c6fb6e09
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7463
diff
changeset
|
3673 %!test test_sparse_slice([1 2], 21, [2,2]); |
7356 | 3674 %!test test_sparse_slice([1 2], 21, 3); |
3675 %!test test_sparse_slice([1 2], 21, 4); | |
3676 %!test test_sparse_slice([1 2], 22, []); | |
3677 %!test test_sparse_slice([1 2], 22, 1); | |
3678 %!test test_sparse_slice([1 2], 22, 2); | |
7546
4249c6fb6e09
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7463
diff
changeset
|
3679 %!test test_sparse_slice([1 2], 22, [2,2]); |
7356 | 3680 %!test test_sparse_slice([1 2], 22, 3); |
3681 %!test test_sparse_slice([1 2], 22, 4); | |
3682 | |
3683 ## size = [2 2] | |
3684 %!test test_sparse_slice([2 2], 21, []); | |
3685 %!test test_sparse_slice([2 2], 21, 1); | |
3686 %!test test_sparse_slice([2 2], 21, 2); | |
7546
4249c6fb6e09
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7463
diff
changeset
|
3687 %!test test_sparse_slice([2 2], 21, [2,2]); |
7356 | 3688 %!test test_sparse_slice([2 2], 21, 3); |
3689 %!test test_sparse_slice([2 2], 21, 4); | |
3690 %!test test_sparse_slice([2 2], 22, []); | |
3691 %!test test_sparse_slice([2 2], 22, 1); | |
3692 %!test test_sparse_slice([2 2], 22, 2); | |
7546
4249c6fb6e09
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7463
diff
changeset
|
3693 %!test test_sparse_slice([2 2], 22, [2,2]); |
7356 | 3694 %!test test_sparse_slice([2 2], 22, 3); |
3695 %!test test_sparse_slice([2 2], 22, 4); | |
3696 | |
3697 */ | |
3698 | |
5164 | 3699 template <class T> |
3700 void | |
3701 Sparse<T>::print_info (std::ostream& os, const std::string& prefix) const | |
3702 { | |
3703 os << prefix << "rep address: " << rep << "\n" | |
5604 | 3704 << prefix << "rep->nzmx: " << rep->nzmx << "\n" |
5164 | 3705 << prefix << "rep->nrows: " << rep->nrows << "\n" |
3706 << prefix << "rep->ncols: " << rep->ncols << "\n" | |
3707 << prefix << "rep->data: " << static_cast<void *> (rep->d) << "\n" | |
3708 << prefix << "rep->ridx: " << static_cast<void *> (rep->r) << "\n" | |
3709 << prefix << "rep->cidx: " << static_cast<void *> (rep->c) << "\n" | |
3710 << prefix << "rep->count: " << rep->count << "\n"; | |
3711 } | |
3712 | |
3713 /* | |
3714 ;;; Local Variables: *** | |
3715 ;;; mode: C++ *** | |
3716 ;;; End: *** | |
3717 */ |