Mercurial > hg > octave-nkf
annotate liboctave/Sparse.cc @ 11669:a6e08ecb4050 release-3-0-x
Treat repeated indices in the sparse assignments
| author | David Bateman <dbateman@free.fr> |
|---|---|
| date | Mon, 03 Mar 2008 20:00:17 -0500 |
| parents | fcc6d853df9e |
| children | 72830070a17b |
| 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 *> | |
| 372 sort (octave_sparse_sidxl_comp); | |
| 373 | |
| 5604 | 374 sort.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 *> | |
| 490 sort (octave_sparse_sidxl_comp); | |
| 491 | |
| 5604 | 492 sort.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)); | |
| 1858 octave_sort<octave_idx_type> sort; | |
| 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); | |
| 1874 sort.sort (itmp, n); | |
| 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); | |
| 1886 sort.sort (itmp, m); | |
| 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 } | |
| 1923 sort.sort (ri + retval.xcidx (j), kk - retval.xcidx (j)); | |
| 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 } | |
| 6988 | 2025 sort.sort (ri + retval.xcidx (j), |
| 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 | |
| 5775 | 2056 // FIXME |
| 5164 | 2057 // Unfortunately numel can overflow for very large but very sparse matrices. |
| 2058 // For now just flag an error when this happens. | |
| 2059 template <class LT, class RT> | |
| 2060 int | |
| 2061 assign1 (Sparse<LT>& lhs, const Sparse<RT>& rhs) | |
| 2062 { | |
| 2063 int retval = 1; | |
| 2064 | |
| 2065 idx_vector *idx_tmp = lhs.get_idx (); | |
| 2066 | |
| 2067 idx_vector lhs_idx = idx_tmp[0]; | |
| 2068 | |
| 5275 | 2069 octave_idx_type lhs_len = lhs.numel (); |
| 2070 octave_idx_type rhs_len = rhs.numel (); | |
| 5164 | 2071 |
| 5828 | 2072 uint64_t long_lhs_len = |
| 2073 static_cast<uint64_t> (lhs.rows ()) * | |
| 2074 static_cast<uint64_t> (lhs.cols ()); | |
| 2075 | |
| 2076 uint64_t long_rhs_len = | |
| 2077 static_cast<uint64_t> (rhs.rows ()) * | |
| 2078 static_cast<uint64_t> (rhs.cols ()); | |
| 2079 | |
| 2080 if (long_rhs_len != static_cast<uint64_t>(rhs_len) || | |
| 2081 long_lhs_len != static_cast<uint64_t>(lhs_len)) | |
| 5164 | 2082 { |
| 2083 (*current_liboctave_error_handler) | |
| 2084 ("A(I) = X: Matrix dimensions too large to ensure correct\n", | |
| 2085 "operation. This is an limitation that should be removed\n", | |
| 2086 "in the future."); | |
| 2087 | |
| 2088 lhs.clear_index (); | |
| 2089 return 0; | |
| 2090 } | |
| 2091 | |
| 5275 | 2092 octave_idx_type nr = lhs.rows (); |
| 2093 octave_idx_type nc = lhs.cols (); | |
| 5681 | 2094 octave_idx_type nz = lhs.nnz (); |
| 5275 | 2095 |
| 5781 | 2096 octave_idx_type n = lhs_idx.freeze (lhs_len, "vector", true); |
| 5164 | 2097 |
| 2098 if (n != 0) | |
| 2099 { | |
| 5275 | 2100 octave_idx_type max_idx = lhs_idx.max () + 1; |
| 5164 | 2101 max_idx = max_idx < lhs_len ? lhs_len : max_idx; |
| 2102 | |
| 2103 // Take a constant copy of lhs. This means that elem won't | |
| 2104 // create missing elements. | |
| 2105 const Sparse<LT> c_lhs (lhs); | |
| 2106 | |
| 2107 if (rhs_len == n) | |
| 2108 { | |
| 5681 | 2109 octave_idx_type new_nzmx = lhs.nnz (); |
| 5164 | 2110 |
| 5603 | 2111 OCTAVE_LOCAL_BUFFER (octave_idx_type, rhs_idx, n); |
| 2112 if (! lhs_idx.is_colon ()) | |
| 2113 { | |
| 2114 // Ok here we have to be careful with the indexing, | |
| 2115 // to treat cases like "a([3,2,1]) = b", and still | |
| 2116 // handle the need for strict sorting of the sparse | |
| 2117 // elements. | |
| 2118 OCTAVE_LOCAL_BUFFER (octave_idx_vector_sort *, sidx, n); | |
| 2119 OCTAVE_LOCAL_BUFFER (octave_idx_vector_sort, sidxX, n); | |
| 2120 | |
| 2121 for (octave_idx_type i = 0; i < n; i++) | |
| 2122 { | |
| 2123 sidx[i] = &sidxX[i]; | |
| 2124 sidx[i]->i = lhs_idx.elem(i); | |
| 2125 sidx[i]->idx = i; | |
| 2126 } | |
| 2127 | |
| 2128 OCTAVE_QUIT; | |
| 2129 octave_sort<octave_idx_vector_sort *> | |
| 2130 sort (octave_idx_vector_comp); | |
| 2131 | |
| 2132 sort.sort (sidx, n); | |
| 2133 | |
| 2134 intNDArray<octave_idx_type> new_idx (dim_vector (n,1)); | |
| 2135 | |
| 2136 for (octave_idx_type i = 0; i < n; i++) | |
| 2137 { | |
| 2138 new_idx.xelem(i) = sidx[i]->i + 1; | |
| 2139 rhs_idx[i] = sidx[i]->idx; | |
| 2140 } | |
| 2141 | |
| 2142 lhs_idx = idx_vector (new_idx); | |
| 2143 } | |
| 2144 else | |
| 2145 for (octave_idx_type i = 0; i < n; i++) | |
| 2146 rhs_idx[i] = i; | |
| 2147 | |
| 5164 | 2148 // First count the number of non-zero elements |
| 5275 | 2149 for (octave_idx_type i = 0; i < n; i++) |
| 5164 | 2150 { |
| 2151 OCTAVE_QUIT; | |
| 2152 | |
| 5275 | 2153 octave_idx_type ii = lhs_idx.elem (i); |
|
11669
a6e08ecb4050
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7326
diff
changeset
|
2154 if (i < n - 1 && lhs_idx.elem (i + 1) == ii) |
|
a6e08ecb4050
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7326
diff
changeset
|
2155 continue; |
| 5164 | 2156 if (ii < lhs_len && c_lhs.elem(ii) != LT ()) |
| 5604 | 2157 new_nzmx--; |
| 5603 | 2158 if (rhs.elem(rhs_idx[i]) != RT ()) |
| 5604 | 2159 new_nzmx++; |
| 5164 | 2160 } |
| 2161 | |
| 2162 if (nr > 1) | |
| 2163 { | |
| 7246 | 2164 Sparse<LT> tmp ((max_idx > nr ? max_idx : nr), 1, new_nzmx); |
| 5164 | 2165 tmp.cidx(0) = 0; |
| 5681 | 2166 tmp.cidx(1) = new_nzmx; |
| 5164 | 2167 |
| 5275 | 2168 octave_idx_type i = 0; |
| 2169 octave_idx_type ii = 0; | |
| 5164 | 2170 if (i < nz) |
| 2171 ii = c_lhs.ridx(i); | |
| 2172 | |
| 5275 | 2173 octave_idx_type j = 0; |
| 2174 octave_idx_type jj = lhs_idx.elem(j); | |
| 2175 | |
| 2176 octave_idx_type kk = 0; | |
| 5164 | 2177 |
| 2178 while (j < n || i < nz) | |
| 2179 { | |
| 2180 if (j == n || (i < nz && ii < jj)) | |
| 2181 { | |
| 2182 tmp.xdata (kk) = c_lhs.data (i); | |
| 2183 tmp.xridx (kk++) = ii; | |
| 2184 if (++i < nz) | |
| 2185 ii = c_lhs.ridx(i); | |
| 2186 } | |
| 2187 else | |
| 2188 { | |
| 5603 | 2189 RT rtmp = rhs.elem (rhs_idx[j]); |
| 5164 | 2190 if (rtmp != RT ()) |
| 2191 { | |
| 2192 tmp.xdata (kk) = rtmp; | |
| 2193 tmp.xridx (kk++) = jj; | |
| 2194 } | |
| 2195 | |
| 2196 if (ii == jj && i < nz) | |
| 2197 if (++i < nz) | |
| 2198 ii = c_lhs.ridx(i); | |
| 2199 if (++j < n) | |
| 2200 jj = lhs_idx.elem(j); | |
| 2201 } | |
| 2202 } | |
| 2203 | |
| 2204 lhs = tmp; | |
| 2205 } | |
| 2206 else | |
| 2207 { | |
| 7246 | 2208 Sparse<LT> tmp (1, (max_idx > nc ? max_idx : nc), new_nzmx); |
| 5164 | 2209 |
| 5275 | 2210 octave_idx_type i = 0; |
| 2211 octave_idx_type ii = 0; | |
| 5164 | 2212 while (ii < nc && c_lhs.cidx(ii+1) <= i) |
| 2213 ii++; | |
| 2214 | |
| 5275 | 2215 octave_idx_type j = 0; |
| 2216 octave_idx_type jj = lhs_idx.elem(j); | |
| 2217 | |
| 2218 octave_idx_type kk = 0; | |
| 2219 octave_idx_type ic = 0; | |
| 5164 | 2220 |
| 2221 while (j < n || i < nz) | |
| 2222 { | |
| 2223 if (j == n || (i < nz && ii < jj)) | |
| 2224 { | |
| 2225 while (ic <= ii) | |
| 2226 tmp.xcidx (ic++) = kk; | |
| 2227 tmp.xdata (kk) = c_lhs.data (i); | |
| 5603 | 2228 tmp.xridx (kk++) = 0; |
| 5164 | 2229 i++; |
| 2230 while (ii < nc && c_lhs.cidx(ii+1) <= i) | |
| 2231 ii++; | |
| 2232 } | |
| 2233 else | |
| 2234 { | |
| 2235 while (ic <= jj) | |
| 2236 tmp.xcidx (ic++) = kk; | |
| 2237 | |
| 5603 | 2238 RT rtmp = rhs.elem (rhs_idx[j]); |
| 5164 | 2239 if (rtmp != RT ()) |
| 5603 | 2240 { |
| 2241 tmp.xdata (kk) = rtmp; | |
| 2242 tmp.xridx (kk++) = 0; | |
| 2243 } | |
| 5164 | 2244 if (ii == jj) |
| 2245 { | |
| 2246 i++; | |
| 2247 while (ii < nc && c_lhs.cidx(ii+1) <= i) | |
| 2248 ii++; | |
| 2249 } | |
| 2250 j++; | |
| 2251 if (j < n) | |
| 2252 jj = lhs_idx.elem(j); | |
| 2253 } | |
| 2254 } | |
| 2255 | |
| 5275 | 2256 for (octave_idx_type iidx = ic; iidx < max_idx+1; iidx++) |
| 5164 | 2257 tmp.xcidx(iidx) = kk; |
| 2258 | |
| 2259 lhs = tmp; | |
| 2260 } | |
| 2261 } | |
| 2262 else if (rhs_len == 1) | |
| 2263 { | |
| 5681 | 2264 octave_idx_type new_nzmx = lhs.nnz (); |
| 5164 | 2265 RT scalar = rhs.elem (0); |
| 2266 bool scalar_non_zero = (scalar != RT ()); | |
| 5603 | 2267 lhs_idx.sort (true); |
|
11669
a6e08ecb4050
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7326
diff
changeset
|
2268 n = lhs_idx.length (n); |
| 5164 | 2269 |
| 2270 // First count the number of non-zero elements | |
| 2271 if (scalar != RT ()) | |
| 5604 | 2272 new_nzmx += n; |
| 5275 | 2273 for (octave_idx_type i = 0; i < n; i++) |
| 5164 | 2274 { |
| 2275 OCTAVE_QUIT; | |
| 2276 | |
| 5275 | 2277 octave_idx_type ii = lhs_idx.elem (i); |
| 5164 | 2278 if (ii < lhs_len && c_lhs.elem(ii) != LT ()) |
| 5604 | 2279 new_nzmx--; |
| 5164 | 2280 } |
| 2281 | |
| 2282 if (nr > 1) | |
| 2283 { | |
| 7246 | 2284 Sparse<LT> tmp ((max_idx > nr ? max_idx : nr), 1, new_nzmx); |
| 5164 | 2285 tmp.cidx(0) = 0; |
| 5681 | 2286 tmp.cidx(1) = new_nzmx; |
| 5164 | 2287 |
| 5275 | 2288 octave_idx_type i = 0; |
| 2289 octave_idx_type ii = 0; | |
| 5164 | 2290 if (i < nz) |
| 2291 ii = c_lhs.ridx(i); | |
| 2292 | |
| 5275 | 2293 octave_idx_type j = 0; |
| 2294 octave_idx_type jj = lhs_idx.elem(j); | |
| 2295 | |
| 2296 octave_idx_type kk = 0; | |
| 5164 | 2297 |
| 2298 while (j < n || i < nz) | |
| 2299 { | |
|
11669
a6e08ecb4050
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7326
diff
changeset
|
2300 if (j < n - 1 && lhs_idx.elem (j + 1) == jj) |
|
a6e08ecb4050
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7326
diff
changeset
|
2301 { |
|
a6e08ecb4050
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7326
diff
changeset
|
2302 j++; |
|
a6e08ecb4050
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7326
diff
changeset
|
2303 jj = lhs_idx.elem (j); |
|
a6e08ecb4050
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7326
diff
changeset
|
2304 continue; |
|
a6e08ecb4050
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7326
diff
changeset
|
2305 } |
| 5164 | 2306 if (j == n || (i < nz && ii < jj)) |
| 2307 { | |
| 2308 tmp.xdata (kk) = c_lhs.data (i); | |
| 2309 tmp.xridx (kk++) = ii; | |
| 2310 if (++i < nz) | |
| 2311 ii = c_lhs.ridx(i); | |
| 2312 } | |
| 2313 else | |
| 2314 { | |
| 2315 if (scalar_non_zero) | |
| 2316 { | |
| 2317 tmp.xdata (kk) = scalar; | |
| 2318 tmp.xridx (kk++) = jj; | |
| 2319 } | |
| 2320 | |
| 2321 if (ii == jj && i < nz) | |
| 2322 if (++i < nz) | |
| 2323 ii = c_lhs.ridx(i); | |
| 2324 if (++j < n) | |
| 2325 jj = lhs_idx.elem(j); | |
| 2326 } | |
| 2327 } | |
| 2328 | |
| 2329 lhs = tmp; | |
| 2330 } | |
| 2331 else | |
| 2332 { | |
| 7246 | 2333 Sparse<LT> tmp (1, (max_idx > nc ? max_idx : nc), new_nzmx); |
| 5164 | 2334 |
| 5275 | 2335 octave_idx_type i = 0; |
| 2336 octave_idx_type ii = 0; | |
| 5164 | 2337 while (ii < nc && c_lhs.cidx(ii+1) <= i) |
| 2338 ii++; | |
| 2339 | |
| 5275 | 2340 octave_idx_type j = 0; |
| 2341 octave_idx_type jj = lhs_idx.elem(j); | |
| 2342 | |
| 2343 octave_idx_type kk = 0; | |
| 2344 octave_idx_type ic = 0; | |
| 5164 | 2345 |
| 2346 while (j < n || i < nz) | |
| 2347 { | |
|
11669
a6e08ecb4050
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7326
diff
changeset
|
2348 if (j < n - 1 && lhs_idx.elem (j + 1) == jj) |
|
a6e08ecb4050
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7326
diff
changeset
|
2349 { |
|
a6e08ecb4050
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7326
diff
changeset
|
2350 j++; |
|
a6e08ecb4050
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7326
diff
changeset
|
2351 jj = lhs_idx.elem (j); |
|
a6e08ecb4050
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7326
diff
changeset
|
2352 continue; |
|
a6e08ecb4050
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7326
diff
changeset
|
2353 } |
| 5164 | 2354 if (j == n || (i < nz && ii < jj)) |
| 2355 { | |
| 2356 while (ic <= ii) | |
| 2357 tmp.xcidx (ic++) = kk; | |
| 2358 tmp.xdata (kk) = c_lhs.data (i); | |
| 2359 i++; | |
| 2360 while (ii < nc && c_lhs.cidx(ii+1) <= i) | |
| 2361 ii++; | |
| 2362 } | |
| 2363 else | |
| 2364 { | |
| 2365 while (ic <= jj) | |
| 2366 tmp.xcidx (ic++) = kk; | |
| 2367 if (scalar_non_zero) | |
| 2368 tmp.xdata (kk) = scalar; | |
| 2369 if (ii == jj) | |
| 2370 { | |
| 2371 i++; | |
| 2372 while (ii < nc && c_lhs.cidx(ii+1) <= i) | |
| 2373 ii++; | |
| 2374 } | |
| 2375 j++; | |
| 2376 if (j < n) | |
| 2377 jj = lhs_idx.elem(j); | |
| 2378 } | |
| 2379 tmp.xridx (kk++) = 0; | |
| 2380 } | |
| 2381 | |
| 5275 | 2382 for (octave_idx_type iidx = ic; iidx < max_idx+1; iidx++) |
| 5164 | 2383 tmp.xcidx(iidx) = kk; |
| 2384 | |
| 2385 lhs = tmp; | |
| 2386 } | |
| 2387 } | |
| 2388 else | |
| 2389 { | |
| 2390 (*current_liboctave_error_handler) | |
| 2391 ("A(I) = X: X must be a scalar or a vector with same length as I"); | |
| 2392 | |
| 2393 retval = 0; | |
| 2394 } | |
| 2395 } | |
| 2396 else if (lhs_idx.is_colon ()) | |
| 2397 { | |
| 2398 if (lhs_len == 0) | |
| 2399 { | |
| 2400 | |
| 5681 | 2401 octave_idx_type new_nzmx = rhs.nnz (); |
| 5604 | 2402 Sparse<LT> tmp (1, rhs_len, new_nzmx); |
| 5164 | 2403 |
| 5275 | 2404 octave_idx_type ii = 0; |
| 2405 octave_idx_type jj = 0; | |
| 2406 for (octave_idx_type i = 0; i < rhs.cols(); i++) | |
| 2407 for (octave_idx_type j = rhs.cidx(i); j < rhs.cidx(i+1); j++) | |
| 5164 | 2408 { |
| 2409 OCTAVE_QUIT; | |
| 5275 | 2410 for (octave_idx_type k = jj; k <= i * rhs.rows() + rhs.ridx(j); k++) |
| 5164 | 2411 tmp.cidx(jj++) = ii; |
| 2412 | |
| 2413 tmp.data(ii) = rhs.data(j); | |
| 2414 tmp.ridx(ii++) = 0; | |
| 2415 } | |
| 2416 | |
| 5275 | 2417 for (octave_idx_type i = jj; i < rhs_len + 1; i++) |
| 5164 | 2418 tmp.cidx(i) = ii; |
| 2419 | |
| 2420 lhs = tmp; | |
| 2421 } | |
| 2422 else | |
| 2423 (*current_liboctave_error_handler) | |
| 2424 ("A(:) = X: A must be the same size as X"); | |
| 2425 } | |
| 2426 else if (! (rhs_len == 1 || rhs_len == 0)) | |
| 2427 { | |
| 2428 (*current_liboctave_error_handler) | |
| 2429 ("A([]) = X: X must also be an empty matrix or a scalar"); | |
| 2430 | |
| 2431 retval = 0; | |
| 2432 } | |
| 2433 | |
| 2434 lhs.clear_index (); | |
| 2435 | |
| 2436 return retval; | |
| 2437 } | |
| 2438 | |
| 2439 template <class LT, class RT> | |
| 2440 int | |
| 2441 assign (Sparse<LT>& lhs, const Sparse<RT>& rhs) | |
| 2442 { | |
| 2443 int retval = 1; | |
| 2444 | |
| 2445 int n_idx = lhs.index_count (); | |
| 2446 | |
| 5275 | 2447 octave_idx_type lhs_nr = lhs.rows (); |
| 2448 octave_idx_type lhs_nc = lhs.cols (); | |
| 5681 | 2449 octave_idx_type lhs_nz = lhs.nnz (); |
| 5275 | 2450 |
| 2451 octave_idx_type rhs_nr = rhs.rows (); | |
| 2452 octave_idx_type rhs_nc = rhs.cols (); | |
| 5164 | 2453 |
| 2454 idx_vector *tmp = lhs.get_idx (); | |
| 2455 | |
| 2456 idx_vector idx_i; | |
| 2457 idx_vector idx_j; | |
| 2458 | |
| 2459 if (n_idx > 2) | |
| 2460 { | |
| 2461 (*current_liboctave_error_handler) | |
| 2462 ("A(I, J) = X: can only have 1 or 2 indexes for sparse matrices"); | |
| 6092 | 2463 |
| 2464 lhs.clear_index (); | |
| 5164 | 2465 return 0; |
| 2466 } | |
| 2467 | |
| 2468 if (n_idx > 1) | |
| 2469 idx_j = tmp[1]; | |
| 2470 | |
| 2471 if (n_idx > 0) | |
| 2472 idx_i = tmp[0]; | |
| 2473 | |
| 6988 | 2474 // Take a constant copy of lhs. This means that ridx and family won't |
| 2475 // call make_unique. | |
| 2476 const Sparse<LT> c_lhs (lhs); | |
| 2477 | |
| 5164 | 2478 if (n_idx == 2) |
| 2479 { | |
| 5781 | 2480 octave_idx_type n = idx_i.freeze (lhs_nr, "row", true); |
| 2481 octave_idx_type m = idx_j.freeze (lhs_nc, "column", true); | |
| 5164 | 2482 |
| 2483 int idx_i_is_colon = idx_i.is_colon (); | |
| 2484 int idx_j_is_colon = idx_j.is_colon (); | |
| 2485 | |
| 7238 | 2486 if (lhs_nr == 0 && lhs_nc == 0) |
| 2487 { | |
| 2488 if (idx_i_is_colon) | |
| 2489 n = rhs_nr; | |
| 2490 | |
| 2491 if (idx_j_is_colon) | |
| 2492 m = rhs_nc; | |
| 2493 } | |
| 5164 | 2494 |
| 2495 if (idx_i && idx_j) | |
| 2496 { | |
| 2497 if (rhs_nr == 0 && rhs_nc == 0) | |
| 2498 { | |
| 2499 lhs.maybe_delete_elements (idx_i, idx_j); | |
| 2500 } | |
| 2501 else | |
| 2502 { | |
| 2503 if (rhs_nr == 1 && rhs_nc == 1 && n >= 0 && m >= 0) | |
| 2504 { | |
| 2505 if (n > 0 && m > 0) | |
| 2506 { | |
| 5603 | 2507 idx_i.sort (true); |
|
11669
a6e08ecb4050
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7326
diff
changeset
|
2508 n = idx_i.length (n); |
| 5603 | 2509 idx_j.sort (true); |
|
11669
a6e08ecb4050
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7326
diff
changeset
|
2510 m = idx_j.length (m); |
| 5603 | 2511 |
| 5275 | 2512 octave_idx_type max_row_idx = idx_i_is_colon ? rhs_nr : |
| 5164 | 2513 idx_i.max () + 1; |
| 5275 | 2514 octave_idx_type max_col_idx = idx_j_is_colon ? rhs_nc : |
| 5164 | 2515 idx_j.max () + 1; |
| 5603 | 2516 octave_idx_type new_nr = max_row_idx > lhs_nr ? |
| 2517 max_row_idx : lhs_nr; | |
| 2518 octave_idx_type new_nc = max_col_idx > lhs_nc ? | |
| 2519 max_col_idx : lhs_nc; | |
| 5164 | 2520 RT scalar = rhs.elem (0, 0); |
| 2521 | |
| 2522 // Count the number of non-zero terms | |
| 5681 | 2523 octave_idx_type new_nzmx = lhs.nnz (); |
| 5275 | 2524 for (octave_idx_type j = 0; j < m; j++) |
| 5164 | 2525 { |
| 5275 | 2526 octave_idx_type jj = idx_j.elem (j); |
| 5164 | 2527 if (jj < lhs_nc) |
| 2528 { | |
| 5275 | 2529 for (octave_idx_type i = 0; i < n; i++) |
| 5164 | 2530 { |
| 2531 OCTAVE_QUIT; | |
| 2532 | |
| 5275 | 2533 octave_idx_type ii = idx_i.elem (i); |
| 5164 | 2534 |
| 2535 if (ii < lhs_nr) | |
| 2536 { | |
| 6988 | 2537 for (octave_idx_type k = c_lhs.cidx(jj); |
| 2538 k < c_lhs.cidx(jj+1); k++) | |
| 5164 | 2539 { |
| 6988 | 2540 if (c_lhs.ridx(k) == ii) |
| 5604 | 2541 new_nzmx--; |
| 6988 | 2542 if (c_lhs.ridx(k) >= ii) |
| 5164 | 2543 break; |
| 2544 } | |
| 2545 } | |
| 2546 } | |
| 2547 } | |
| 2548 } | |
| 2549 | |
| 2550 if (scalar != RT()) | |
| 5604 | 2551 new_nzmx += m * n; |
| 2552 | |
| 2553 Sparse<LT> stmp (new_nr, new_nc, new_nzmx); | |
| 5164 | 2554 |
| 5275 | 2555 octave_idx_type jji = 0; |
| 2556 octave_idx_type jj = idx_j.elem (jji); | |
| 2557 octave_idx_type kk = 0; | |
| 5164 | 2558 stmp.cidx(0) = 0; |
| 5275 | 2559 for (octave_idx_type j = 0; j < new_nc; j++) |
| 5164 | 2560 { |
| 2561 if (jji < m && jj == j) | |
| 2562 { | |
| 5275 | 2563 octave_idx_type iii = 0; |
| 2564 octave_idx_type ii = idx_i.elem (iii); | |
| 5760 | 2565 octave_idx_type ppp = 0; |
| 6092 | 2566 octave_idx_type ppi = (j >= lhs_nc ? 0 : |
| 6988 | 2567 c_lhs.cidx(j+1) - |
| 2568 c_lhs.cidx(j)); | |
| 5760 | 2569 octave_idx_type pp = (ppp < ppi ? |
| 6988 | 2570 c_lhs.ridx(c_lhs.cidx(j)+ppp) : |
| 2571 new_nr); | |
| 5760 | 2572 while (ppp < ppi || iii < n) |
| 5164 | 2573 { |
| 5760 | 2574 if (iii < n && ii <= pp) |
| 5164 | 2575 { |
| 2576 if (scalar != RT ()) | |
| 2577 { | |
| 2578 stmp.data(kk) = scalar; | |
| 5760 | 2579 stmp.ridx(kk++) = ii; |
| 5164 | 2580 } |
| 5760 | 2581 if (ii == pp) |
| 6988 | 2582 pp = (++ppp < ppi ? c_lhs.ridx(c_lhs.cidx(j)+ppp) : new_nr); |
| 5164 | 2583 if (++iii < n) |
| 2584 ii = idx_i.elem(iii); | |
| 2585 } | |
| 5760 | 2586 else |
| 5164 | 2587 { |
| 5760 | 2588 stmp.data(kk) = |
| 6988 | 2589 c_lhs.data(c_lhs.cidx(j)+ppp); |
| 5760 | 2590 stmp.ridx(kk++) = pp; |
| 6988 | 2591 pp = (++ppp < ppi ? c_lhs.ridx(c_lhs.cidx(j)+ppp) : new_nr); |
| 5164 | 2592 } |
| 2593 } | |
| 2594 if (++jji < m) | |
| 2595 jj = idx_j.elem(jji); | |
| 2596 } | |
| 6988 | 2597 else if (j < lhs_nc) |
| 5164 | 2598 { |
| 6988 | 2599 for (octave_idx_type i = c_lhs.cidx(j); |
| 2600 i < c_lhs.cidx(j+1); i++) | |
| 5164 | 2601 { |
| 6988 | 2602 stmp.data(kk) = c_lhs.data(i); |
| 2603 stmp.ridx(kk++) = c_lhs.ridx(i); | |
| 5164 | 2604 } |
| 2605 } | |
| 2606 stmp.cidx(j+1) = kk; | |
| 2607 } | |
| 2608 | |
| 2609 lhs = stmp; | |
| 2610 } | |
| 7246 | 2611 else |
| 2612 { | |
| 7253 | 2613 #if 0 |
| 2614 // FIXME -- the following code will make this | |
| 2615 // function behave the same as the full matrix | |
| 2616 // case for things like | |
| 2617 // | |
| 2618 // x = sparse (ones (2)); | |
| 2619 // x([],3) = 2; | |
| 2620 // | |
| 2621 // x = | |
| 2622 // | |
| 2623 // Compressed Column Sparse (rows = 2, cols = 3, nnz = 4) | |
| 2624 // | |
| 2625 // (1, 1) -> 1 | |
| 2626 // (2, 1) -> 1 | |
| 2627 // (1, 2) -> 1 | |
| 2628 // (2, 2) -> 1 | |
| 2629 // | |
| 2630 // However, Matlab doesn't resize in this case | |
| 2631 // even though it does in the full matrix case. | |
| 2632 | |
| 7246 | 2633 if (n > 0) |
| 2634 { | |
| 2635 octave_idx_type max_row_idx = idx_i_is_colon ? | |
| 2636 rhs_nr : idx_i.max () + 1; | |
| 2637 octave_idx_type new_nr = max_row_idx > lhs_nr ? | |
| 2638 max_row_idx : lhs_nr; | |
| 2639 octave_idx_type new_nc = lhs_nc; | |
| 2640 | |
| 7253 | 2641 lhs.resize (new_nr, new_nc); |
| 7246 | 2642 } |
| 2643 else if (m > 0) | |
| 2644 { | |
| 2645 octave_idx_type max_col_idx = idx_j_is_colon ? | |
| 2646 rhs_nc : idx_j.max () + 1; | |
| 2647 octave_idx_type new_nr = lhs_nr; | |
| 2648 octave_idx_type new_nc = max_col_idx > lhs_nc ? | |
| 2649 max_col_idx : lhs_nc; | |
| 2650 | |
| 7253 | 2651 lhs.resize (new_nr, new_nc); |
| 7246 | 2652 } |
| 7253 | 2653 #endif |
| 7246 | 2654 } |
| 5164 | 2655 } |
| 2656 else if (n == rhs_nr && m == rhs_nc) | |
| 2657 { | |
| 2658 if (n > 0 && m > 0) | |
| 2659 { | |
| 5275 | 2660 octave_idx_type max_row_idx = idx_i_is_colon ? rhs_nr : |
| 5164 | 2661 idx_i.max () + 1; |
| 5275 | 2662 octave_idx_type max_col_idx = idx_j_is_colon ? rhs_nc : |
| 5164 | 2663 idx_j.max () + 1; |
| 5603 | 2664 octave_idx_type new_nr = max_row_idx > lhs_nr ? |
| 2665 max_row_idx : lhs_nr; | |
| 2666 octave_idx_type new_nc = max_col_idx > lhs_nc ? | |
| 2667 max_col_idx : lhs_nc; | |
| 2668 | |
| 2669 OCTAVE_LOCAL_BUFFER (octave_idx_type, rhs_idx_i, n); | |
| 2670 if (! idx_i.is_colon ()) | |
| 2671 { | |
| 2672 // Ok here we have to be careful with the indexing, | |
| 2673 // to treat cases like "a([3,2,1],:) = b", and still | |
| 2674 // handle the need for strict sorting of the sparse | |
| 2675 // elements. | |
| 2676 OCTAVE_LOCAL_BUFFER (octave_idx_vector_sort *, | |
| 2677 sidx, n); | |
| 2678 OCTAVE_LOCAL_BUFFER (octave_idx_vector_sort, | |
| 2679 sidxX, n); | |
| 2680 | |
| 2681 for (octave_idx_type i = 0; i < n; i++) | |
| 2682 { | |
| 2683 sidx[i] = &sidxX[i]; | |
| 2684 sidx[i]->i = idx_i.elem(i); | |
| 2685 sidx[i]->idx = i; | |
| 2686 } | |
| 2687 | |
| 2688 OCTAVE_QUIT; | |
| 2689 octave_sort<octave_idx_vector_sort *> | |
| 2690 sort (octave_idx_vector_comp); | |
| 2691 | |
| 2692 sort.sort (sidx, n); | |
| 2693 | |
| 2694 intNDArray<octave_idx_type> new_idx (dim_vector (n,1)); | |
| 2695 | |
| 2696 for (octave_idx_type i = 0; i < n; i++) | |
| 2697 { | |
| 2698 new_idx.xelem(i) = sidx[i]->i + 1; | |
| 2699 rhs_idx_i[i] = sidx[i]->idx; | |
| 2700 } | |
| 2701 | |
| 2702 idx_i = idx_vector (new_idx); | |
| 2703 } | |
| 2704 else | |
| 2705 for (octave_idx_type i = 0; i < n; i++) | |
| 2706 rhs_idx_i[i] = i; | |
| 2707 | |
| 2708 OCTAVE_LOCAL_BUFFER (octave_idx_type, rhs_idx_j, m); | |
| 2709 if (! idx_j.is_colon ()) | |
| 2710 { | |
| 2711 // Ok here we have to be careful with the indexing, | |
| 2712 // to treat cases like "a([3,2,1],:) = b", and still | |
| 2713 // handle the need for strict sorting of the sparse | |
| 2714 // elements. | |
| 2715 OCTAVE_LOCAL_BUFFER (octave_idx_vector_sort *, | |
| 2716 sidx, m); | |
| 2717 OCTAVE_LOCAL_BUFFER (octave_idx_vector_sort, | |
| 2718 sidxX, m); | |
| 2719 | |
| 2720 for (octave_idx_type i = 0; i < m; i++) | |
| 2721 { | |
| 2722 sidx[i] = &sidxX[i]; | |
| 2723 sidx[i]->i = idx_j.elem(i); | |
| 2724 sidx[i]->idx = i; | |
| 2725 } | |
| 2726 | |
| 2727 OCTAVE_QUIT; | |
| 2728 octave_sort<octave_idx_vector_sort *> | |
| 2729 sort (octave_idx_vector_comp); | |
| 2730 | |
| 2731 sort.sort (sidx, m); | |
| 2732 | |
| 2733 intNDArray<octave_idx_type> new_idx (dim_vector (m,1)); | |
| 2734 | |
| 2735 for (octave_idx_type i = 0; i < m; i++) | |
| 2736 { | |
| 2737 new_idx.xelem(i) = sidx[i]->i + 1; | |
| 2738 rhs_idx_j[i] = sidx[i]->idx; | |
| 2739 } | |
| 2740 | |
| 2741 idx_j = idx_vector (new_idx); | |
| 2742 } | |
| 2743 else | |
| 2744 for (octave_idx_type i = 0; i < m; i++) | |
| 2745 rhs_idx_j[i] = i; | |
| 5164 | 2746 |
| 5760 | 2747 // Maximum number of non-zero elements |
| 2748 octave_idx_type new_nzmx = lhs.nnz() + rhs.nnz(); | |
| 5164 | 2749 |
| 5604 | 2750 Sparse<LT> stmp (new_nr, new_nc, new_nzmx); |
| 5164 | 2751 |
| 5275 | 2752 octave_idx_type jji = 0; |
| 2753 octave_idx_type jj = idx_j.elem (jji); | |
| 2754 octave_idx_type kk = 0; | |
| 5164 | 2755 stmp.cidx(0) = 0; |
| 5275 | 2756 for (octave_idx_type j = 0; j < new_nc; j++) |
| 5164 | 2757 { |
| 2758 if (jji < m && jj == j) | |
| 2759 { | |
| 5275 | 2760 octave_idx_type iii = 0; |
| 2761 octave_idx_type ii = idx_i.elem (iii); | |
| 5760 | 2762 octave_idx_type ppp = 0; |
| 6092 | 2763 octave_idx_type ppi = (j >= lhs_nc ? 0 : |
| 6988 | 2764 c_lhs.cidx(j+1) - |
| 2765 c_lhs.cidx(j)); | |
| 5760 | 2766 octave_idx_type pp = (ppp < ppi ? |
| 6988 | 2767 c_lhs.ridx(c_lhs.cidx(j)+ppp) : |
| 2768 new_nr); | |
| 5760 | 2769 while (ppp < ppi || iii < n) |
| 5164 | 2770 { |
| 5760 | 2771 if (iii < n && ii <= pp) |
| 5164 | 2772 { |
|
11669
a6e08ecb4050
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7326
diff
changeset
|
2773 if (iii < n - 1 && |
|
a6e08ecb4050
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7326
diff
changeset
|
2774 idx_i.elem (iii + 1) == ii) |
|
a6e08ecb4050
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7326
diff
changeset
|
2775 { |
|
a6e08ecb4050
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7326
diff
changeset
|
2776 iii++; |
|
a6e08ecb4050
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7326
diff
changeset
|
2777 ii = idx_i.elem(iii); |
|
a6e08ecb4050
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7326
diff
changeset
|
2778 continue; |
|
a6e08ecb4050
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7326
diff
changeset
|
2779 } |
|
a6e08ecb4050
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7326
diff
changeset
|
2780 |
| 5603 | 2781 RT rtmp = rhs.elem (rhs_idx_i[iii], |
| 2782 rhs_idx_j[jji]); | |
| 5164 | 2783 if (rtmp != RT ()) |
| 2784 { | |
| 2785 stmp.data(kk) = rtmp; | |
| 5760 | 2786 stmp.ridx(kk++) = ii; |
| 5164 | 2787 } |
| 5760 | 2788 if (ii == pp) |
| 6988 | 2789 pp = (++ppp < ppi ? c_lhs.ridx(c_lhs.cidx(j)+ppp) : new_nr); |
| 5164 | 2790 if (++iii < n) |
| 2791 ii = idx_i.elem(iii); | |
| 2792 } | |
| 5760 | 2793 else |
| 5164 | 2794 { |
| 5760 | 2795 stmp.data(kk) = |
| 6988 | 2796 c_lhs.data(c_lhs.cidx(j)+ppp); |
| 5760 | 2797 stmp.ridx(kk++) = pp; |
| 6988 | 2798 pp = (++ppp < ppi ? c_lhs.ridx(c_lhs.cidx(j)+ppp) : new_nr); |
| 5164 | 2799 } |
| 2800 } | |
| 2801 if (++jji < m) | |
| 2802 jj = idx_j.elem(jji); | |
| 2803 } | |
| 6988 | 2804 else if (j < lhs_nc) |
| 5164 | 2805 { |
| 6988 | 2806 for (octave_idx_type i = c_lhs.cidx(j); |
| 2807 i < c_lhs.cidx(j+1); i++) | |
| 5164 | 2808 { |
| 6988 | 2809 stmp.data(kk) = c_lhs.data(i); |
| 2810 stmp.ridx(kk++) = c_lhs.ridx(i); | |
| 5164 | 2811 } |
| 2812 } | |
| 2813 stmp.cidx(j+1) = kk; | |
| 2814 } | |
| 2815 | |
| 5760 | 2816 stmp.maybe_compress(); |
| 5164 | 2817 lhs = stmp; |
| 2818 } | |
| 2819 } | |
| 2820 else if (n == 0 && m == 0) | |
| 2821 { | |
| 2822 if (! ((rhs_nr == 1 && rhs_nc == 1) | |
| 2823 || (rhs_nr == 0 || rhs_nc == 0))) | |
| 2824 { | |
| 2825 (*current_liboctave_error_handler) | |
| 2826 ("A([], []) = X: X must be an empty matrix or a scalar"); | |
| 2827 | |
| 2828 retval = 0; | |
| 2829 } | |
| 2830 } | |
| 2831 else | |
| 2832 { | |
| 2833 (*current_liboctave_error_handler) | |
| 2834 ("A(I, J) = X: X must be a scalar or the number of elements in I must"); | |
| 2835 (*current_liboctave_error_handler) | |
| 2836 ("match the number of rows in X and the number of elements in J must"); | |
| 2837 (*current_liboctave_error_handler) | |
| 2838 ("match the number of columns in X"); | |
| 2839 | |
| 2840 retval = 0; | |
| 2841 } | |
| 2842 } | |
| 2843 } | |
| 2844 // idx_vector::freeze() printed an error message for us. | |
| 2845 } | |
| 2846 else if (n_idx == 1) | |
| 2847 { | |
| 2848 int lhs_is_empty = lhs_nr == 0 || lhs_nc == 0; | |
| 2849 | |
| 2850 if (lhs_is_empty || (lhs_nr == 1 && lhs_nc == 1)) | |
| 2851 { | |
| 5275 | 2852 octave_idx_type lhs_len = lhs.length (); |
| 2853 | |
| 5781 | 2854 octave_idx_type n = idx_i.freeze (lhs_len, 0, true); |
| 5164 | 2855 |
| 2856 if (idx_i) | |
| 2857 { | |
| 2858 if (rhs_nr == 0 && rhs_nc == 0) | |
| 2859 { | |
| 2860 if (n != 0 && (lhs_nr != 0 || lhs_nc != 0)) | |
| 2861 lhs.maybe_delete_elements (idx_i); | |
| 2862 } | |
| 2863 else | |
| 2864 { | |
| 5781 | 2865 if (lhs_is_empty |
| 2866 && idx_i.is_colon () | |
| 2867 && ! (rhs_nr == 1 || rhs_nc == 1)) | |
| 5164 | 2868 { |
| 5781 | 2869 (*current_liboctave_warning_with_id_handler) |
| 2870 ("Octave:fortran-indexing", | |
| 2871 "A(:) = X: X is not a vector or scalar"); | |
| 2872 } | |
| 2873 else | |
| 2874 { | |
| 2875 octave_idx_type idx_nr = idx_i.orig_rows (); | |
| 2876 octave_idx_type idx_nc = idx_i.orig_columns (); | |
| 2877 | |
| 2878 if (! (rhs_nr == idx_nr && rhs_nc == idx_nc)) | |
| 2879 (*current_liboctave_warning_with_id_handler) | |
| 2880 ("Octave:fortran-indexing", | |
| 2881 "A(I) = X: X does not have same shape as I"); | |
| 5164 | 2882 } |
| 2883 | |
| 5760 | 2884 if (! assign1 (lhs, rhs)) |
| 5164 | 2885 retval = 0; |
| 2886 } | |
| 2887 } | |
| 2888 // idx_vector::freeze() printed an error message for us. | |
| 2889 } | |
| 2890 else if (lhs_nr == 1) | |
| 2891 { | |
| 5781 | 2892 idx_i.freeze (lhs_nc, "vector", true); |
| 5164 | 2893 |
| 2894 if (idx_i) | |
| 2895 { | |
| 2896 if (rhs_nr == 0 && rhs_nc == 0) | |
| 2897 lhs.maybe_delete_elements (idx_i); | |
| 5760 | 2898 else if (! assign1 (lhs, rhs)) |
| 5164 | 2899 retval = 0; |
| 2900 } | |
| 2901 // idx_vector::freeze() printed an error message for us. | |
| 2902 } | |
| 2903 else if (lhs_nc == 1) | |
| 2904 { | |
| 5781 | 2905 idx_i.freeze (lhs_nr, "vector", true); |
| 5164 | 2906 |
| 2907 if (idx_i) | |
| 2908 { | |
| 2909 if (rhs_nr == 0 && rhs_nc == 0) | |
| 2910 lhs.maybe_delete_elements (idx_i); | |
| 5760 | 2911 else if (! assign1 (lhs, rhs)) |
| 5164 | 2912 retval = 0; |
| 2913 } | |
| 2914 // idx_vector::freeze() printed an error message for us. | |
| 2915 } | |
| 2916 else | |
| 2917 { | |
| 5781 | 2918 if (! (idx_i.is_colon () |
| 2919 || (idx_i.one_zero_only () | |
| 2920 && idx_i.orig_rows () == lhs_nr | |
| 2921 && idx_i.orig_columns () == lhs_nc))) | |
| 2922 (*current_liboctave_warning_with_id_handler) | |
| 2923 ("Octave:fortran-indexing", "single index used for matrix"); | |
| 5164 | 2924 |
| 5275 | 2925 octave_idx_type lhs_len = lhs.length (); |
| 2926 | |
| 2927 octave_idx_type len = idx_i.freeze (lhs_nr * lhs_nc, "matrix"); | |
| 5164 | 2928 |
| 2929 if (idx_i) | |
| 2930 { | |
| 2931 if (rhs_nr == 0 && rhs_nc == 0) | |
| 2932 lhs.maybe_delete_elements (idx_i); | |
| 2933 else if (len == 0) | |
| 2934 { | |
| 2935 if (! ((rhs_nr == 1 && rhs_nc == 1) | |
| 2936 || (rhs_nr == 0 || rhs_nc == 0))) | |
| 2937 (*current_liboctave_error_handler) | |
| 2938 ("A([]) = X: X must be an empty matrix or scalar"); | |
| 2939 } | |
| 2940 else if (len == rhs_nr * rhs_nc) | |
| 2941 { | |
| 5604 | 2942 octave_idx_type new_nzmx = lhs_nz; |
| 5603 | 2943 OCTAVE_LOCAL_BUFFER (octave_idx_type, rhs_idx, len); |
| 2944 | |
| 2945 if (! idx_i.is_colon ()) | |
| 2946 { | |
| 2947 // Ok here we have to be careful with the indexing, to | |
| 2948 // treat cases like "a([3,2,1]) = b", and still handle | |
| 2949 // the need for strict sorting of the sparse elements. | |
| 2950 | |
| 2951 OCTAVE_LOCAL_BUFFER (octave_idx_vector_sort *, sidx, | |
| 2952 len); | |
| 2953 OCTAVE_LOCAL_BUFFER (octave_idx_vector_sort, sidxX, | |
| 2954 len); | |
| 2955 | |
| 2956 for (octave_idx_type i = 0; i < len; i++) | |
| 2957 { | |
| 2958 sidx[i] = &sidxX[i]; | |
| 2959 sidx[i]->i = idx_i.elem(i); | |
| 2960 sidx[i]->idx = i; | |
| 2961 } | |
| 2962 | |
| 2963 OCTAVE_QUIT; | |
| 2964 octave_sort<octave_idx_vector_sort *> | |
| 2965 sort (octave_idx_vector_comp); | |
| 2966 | |
| 2967 sort.sort (sidx, len); | |
| 2968 | |
| 2969 intNDArray<octave_idx_type> new_idx (dim_vector (len,1)); | |
| 2970 | |
| 2971 for (octave_idx_type i = 0; i < len; i++) | |
| 2972 { | |
| 2973 new_idx.xelem(i) = sidx[i]->i + 1; | |
| 2974 rhs_idx[i] = sidx[i]->idx; | |
| 2975 } | |
| 2976 | |
| 2977 idx_i = idx_vector (new_idx); | |
| 2978 } | |
| 2979 else | |
| 2980 for (octave_idx_type i = 0; i < len; i++) | |
| 2981 rhs_idx[i] = i; | |
| 5164 | 2982 |
| 2983 // First count the number of non-zero elements | |
| 5275 | 2984 for (octave_idx_type i = 0; i < len; i++) |
| 5164 | 2985 { |
| 2986 OCTAVE_QUIT; | |
| 2987 | |
| 5275 | 2988 octave_idx_type ii = idx_i.elem (i); |
|
11669
a6e08ecb4050
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7326
diff
changeset
|
2989 if (i < len - 1 && idx_i.elem (i + 1) == ii) |
|
a6e08ecb4050
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7326
diff
changeset
|
2990 continue; |
| 5164 | 2991 if (ii < lhs_len && c_lhs.elem(ii) != LT ()) |
| 5604 | 2992 new_nzmx--; |
| 5603 | 2993 if (rhs.elem(rhs_idx[i]) != RT ()) |
| 5604 | 2994 new_nzmx++; |
| 5164 | 2995 } |
| 2996 | |
| 5604 | 2997 Sparse<LT> stmp (lhs_nr, lhs_nc, new_nzmx); |
| 5164 | 2998 |
| 5275 | 2999 octave_idx_type i = 0; |
| 3000 octave_idx_type ii = 0; | |
| 3001 octave_idx_type ic = 0; | |
| 5164 | 3002 if (i < lhs_nz) |
| 3003 { | |
| 3004 while (ic < lhs_nc && i >= c_lhs.cidx(ic+1)) | |
| 3005 ic++; | |
| 3006 ii = ic * lhs_nr + c_lhs.ridx(i); | |
| 3007 } | |
| 3008 | |
| 5275 | 3009 octave_idx_type j = 0; |
| 3010 octave_idx_type jj = idx_i.elem (j); | |
| 3011 octave_idx_type jr = jj % lhs_nr; | |
| 3012 octave_idx_type jc = (jj - jr) / lhs_nr; | |
| 3013 | |
| 3014 octave_idx_type kk = 0; | |
| 3015 octave_idx_type kc = 0; | |
| 5164 | 3016 |
| 3017 while (j < len || i < lhs_nz) | |
| 3018 { | |
| 3019 if (j == len || (i < lhs_nz && ii < jj)) | |
| 3020 { | |
| 3021 while (kc <= ic) | |
| 3022 stmp.xcidx (kc++) = kk; | |
| 3023 stmp.xdata (kk) = c_lhs.data (i); | |
| 3024 stmp.xridx (kk++) = c_lhs.ridx (i); | |
| 3025 i++; | |
| 3026 while (ic < lhs_nc && i >= c_lhs.cidx(ic+1)) | |
| 3027 ic++; | |
| 3028 if (i < lhs_nz) | |
| 3029 ii = ic * lhs_nr + c_lhs.ridx(i); | |
| 3030 } | |
| 3031 else | |
| 3032 { | |
| 3033 while (kc <= jc) | |
| 3034 stmp.xcidx (kc++) = kk; | |
| 5603 | 3035 RT rtmp = rhs.elem (rhs_idx[j]); |
| 5164 | 3036 if (rtmp != RT ()) |
| 3037 { | |
| 3038 stmp.xdata (kk) = rtmp; | |
| 3039 stmp.xridx (kk++) = jr; | |
| 3040 } | |
| 3041 if (ii == jj) | |
| 3042 { | |
| 3043 i++; | |
| 3044 while (ic < lhs_nc && i >= c_lhs.cidx(ic+1)) | |
| 3045 ic++; | |
| 3046 if (i < lhs_nz) | |
| 3047 ii = ic * lhs_nr + c_lhs.ridx(i); | |
| 3048 } | |
| 3049 j++; | |
| 3050 if (j < len) | |
| 3051 { | |
| 3052 jj = idx_i.elem (j); | |
| 3053 jr = jj % lhs_nr; | |
| 3054 jc = (jj - jr) / lhs_nr; | |
| 3055 } | |
| 3056 } | |
| 3057 } | |
| 3058 | |
| 5275 | 3059 for (octave_idx_type iidx = kc; iidx < lhs_nc+1; iidx++) |
| 5603 | 3060 stmp.xcidx(iidx) = kk; |
| 5164 | 3061 |
| 3062 lhs = stmp; | |
| 3063 } | |
| 3064 else if (rhs_nr == 1 && rhs_nc == 1) | |
| 3065 { | |
| 3066 RT scalar = rhs.elem (0, 0); | |
| 5604 | 3067 octave_idx_type new_nzmx = lhs_nz; |
| 5603 | 3068 idx_i.sort (true); |
|
11669
a6e08ecb4050
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7326
diff
changeset
|
3069 len = idx_i.length (len); |
| 5164 | 3070 |
| 3071 // First count the number of non-zero elements | |
| 3072 if (scalar != RT ()) | |
| 5604 | 3073 new_nzmx += len; |
| 5275 | 3074 for (octave_idx_type i = 0; i < len; i++) |
| 5164 | 3075 { |
| 3076 OCTAVE_QUIT; | |
| 5275 | 3077 octave_idx_type ii = idx_i.elem (i); |
| 5164 | 3078 if (ii < lhs_len && c_lhs.elem(ii) != LT ()) |
| 5604 | 3079 new_nzmx--; |
| 5164 | 3080 } |
| 3081 | |
| 5604 | 3082 Sparse<LT> stmp (lhs_nr, lhs_nc, new_nzmx); |
| 5164 | 3083 |
| 5275 | 3084 octave_idx_type i = 0; |
| 3085 octave_idx_type ii = 0; | |
| 3086 octave_idx_type ic = 0; | |
| 5164 | 3087 if (i < lhs_nz) |
| 3088 { | |
| 3089 while (ic < lhs_nc && i >= c_lhs.cidx(ic+1)) | |
| 3090 ic++; | |
| 3091 ii = ic * lhs_nr + c_lhs.ridx(i); | |
| 3092 } | |
| 3093 | |
| 5275 | 3094 octave_idx_type j = 0; |
| 3095 octave_idx_type jj = idx_i.elem (j); | |
| 3096 octave_idx_type jr = jj % lhs_nr; | |
| 3097 octave_idx_type jc = (jj - jr) / lhs_nr; | |
| 3098 | |
| 3099 octave_idx_type kk = 0; | |
| 3100 octave_idx_type kc = 0; | |
| 5164 | 3101 |
| 3102 while (j < len || i < lhs_nz) | |
| 3103 { | |
|
11669
a6e08ecb4050
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7326
diff
changeset
|
3104 if (j < len - 1 && idx_i.elem (j + 1) == jj) |
|
a6e08ecb4050
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7326
diff
changeset
|
3105 { |
|
a6e08ecb4050
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7326
diff
changeset
|
3106 j++; |
|
a6e08ecb4050
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7326
diff
changeset
|
3107 jj = idx_i.elem (j); |
|
a6e08ecb4050
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7326
diff
changeset
|
3108 jr = jj % lhs_nr; |
|
a6e08ecb4050
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7326
diff
changeset
|
3109 jc = (jj - jr) / lhs_nr; |
|
a6e08ecb4050
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7326
diff
changeset
|
3110 continue; |
|
a6e08ecb4050
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7326
diff
changeset
|
3111 } |
|
a6e08ecb4050
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7326
diff
changeset
|
3112 |
| 5164 | 3113 if (j == len || (i < lhs_nz && ii < jj)) |
| 3114 { | |
| 3115 while (kc <= ic) | |
| 3116 stmp.xcidx (kc++) = kk; | |
| 3117 stmp.xdata (kk) = c_lhs.data (i); | |
| 3118 stmp.xridx (kk++) = c_lhs.ridx (i); | |
| 3119 i++; | |
| 3120 while (ic < lhs_nc && i >= c_lhs.cidx(ic+1)) | |
| 3121 ic++; | |
| 3122 if (i < lhs_nz) | |
| 3123 ii = ic * lhs_nr + c_lhs.ridx(i); | |
| 3124 } | |
| 3125 else | |
| 3126 { | |
| 3127 while (kc <= jc) | |
| 3128 stmp.xcidx (kc++) = kk; | |
| 3129 if (scalar != RT ()) | |
| 3130 { | |
| 3131 stmp.xdata (kk) = scalar; | |
| 3132 stmp.xridx (kk++) = jr; | |
| 3133 } | |
| 3134 if (ii == jj) | |
| 3135 { | |
| 3136 i++; | |
| 3137 while (ic < lhs_nc && i >= c_lhs.cidx(ic+1)) | |
| 3138 ic++; | |
| 3139 if (i < lhs_nz) | |
| 3140 ii = ic * lhs_nr + c_lhs.ridx(i); | |
| 3141 } | |
| 3142 j++; | |
| 3143 if (j < len) | |
| 3144 { | |
| 3145 jj = idx_i.elem (j); | |
| 3146 jr = jj % lhs_nr; | |
| 3147 jc = (jj - jr) / lhs_nr; | |
| 3148 } | |
| 3149 } | |
| 3150 } | |
| 3151 | |
| 5275 | 3152 for (octave_idx_type iidx = kc; iidx < lhs_nc+1; iidx++) |
| 5164 | 3153 stmp.xcidx(iidx) = kk; |
| 3154 | |
| 3155 lhs = stmp; | |
| 3156 } | |
| 3157 else | |
| 3158 { | |
| 3159 (*current_liboctave_error_handler) | |
| 3160 ("A(I) = X: X must be a scalar or a matrix with the same size as I"); | |
| 3161 | |
| 3162 retval = 0; | |
| 3163 } | |
| 3164 } | |
| 3165 // idx_vector::freeze() printed an error message for us. | |
| 3166 } | |
| 3167 } | |
| 3168 else | |
| 3169 { | |
| 3170 (*current_liboctave_error_handler) | |
| 3171 ("invalid number of indices for matrix expression"); | |
| 3172 | |
| 3173 retval = 0; | |
| 3174 } | |
| 3175 | |
| 3176 lhs.clear_index (); | |
| 3177 | |
| 3178 return retval; | |
| 3179 } | |
| 3180 | |
| 3181 template <class T> | |
| 3182 void | |
| 3183 Sparse<T>::print_info (std::ostream& os, const std::string& prefix) const | |
| 3184 { | |
| 3185 os << prefix << "rep address: " << rep << "\n" | |
| 5604 | 3186 << prefix << "rep->nzmx: " << rep->nzmx << "\n" |
| 5164 | 3187 << prefix << "rep->nrows: " << rep->nrows << "\n" |
| 3188 << prefix << "rep->ncols: " << rep->ncols << "\n" | |
| 3189 << prefix << "rep->data: " << static_cast<void *> (rep->d) << "\n" | |
| 3190 << prefix << "rep->ridx: " << static_cast<void *> (rep->r) << "\n" | |
| 3191 << prefix << "rep->cidx: " << static_cast<void *> (rep->c) << "\n" | |
| 3192 << prefix << "rep->count: " << rep->count << "\n"; | |
| 3193 } | |
| 3194 | |
| 3195 /* | |
| 3196 ;;; Local Variables: *** | |
| 3197 ;;; mode: C++ *** | |
| 3198 ;;; End: *** | |
| 3199 */ |