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