Mercurial > hg > octave-nkf
annotate liboctave/Sparse.cc @ 11951:9cfbc1a1bf0b release-3-0-x
this branch is no longer maintained and is closed for further development
| author | John W. Eaton <jwe@octave.org> |
|---|---|
| date | Sat, 22 Jan 2011 00:59:43 -0500 |
| parents | ccab9d3d1d21 |
| children |
| rev | line source |
|---|---|
| 5164 | 1 // Template sparse array class |
| 2 /* | |
| 3 | |
| 11740 | 4 Copyright (C) 2004, 2005, 2006, 2007, 2008 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> | |
|
11755
ccab9d3d1d21
Delete idx in Sparse<T> and Array<T> operator =
David Bateman <dbateman@free.fr>
parents:
11740
diff
changeset
|
618 Sparse<T>& |
|
ccab9d3d1d21
Delete idx in Sparse<T> and Array<T> operator =
David Bateman <dbateman@free.fr>
parents:
11740
diff
changeset
|
619 Sparse<T>::operator = (const Sparse<T>& a) |
|
ccab9d3d1d21
Delete idx in Sparse<T> and Array<T> operator =
David Bateman <dbateman@free.fr>
parents:
11740
diff
changeset
|
620 { |
|
ccab9d3d1d21
Delete idx in Sparse<T> and Array<T> operator =
David Bateman <dbateman@free.fr>
parents:
11740
diff
changeset
|
621 if (this != &a) |
|
ccab9d3d1d21
Delete idx in Sparse<T> and Array<T> operator =
David Bateman <dbateman@free.fr>
parents:
11740
diff
changeset
|
622 { |
|
ccab9d3d1d21
Delete idx in Sparse<T> and Array<T> operator =
David Bateman <dbateman@free.fr>
parents:
11740
diff
changeset
|
623 if (--rep->count <= 0) |
|
ccab9d3d1d21
Delete idx in Sparse<T> and Array<T> operator =
David Bateman <dbateman@free.fr>
parents:
11740
diff
changeset
|
624 delete rep; |
|
ccab9d3d1d21
Delete idx in Sparse<T> and Array<T> operator =
David Bateman <dbateman@free.fr>
parents:
11740
diff
changeset
|
625 |
|
ccab9d3d1d21
Delete idx in Sparse<T> and Array<T> operator =
David Bateman <dbateman@free.fr>
parents:
11740
diff
changeset
|
626 rep = a.rep; |
|
ccab9d3d1d21
Delete idx in Sparse<T> and Array<T> operator =
David Bateman <dbateman@free.fr>
parents:
11740
diff
changeset
|
627 rep->count++; |
|
ccab9d3d1d21
Delete idx in Sparse<T> and Array<T> operator =
David Bateman <dbateman@free.fr>
parents:
11740
diff
changeset
|
628 |
|
ccab9d3d1d21
Delete idx in Sparse<T> and Array<T> operator =
David Bateman <dbateman@free.fr>
parents:
11740
diff
changeset
|
629 dimensions = a.dimensions; |
|
ccab9d3d1d21
Delete idx in Sparse<T> and Array<T> operator =
David Bateman <dbateman@free.fr>
parents:
11740
diff
changeset
|
630 |
|
ccab9d3d1d21
Delete idx in Sparse<T> and Array<T> operator =
David Bateman <dbateman@free.fr>
parents:
11740
diff
changeset
|
631 delete [] idx; |
|
ccab9d3d1d21
Delete idx in Sparse<T> and Array<T> operator =
David Bateman <dbateman@free.fr>
parents:
11740
diff
changeset
|
632 idx_count = 0; |
|
ccab9d3d1d21
Delete idx in Sparse<T> and Array<T> operator =
David Bateman <dbateman@free.fr>
parents:
11740
diff
changeset
|
633 idx = 0; |
|
ccab9d3d1d21
Delete idx in Sparse<T> and Array<T> operator =
David Bateman <dbateman@free.fr>
parents:
11740
diff
changeset
|
634 } |
|
ccab9d3d1d21
Delete idx in Sparse<T> and Array<T> operator =
David Bateman <dbateman@free.fr>
parents:
11740
diff
changeset
|
635 |
|
ccab9d3d1d21
Delete idx in Sparse<T> and Array<T> operator =
David Bateman <dbateman@free.fr>
parents:
11740
diff
changeset
|
636 return *this; |
|
ccab9d3d1d21
Delete idx in Sparse<T> and Array<T> operator =
David Bateman <dbateman@free.fr>
parents:
11740
diff
changeset
|
637 } |
|
ccab9d3d1d21
Delete idx in Sparse<T> and Array<T> operator =
David Bateman <dbateman@free.fr>
parents:
11740
diff
changeset
|
638 |
|
ccab9d3d1d21
Delete idx in Sparse<T> and Array<T> operator =
David Bateman <dbateman@free.fr>
parents:
11740
diff
changeset
|
639 template <class T> |
| 5275 | 640 octave_idx_type |
| 641 Sparse<T>::compute_index (const Array<octave_idx_type>& ra_idx) const | |
| 5164 | 642 { |
| 5275 | 643 octave_idx_type retval = -1; |
| 644 | |
| 645 octave_idx_type n = dimensions.length (); | |
| 5164 | 646 |
| 647 if (n > 0 && n == ra_idx.length ()) | |
| 648 { | |
| 649 retval = ra_idx(--n); | |
| 650 | |
| 651 while (--n >= 0) | |
| 652 { | |
| 653 retval *= dimensions(n); | |
| 654 retval += ra_idx(n); | |
| 655 } | |
| 656 } | |
| 657 else | |
| 658 (*current_liboctave_error_handler) | |
| 659 ("Sparse<T>::compute_index: invalid ra_idxing operation"); | |
| 660 | |
| 661 return retval; | |
| 662 } | |
| 663 | |
| 664 template <class T> | |
| 665 T | |
| 5275 | 666 Sparse<T>::range_error (const char *fcn, octave_idx_type n) const |
| 5164 | 667 { |
| 668 (*current_liboctave_error_handler) ("%s (%d): range error", fcn, n); | |
| 669 return T (); | |
| 670 } | |
| 671 | |
| 672 template <class T> | |
| 673 T& | |
| 5275 | 674 Sparse<T>::range_error (const char *fcn, octave_idx_type n) |
| 5164 | 675 { |
| 676 (*current_liboctave_error_handler) ("%s (%d): range error", fcn, n); | |
| 677 static T foo; | |
| 678 return foo; | |
| 679 } | |
| 680 | |
| 681 template <class T> | |
| 682 T | |
| 5275 | 683 Sparse<T>::range_error (const char *fcn, octave_idx_type i, octave_idx_type j) const |
| 5164 | 684 { |
| 685 (*current_liboctave_error_handler) | |
| 686 ("%s (%d, %d): range error", fcn, i, j); | |
| 687 return T (); | |
| 688 } | |
| 689 | |
| 690 template <class T> | |
| 691 T& | |
| 5275 | 692 Sparse<T>::range_error (const char *fcn, octave_idx_type i, octave_idx_type j) |
| 5164 | 693 { |
| 694 (*current_liboctave_error_handler) | |
| 695 ("%s (%d, %d): range error", fcn, i, j); | |
| 696 static T foo; | |
| 697 return foo; | |
| 698 } | |
| 699 | |
| 700 template <class T> | |
| 701 T | |
| 5275 | 702 Sparse<T>::range_error (const char *fcn, const Array<octave_idx_type>& ra_idx) const |
| 5164 | 703 { |
| 5765 | 704 std::ostringstream buf; |
| 5164 | 705 |
| 706 buf << fcn << " ("; | |
| 707 | |
| 5275 | 708 octave_idx_type n = ra_idx.length (); |
| 5164 | 709 |
| 710 if (n > 0) | |
| 711 buf << ra_idx(0); | |
| 712 | |
| 5275 | 713 for (octave_idx_type i = 1; i < n; i++) |
| 5164 | 714 buf << ", " << ra_idx(i); |
| 715 | |
| 716 buf << "): range error"; | |
| 5765 | 717 |
| 718 std::string buf_str = buf.str (); | |
| 719 | |
| 720 (*current_liboctave_error_handler) (buf_str.c_str ()); | |
| 5164 | 721 |
| 722 return T (); | |
| 723 } | |
| 724 | |
| 725 template <class T> | |
| 726 T& | |
| 5275 | 727 Sparse<T>::range_error (const char *fcn, const Array<octave_idx_type>& ra_idx) |
| 5164 | 728 { |
| 5765 | 729 std::ostringstream buf; |
| 5164 | 730 |
| 731 buf << fcn << " ("; | |
| 732 | |
| 5275 | 733 octave_idx_type n = ra_idx.length (); |
| 5164 | 734 |
| 735 if (n > 0) | |
| 736 buf << ra_idx(0); | |
| 737 | |
| 5275 | 738 for (octave_idx_type i = 1; i < n; i++) |
| 5164 | 739 buf << ", " << ra_idx(i); |
| 740 | |
| 741 buf << "): range error"; | |
| 742 | |
| 5765 | 743 std::string buf_str = buf.str (); |
| 744 | |
| 745 (*current_liboctave_error_handler) (buf_str.c_str ()); | |
| 5164 | 746 |
| 747 static T foo; | |
| 748 return foo; | |
| 749 } | |
| 750 | |
| 751 template <class T> | |
| 752 Sparse<T> | |
| 753 Sparse<T>::reshape (const dim_vector& new_dims) const | |
| 754 { | |
| 755 Sparse<T> retval; | |
| 6689 | 756 dim_vector dims2 = new_dims; |
| 757 | |
| 758 if (dims2.length () > 2) | |
| 5164 | 759 { |
| 6814 | 760 (*current_liboctave_warning_handler) |
| 761 ("reshape: sparse reshape to N-d array smashes dims"); | |
| 762 | |
| 6689 | 763 for (octave_idx_type i = 2; i < dims2.length(); i++) |
| 6814 | 764 dims2(1) *= dims2(i); |
| 765 | |
| 6689 | 766 dims2.resize (2); |
| 767 } | |
| 768 | |
| 769 if (dimensions != dims2) | |
| 770 { | |
| 771 if (dimensions.numel () == dims2.numel ()) | |
| 5164 | 772 { |
| 5681 | 773 octave_idx_type new_nnz = nnz (); |
| 6689 | 774 octave_idx_type new_nr = dims2 (0); |
| 775 octave_idx_type new_nc = dims2 (1); | |
| 5275 | 776 octave_idx_type old_nr = rows (); |
| 777 octave_idx_type old_nc = cols (); | |
| 5681 | 778 retval = Sparse<T> (new_nr, new_nc, new_nnz); |
| 5164 | 779 |
| 5275 | 780 octave_idx_type kk = 0; |
| 5164 | 781 retval.xcidx(0) = 0; |
| 5275 | 782 for (octave_idx_type i = 0; i < old_nc; i++) |
| 783 for (octave_idx_type j = cidx(i); j < cidx(i+1); j++) | |
| 5164 | 784 { |
| 5275 | 785 octave_idx_type tmp = i * old_nr + ridx(j); |
| 786 octave_idx_type ii = tmp % new_nr; | |
| 787 octave_idx_type jj = (tmp - ii) / new_nr; | |
| 788 for (octave_idx_type k = kk; k < jj; k++) | |
| 5164 | 789 retval.xcidx(k+1) = j; |
| 790 kk = jj; | |
| 791 retval.xdata(j) = data(j); | |
| 792 retval.xridx(j) = ii; | |
| 793 } | |
| 5275 | 794 for (octave_idx_type k = kk; k < new_nc; k++) |
| 5681 | 795 retval.xcidx(k+1) = new_nnz; |
| 5164 | 796 } |
| 797 else | |
| 798 (*current_liboctave_error_handler) ("reshape: size mismatch"); | |
| 799 } | |
| 800 else | |
| 801 retval = *this; | |
| 802 | |
| 803 return retval; | |
| 804 } | |
| 805 | |
| 806 template <class T> | |
| 807 Sparse<T> | |
| 5275 | 808 Sparse<T>::permute (const Array<octave_idx_type>& perm_vec, bool) const |
| 5164 | 809 { |
| 6813 | 810 // The only valid permutations of a sparse array are [1, 2] and [2, 1]. |
| 811 | |
| 812 bool fail = false; | |
| 6817 | 813 bool trans = false; |
| 6813 | 814 |
| 815 if (perm_vec.length () == 2) | |
| 5164 | 816 { |
| 6813 | 817 if (perm_vec(0) == 0 && perm_vec(1) == 1) |
| 818 /* do nothing */; | |
| 819 else if (perm_vec(0) == 1 && perm_vec(1) == 0) | |
| 6817 | 820 trans = true; |
| 5164 | 821 else |
| 6813 | 822 fail = true; |
| 5164 | 823 } |
| 824 else | |
| 6813 | 825 fail = true; |
| 826 | |
| 827 if (fail) | |
| 828 (*current_liboctave_error_handler) | |
| 829 ("permutation vector contains an invalid element"); | |
| 830 | |
| 6817 | 831 return trans ? this->transpose () : *this; |
| 5164 | 832 } |
| 833 | |
| 834 template <class T> | |
| 835 void | |
| 836 Sparse<T>::resize_no_fill (const dim_vector& dv) | |
| 837 { | |
| 5275 | 838 octave_idx_type n = dv.length (); |
| 5164 | 839 |
| 840 if (n != 2) | |
| 841 { | |
| 842 (*current_liboctave_error_handler) ("sparse array must be 2-D"); | |
| 843 return; | |
| 844 } | |
| 845 | |
| 846 resize_no_fill (dv(0), dv(1)); | |
| 847 } | |
| 848 | |
| 849 template <class T> | |
| 850 void | |
| 5275 | 851 Sparse<T>::resize_no_fill (octave_idx_type r, octave_idx_type c) |
| 5164 | 852 { |
| 853 if (r < 0 || c < 0) | |
| 854 { | |
| 855 (*current_liboctave_error_handler) | |
| 856 ("can't resize to negative dimension"); | |
| 857 return; | |
| 858 } | |
| 859 | |
| 860 if (ndims () == 0) | |
| 861 dimensions = dim_vector (0, 0); | |
| 862 | |
| 863 if (r == dim1 () && c == dim2 ()) | |
| 864 return; | |
| 865 | |
| 5731 | 866 typename Sparse<T>::SparseRep *old_rep = rep; |
| 867 | |
| 5275 | 868 octave_idx_type nc = cols (); |
| 869 octave_idx_type nr = rows (); | |
| 5164 | 870 |
| 5681 | 871 if (nnz () == 0 || r == 0 || c == 0) |
| 5164 | 872 // Special case of redimensioning to/from a sparse matrix with |
| 873 // no elements | |
| 874 rep = new typename Sparse<T>::SparseRep (r, c); | |
| 875 else | |
| 876 { | |
| 5275 | 877 octave_idx_type n = 0; |
| 5164 | 878 Sparse<T> tmpval; |
| 879 if (r >= nr) | |
| 880 { | |
| 881 if (c > nc) | |
| 5731 | 882 n = xcidx(nc); |
| 5164 | 883 else |
| 5731 | 884 n = xcidx(c); |
| 5164 | 885 |
| 886 tmpval = Sparse<T> (r, c, n); | |
| 887 | |
| 888 if (c > nc) | |
| 889 { | |
| 6101 | 890 for (octave_idx_type i = 0; i < nc + 1; i++) |
| 5731 | 891 tmpval.cidx(i) = xcidx(i); |
| 6101 | 892 for (octave_idx_type i = nc + 1; i < c + 1; i++) |
| 5164 | 893 tmpval.cidx(i) = tmpval.cidx(i-1); |
| 894 } | |
| 895 else if (c <= nc) | |
| 6101 | 896 for (octave_idx_type i = 0; i < c + 1; i++) |
| 5731 | 897 tmpval.cidx(i) = xcidx(i); |
| 5164 | 898 |
| 5275 | 899 for (octave_idx_type i = 0; i < n; i++) |
| 5164 | 900 { |
| 5731 | 901 tmpval.data(i) = xdata(i); |
| 902 tmpval.ridx(i) = xridx(i); | |
| 5164 | 903 } |
| 904 } | |
| 905 else | |
| 906 { | |
| 907 // Count how many non zero terms before we do anything | |
| 6101 | 908 octave_idx_type min_nc = (c < nc ? c : nc); |
| 909 for (octave_idx_type i = 0; i < min_nc; i++) | |
| 5731 | 910 for (octave_idx_type j = xcidx(i); j < xcidx(i+1); j++) |
| 911 if (xridx(j) < r) | |
| 5164 | 912 n++; |
| 913 | |
| 914 if (n) | |
| 915 { | |
| 916 // Now that we know the size we can do something | |
| 917 tmpval = Sparse<T> (r, c, n); | |
| 918 | |
| 919 tmpval.cidx(0); | |
| 6101 | 920 for (octave_idx_type i = 0, ii = 0; i < min_nc; i++) |
| 5164 | 921 { |
| 5731 | 922 for (octave_idx_type j = xcidx(i); j < xcidx(i+1); j++) |
| 923 if (xridx(j) < r) | |
| 5164 | 924 { |
| 5731 | 925 tmpval.data(ii) = xdata(j); |
| 926 tmpval.ridx(ii++) = xridx(j); | |
| 5164 | 927 } |
| 928 tmpval.cidx(i+1) = ii; | |
| 929 } | |
| 6101 | 930 if (c > min_nc) |
| 931 for (octave_idx_type i = nc; i < c; i++) | |
| 932 tmpval.cidx(i+1) = tmpval.cidx(i); | |
| 5164 | 933 } |
| 934 else | |
| 935 tmpval = Sparse<T> (r, c); | |
| 936 } | |
| 937 | |
| 938 rep = tmpval.rep; | |
| 939 rep->count++; | |
| 940 } | |
| 941 | |
| 942 dimensions = dim_vector (r, c); | |
| 943 | |
| 944 if (--old_rep->count <= 0) | |
| 945 delete old_rep; | |
| 946 } | |
| 947 | |
| 948 template <class T> | |
| 949 Sparse<T>& | |
| 5275 | 950 Sparse<T>::insert (const Sparse<T>& a, octave_idx_type r, octave_idx_type c) |
| 5164 | 951 { |
| 5275 | 952 octave_idx_type a_rows = a.rows (); |
| 953 octave_idx_type a_cols = a.cols (); | |
| 954 octave_idx_type nr = rows (); | |
| 955 octave_idx_type nc = cols (); | |
| 5164 | 956 |
| 957 if (r < 0 || r + a_rows > rows () || c < 0 || c + a_cols > cols ()) | |
| 958 { | |
| 959 (*current_liboctave_error_handler) ("range error for insert"); | |
| 960 return *this; | |
| 961 } | |
| 962 | |
| 963 // First count the number of elements in the final array | |
| 5681 | 964 octave_idx_type nel = cidx(c) + a.nnz (); |
| 5164 | 965 |
| 966 if (c + a_cols < nc) | |
| 967 nel += cidx(nc) - cidx(c + a_cols); | |
| 968 | |
| 5275 | 969 for (octave_idx_type i = c; i < c + a_cols; i++) |
| 970 for (octave_idx_type j = cidx(i); j < cidx(i+1); j++) | |
| 5164 | 971 if (ridx(j) < r || ridx(j) >= r + a_rows) |
| 972 nel++; | |
| 973 | |
| 974 Sparse<T> tmp (*this); | |
| 975 --rep->count; | |
| 976 rep = new typename Sparse<T>::SparseRep (nr, nc, nel); | |
| 977 | |
| 5275 | 978 for (octave_idx_type i = 0; i < tmp.cidx(c); i++) |
| 5164 | 979 { |
| 980 data(i) = tmp.data(i); | |
| 981 ridx(i) = tmp.ridx(i); | |
| 982 } | |
| 5275 | 983 for (octave_idx_type i = 0; i < c + 1; i++) |
| 5164 | 984 cidx(i) = tmp.cidx(i); |
| 985 | |
| 5275 | 986 octave_idx_type ii = cidx(c); |
| 987 | |
| 988 for (octave_idx_type i = c; i < c + a_cols; i++) | |
| 5164 | 989 { |
| 990 OCTAVE_QUIT; | |
| 991 | |
| 5275 | 992 for (octave_idx_type j = tmp.cidx(i); j < tmp.cidx(i+1); j++) |
| 5164 | 993 if (tmp.ridx(j) < r) |
| 994 { | |
| 995 data(ii) = tmp.data(j); | |
| 996 ridx(ii++) = tmp.ridx(j); | |
| 997 } | |
| 998 | |
| 999 OCTAVE_QUIT; | |
| 1000 | |
| 5275 | 1001 for (octave_idx_type j = a.cidx(i-c); j < a.cidx(i-c+1); j++) |
| 5164 | 1002 { |
| 1003 data(ii) = a.data(j); | |
| 1004 ridx(ii++) = r + a.ridx(j); | |
| 1005 } | |
| 1006 | |
| 1007 OCTAVE_QUIT; | |
| 1008 | |
| 5275 | 1009 for (octave_idx_type j = tmp.cidx(i); j < tmp.cidx(i+1); j++) |
| 5164 | 1010 if (tmp.ridx(j) >= r + a_rows) |
| 1011 { | |
| 1012 data(ii) = tmp.data(j); | |
| 1013 ridx(ii++) = tmp.ridx(j); | |
| 1014 } | |
| 1015 | |
| 1016 cidx(i+1) = ii; | |
| 1017 } | |
| 1018 | |
| 5275 | 1019 for (octave_idx_type i = c + a_cols; i < nc; i++) |
| 5164 | 1020 { |
| 5275 | 1021 for (octave_idx_type j = tmp.cidx(i); j < tmp.cidx(i+1); j++) |
| 5164 | 1022 { |
| 1023 data(ii) = tmp.data(j); | |
| 1024 ridx(ii++) = tmp.ridx(j); | |
| 1025 } | |
| 1026 cidx(i+1) = ii; | |
| 1027 } | |
| 1028 | |
| 1029 return *this; | |
| 1030 } | |
| 1031 | |
| 1032 template <class T> | |
| 1033 Sparse<T>& | |
| 5275 | 1034 Sparse<T>::insert (const Sparse<T>& a, const Array<octave_idx_type>& ra_idx) |
| 5164 | 1035 { |
| 1036 | |
| 1037 if (ra_idx.length () != 2) | |
| 1038 { | |
| 1039 (*current_liboctave_error_handler) ("range error for insert"); | |
| 1040 return *this; | |
| 1041 } | |
| 1042 | |
| 1043 return insert (a, ra_idx (0), ra_idx (1)); | |
| 1044 } | |
| 1045 | |
| 1046 template <class T> | |
| 1047 Sparse<T> | |
| 1048 Sparse<T>::transpose (void) const | |
| 1049 { | |
| 1050 assert (ndims () == 2); | |
| 1051 | |
| 5275 | 1052 octave_idx_type nr = rows (); |
| 1053 octave_idx_type nc = cols (); | |
| 5648 | 1054 octave_idx_type nz = nnz (); |
| 5164 | 1055 Sparse<T> retval (nc, nr, nz); |
| 1056 | |
| 5648 | 1057 OCTAVE_LOCAL_BUFFER (octave_idx_type, w, nr + 1); |
| 1058 for (octave_idx_type i = 0; i < nr; i++) | |
| 1059 w[i] = 0; | |
| 1060 for (octave_idx_type i = 0; i < nz; i++) | |
| 1061 w[ridx(i)]++; | |
| 1062 nz = 0; | |
| 1063 for (octave_idx_type i = 0; i < nr; i++) | |
| 5164 | 1064 { |
| 5648 | 1065 retval.xcidx(i) = nz; |
| 1066 nz += w[i]; | |
| 1067 w[i] = retval.xcidx(i); | |
| 5164 | 1068 } |
| 5648 | 1069 retval.xcidx(nr) = nz; |
| 1070 w[nr] = nz; | |
| 1071 | |
| 1072 for (octave_idx_type j = 0; j < nc; j++) | |
| 1073 for (octave_idx_type k = cidx(j); k < cidx(j+1); k++) | |
| 1074 { | |
| 1075 octave_idx_type q = w [ridx(k)]++; | |
| 1076 retval.xridx (q) = j; | |
| 1077 retval.xdata (q) = data (k); | |
| 1078 } | |
| 5164 | 1079 |
| 1080 return retval; | |
| 1081 } | |
| 1082 | |
| 1083 template <class T> | |
| 1084 void | |
| 1085 Sparse<T>::clear_index (void) | |
| 1086 { | |
| 1087 delete [] idx; | |
| 1088 idx = 0; | |
| 1089 idx_count = 0; | |
| 1090 } | |
| 1091 | |
| 1092 template <class T> | |
| 1093 void | |
| 1094 Sparse<T>::set_index (const idx_vector& idx_arg) | |
| 1095 { | |
| 5275 | 1096 octave_idx_type nd = ndims (); |
| 5164 | 1097 |
| 1098 if (! idx && nd > 0) | |
| 1099 idx = new idx_vector [nd]; | |
| 1100 | |
| 1101 if (idx_count < nd) | |
| 1102 { | |
| 1103 idx[idx_count++] = idx_arg; | |
| 1104 } | |
| 1105 else | |
| 1106 { | |
| 1107 idx_vector *new_idx = new idx_vector [idx_count+1]; | |
| 1108 | |
| 5275 | 1109 for (octave_idx_type i = 0; i < idx_count; i++) |
| 5164 | 1110 new_idx[i] = idx[i]; |
| 1111 | |
| 1112 new_idx[idx_count++] = idx_arg; | |
| 1113 | |
| 1114 delete [] idx; | |
| 1115 | |
| 1116 idx = new_idx; | |
| 1117 } | |
| 1118 } | |
| 1119 | |
| 1120 template <class T> | |
| 1121 void | |
| 1122 Sparse<T>::maybe_delete_elements (idx_vector& idx_arg) | |
| 1123 { | |
| 5275 | 1124 octave_idx_type nr = dim1 (); |
| 1125 octave_idx_type nc = dim2 (); | |
| 5164 | 1126 |
| 1127 if (nr == 0 && nc == 0) | |
| 1128 return; | |
| 1129 | |
| 5275 | 1130 octave_idx_type n; |
| 5164 | 1131 if (nr == 1) |
| 1132 n = nc; | |
| 1133 else if (nc == 1) | |
| 1134 n = nr; | |
| 1135 else | |
| 1136 { | |
| 1137 // Reshape to row vector for Matlab compatibility. | |
| 1138 | |
| 1139 n = nr * nc; | |
| 1140 nr = 1; | |
| 1141 nc = n; | |
| 1142 } | |
| 1143 | |
| 1144 if (idx_arg.is_colon_equiv (n, 1)) | |
| 1145 { | |
| 1146 // Either A(:) = [] or A(idx) = [] with idx enumerating all | |
| 1147 // elements, so we delete all elements and return [](0x0). To | |
| 1148 // preserve the orientation of the vector, you have to use | |
| 1149 // A(idx,:) = [] (delete rows) or A(:,idx) (delete columns). | |
| 1150 | |
| 1151 resize_no_fill (0, 0); | |
| 1152 return; | |
| 1153 } | |
| 1154 | |
| 1155 idx_arg.sort (true); | |
| 1156 | |
| 5275 | 1157 octave_idx_type num_to_delete = idx_arg.length (n); |
| 5164 | 1158 |
| 1159 if (num_to_delete != 0) | |
| 1160 { | |
| 5275 | 1161 octave_idx_type new_n = n; |
| 5681 | 1162 octave_idx_type new_nnz = nnz (); |
| 5275 | 1163 |
| 1164 octave_idx_type iidx = 0; | |
| 5164 | 1165 |
| 1166 const Sparse<T> tmp (*this); | |
| 1167 | |
| 5275 | 1168 for (octave_idx_type i = 0; i < n; i++) |
| 5164 | 1169 { |
| 1170 OCTAVE_QUIT; | |
| 1171 | |
| 1172 if (i == idx_arg.elem (iidx)) | |
| 1173 { | |
| 1174 iidx++; | |
| 1175 new_n--; | |
| 1176 | |
| 1177 if (tmp.elem (i) != T ()) | |
| 5681 | 1178 new_nnz--; |
| 5164 | 1179 |
| 1180 if (iidx == num_to_delete) | |
| 1181 break; | |
| 1182 } | |
| 1183 } | |
| 1184 | |
| 1185 if (new_n > 0) | |
| 1186 { | |
| 1187 rep->count--; | |
| 1188 | |
| 1189 if (nr == 1) | |
| 5681 | 1190 rep = new typename Sparse<T>::SparseRep (1, new_n, new_nnz); |
| 5164 | 1191 else |
| 5681 | 1192 rep = new typename Sparse<T>::SparseRep (new_n, 1, new_nnz); |
| 5164 | 1193 |
| 5275 | 1194 octave_idx_type ii = 0; |
| 1195 octave_idx_type jj = 0; | |
| 5164 | 1196 iidx = 0; |
| 5275 | 1197 for (octave_idx_type i = 0; i < n; i++) |
| 5164 | 1198 { |
| 1199 OCTAVE_QUIT; | |
| 1200 | |
| 1201 if (iidx < num_to_delete && i == idx_arg.elem (iidx)) | |
| 1202 iidx++; | |
| 1203 else | |
| 1204 { | |
| 1205 T el = tmp.elem (i); | |
| 1206 if (el != T ()) | |
| 1207 { | |
| 1208 data(ii) = el; | |
| 1209 ridx(ii++) = jj; | |
| 1210 } | |
| 1211 jj++; | |
| 1212 } | |
| 1213 } | |
| 1214 | |
| 1215 dimensions.resize (2); | |
| 1216 | |
| 1217 if (nr == 1) | |
| 1218 { | |
| 1219 ii = 0; | |
| 1220 cidx(0) = 0; | |
| 5275 | 1221 for (octave_idx_type i = 0; i < new_n; i++) |
| 5164 | 1222 { |
| 1223 OCTAVE_QUIT; | |
| 1224 if (ridx(ii) == i) | |
| 1225 ridx(ii++) = 0; | |
| 1226 cidx(i+1) = ii; | |
| 1227 } | |
| 1228 | |
| 1229 dimensions(0) = 1; | |
| 1230 dimensions(1) = new_n; | |
| 1231 } | |
| 1232 else | |
| 1233 { | |
| 1234 cidx(0) = 0; | |
| 5681 | 1235 cidx(1) = new_nnz; |
| 5164 | 1236 dimensions(0) = new_n; |
| 1237 dimensions(1) = 1; | |
| 1238 } | |
| 1239 } | |
| 1240 else | |
| 1241 (*current_liboctave_error_handler) | |
| 1242 ("A(idx) = []: index out of range"); | |
| 1243 } | |
| 1244 } | |
| 1245 | |
| 1246 template <class T> | |
| 1247 void | |
| 1248 Sparse<T>::maybe_delete_elements (idx_vector& idx_i, idx_vector& idx_j) | |
| 1249 { | |
| 1250 assert (ndims () == 2); | |
| 1251 | |
| 5275 | 1252 octave_idx_type nr = dim1 (); |
| 1253 octave_idx_type nc = dim2 (); | |
| 5164 | 1254 |
| 1255 if (nr == 0 && nc == 0) | |
| 1256 return; | |
| 1257 | |
| 1258 if (idx_i.is_colon ()) | |
| 1259 { | |
| 1260 if (idx_j.is_colon ()) | |
| 1261 { | |
| 1262 // A(:,:) -- We are deleting columns and rows, so the result | |
| 1263 // is [](0x0). | |
| 1264 | |
| 1265 resize_no_fill (0, 0); | |
| 1266 return; | |
| 1267 } | |
| 1268 | |
| 1269 if (idx_j.is_colon_equiv (nc, 1)) | |
| 1270 { | |
| 1271 // A(:,j) -- We are deleting columns by enumerating them, | |
| 1272 // If we enumerate all of them, we should have zero columns | |
| 1273 // with the same number of rows that we started with. | |
| 1274 | |
| 1275 resize_no_fill (nr, 0); | |
| 1276 return; | |
| 1277 } | |
| 1278 } | |
| 1279 | |
| 1280 if (idx_j.is_colon () && idx_i.is_colon_equiv (nr, 1)) | |
| 1281 { | |
| 1282 // A(i,:) -- We are deleting rows by enumerating them. If we | |
| 1283 // enumerate all of them, we should have zero rows with the | |
| 1284 // same number of columns that we started with. | |
| 1285 | |
| 1286 resize_no_fill (0, nc); | |
| 1287 return; | |
| 1288 } | |
| 1289 | |
| 1290 if (idx_i.is_colon_equiv (nr, 1)) | |
| 1291 { | |
| 1292 if (idx_j.is_colon_equiv (nc, 1)) | |
| 1293 resize_no_fill (0, 0); | |
| 1294 else | |
| 1295 { | |
| 1296 idx_j.sort (true); | |
| 1297 | |
| 5275 | 1298 octave_idx_type num_to_delete = idx_j.length (nc); |
| 5164 | 1299 |
| 1300 if (num_to_delete != 0) | |
| 1301 { | |
| 1302 if (nr == 1 && num_to_delete == nc) | |
| 1303 resize_no_fill (0, 0); | |
| 1304 else | |
| 1305 { | |
| 5275 | 1306 octave_idx_type new_nc = nc; |
| 5681 | 1307 octave_idx_type new_nnz = nnz (); |
| 5275 | 1308 |
| 1309 octave_idx_type iidx = 0; | |
| 1310 | |
| 1311 for (octave_idx_type j = 0; j < nc; j++) | |
| 5164 | 1312 { |
| 1313 OCTAVE_QUIT; | |
| 1314 | |
| 1315 if (j == idx_j.elem (iidx)) | |
| 1316 { | |
| 1317 iidx++; | |
| 1318 new_nc--; | |
| 1319 | |
| 5681 | 1320 new_nnz -= cidx(j+1) - cidx(j); |
| 5164 | 1321 |
| 1322 if (iidx == num_to_delete) | |
| 1323 break; | |
| 1324 } | |
| 1325 } | |
| 1326 | |
| 1327 if (new_nc > 0) | |
| 1328 { | |
| 1329 const Sparse<T> tmp (*this); | |
| 1330 --rep->count; | |
| 1331 rep = new typename Sparse<T>::SparseRep (nr, new_nc, | |
| 5681 | 1332 new_nnz); |
| 5275 | 1333 octave_idx_type ii = 0; |
| 1334 octave_idx_type jj = 0; | |
| 5164 | 1335 iidx = 0; |
| 1336 cidx(0) = 0; | |
| 5275 | 1337 for (octave_idx_type j = 0; j < nc; j++) |
| 5164 | 1338 { |
| 1339 OCTAVE_QUIT; | |
| 1340 | |
| 1341 if (iidx < num_to_delete && j == idx_j.elem (iidx)) | |
| 1342 iidx++; | |
| 1343 else | |
| 1344 { | |
| 5275 | 1345 for (octave_idx_type i = tmp.cidx(j); |
| 5164 | 1346 i < tmp.cidx(j+1); i++) |
| 1347 { | |
| 1348 data(jj) = tmp.data(i); | |
| 1349 ridx(jj++) = tmp.ridx(i); | |
| 1350 } | |
| 1351 cidx(++ii) = jj; | |
| 1352 } | |
| 1353 } | |
| 1354 | |
| 1355 dimensions.resize (2); | |
| 1356 dimensions(1) = new_nc; | |
| 1357 } | |
| 1358 else | |
| 1359 (*current_liboctave_error_handler) | |
| 1360 ("A(idx) = []: index out of range"); | |
| 1361 } | |
| 1362 } | |
| 1363 } | |
| 1364 } | |
| 1365 else if (idx_j.is_colon_equiv (nc, 1)) | |
| 1366 { | |
| 1367 if (idx_i.is_colon_equiv (nr, 1)) | |
| 1368 resize_no_fill (0, 0); | |
| 1369 else | |
| 1370 { | |
| 1371 idx_i.sort (true); | |
| 1372 | |
| 5275 | 1373 octave_idx_type num_to_delete = idx_i.length (nr); |
| 5164 | 1374 |
| 1375 if (num_to_delete != 0) | |
| 1376 { | |
| 1377 if (nc == 1 && num_to_delete == nr) | |
| 1378 resize_no_fill (0, 0); | |
| 1379 else | |
| 1380 { | |
| 5275 | 1381 octave_idx_type new_nr = nr; |
| 5681 | 1382 octave_idx_type new_nnz = nnz (); |
| 5275 | 1383 |
| 1384 octave_idx_type iidx = 0; | |
| 1385 | |
| 1386 for (octave_idx_type i = 0; i < nr; i++) | |
| 5164 | 1387 { |
| 1388 OCTAVE_QUIT; | |
| 1389 | |
| 1390 if (i == idx_i.elem (iidx)) | |
| 1391 { | |
| 1392 iidx++; | |
| 1393 new_nr--; | |
| 1394 | |
| 5681 | 1395 for (octave_idx_type j = 0; j < nnz (); j++) |
| 5164 | 1396 if (ridx(j) == i) |
| 5681 | 1397 new_nnz--; |
| 5164 | 1398 |
| 1399 if (iidx == num_to_delete) | |
| 1400 break; | |
| 1401 } | |
| 1402 } | |
| 1403 | |
| 1404 if (new_nr > 0) | |
| 1405 { | |
| 1406 const Sparse<T> tmp (*this); | |
| 1407 --rep->count; | |
| 1408 rep = new typename Sparse<T>::SparseRep (new_nr, nc, | |
| 5681 | 1409 new_nnz); |
| 5164 | 1410 |
| 5275 | 1411 octave_idx_type jj = 0; |
| 5164 | 1412 cidx(0) = 0; |
| 5275 | 1413 for (octave_idx_type i = 0; i < nc; i++) |
| 5164 | 1414 { |
| 1415 iidx = 0; | |
| 5275 | 1416 for (octave_idx_type j = tmp.cidx(i); j < tmp.cidx(i+1); j++) |
| 5164 | 1417 { |
| 1418 OCTAVE_QUIT; | |
| 1419 | |
| 5275 | 1420 octave_idx_type ri = tmp.ridx(j); |
| 5164 | 1421 |
| 1422 while (iidx < num_to_delete && | |
| 1423 ri > idx_i.elem (iidx)) | |
| 1424 { | |
| 1425 iidx++; | |
| 1426 } | |
| 1427 | |
| 1428 if (iidx == num_to_delete || | |
| 1429 ri != idx_i.elem(iidx)) | |
| 1430 { | |
| 1431 data(jj) = tmp.data(j); | |
| 1432 ridx(jj++) = ri - iidx; | |
| 1433 } | |
| 1434 } | |
| 1435 cidx(i+1) = jj; | |
| 1436 } | |
| 1437 | |
| 1438 dimensions.resize (2); | |
| 1439 dimensions(0) = new_nr; | |
| 1440 } | |
| 1441 else | |
| 1442 (*current_liboctave_error_handler) | |
| 1443 ("A(idx) = []: index out of range"); | |
| 1444 } | |
| 1445 } | |
| 1446 } | |
| 1447 } | |
| 1448 } | |
| 1449 | |
| 1450 template <class T> | |
| 1451 void | |
| 1452 Sparse<T>::maybe_delete_elements (Array<idx_vector>& ra_idx) | |
| 1453 { | |
| 1454 if (ra_idx.length () == 1) | |
| 1455 maybe_delete_elements (ra_idx(0)); | |
| 1456 else if (ra_idx.length () == 2) | |
| 1457 maybe_delete_elements (ra_idx(0), ra_idx(1)); | |
| 1458 else | |
| 1459 (*current_liboctave_error_handler) | |
| 1460 ("range error for maybe_delete_elements"); | |
| 1461 } | |
| 1462 | |
| 1463 template <class T> | |
| 1464 Sparse<T> | |
| 1465 Sparse<T>::value (void) | |
| 1466 { | |
| 1467 Sparse<T> retval; | |
| 1468 | |
| 1469 int n_idx = index_count (); | |
| 1470 | |
| 1471 if (n_idx == 2) | |
| 1472 { | |
| 1473 idx_vector *tmp = get_idx (); | |
| 1474 | |
| 1475 idx_vector idx_i = tmp[0]; | |
| 1476 idx_vector idx_j = tmp[1]; | |
| 1477 | |
| 1478 retval = index (idx_i, idx_j); | |
| 1479 } | |
| 1480 else if (n_idx == 1) | |
| 1481 { | |
| 1482 retval = index (idx[0]); | |
| 1483 } | |
| 1484 else | |
| 1485 (*current_liboctave_error_handler) | |
| 1486 ("Sparse<T>::value: invalid number of indices specified"); | |
| 1487 | |
| 1488 clear_index (); | |
| 1489 | |
| 1490 return retval; | |
| 1491 } | |
| 1492 | |
| 1493 template <class T> | |
| 1494 Sparse<T> | |
| 1495 Sparse<T>::index (idx_vector& idx_arg, int resize_ok) const | |
| 1496 { | |
| 1497 Sparse<T> retval; | |
| 1498 | |
| 1499 assert (ndims () == 2); | |
| 1500 | |
| 5275 | 1501 octave_idx_type nr = dim1 (); |
| 1502 octave_idx_type nc = dim2 (); | |
| 5681 | 1503 octave_idx_type nz = nnz (); |
| 5275 | 1504 |
| 1505 octave_idx_type orig_len = nr * nc; | |
| 5164 | 1506 |
| 1507 dim_vector idx_orig_dims = idx_arg.orig_dimensions (); | |
| 1508 | |
| 5275 | 1509 octave_idx_type idx_orig_rows = idx_arg.orig_rows (); |
| 1510 octave_idx_type idx_orig_columns = idx_arg.orig_columns (); | |
| 5164 | 1511 |
| 1512 if (idx_orig_dims.length () > 2) | |
| 1513 (*current_liboctave_error_handler) | |
| 1514 ("Sparse<T>::index: Can not index Sparse<T> with an N-D Array"); | |
| 1515 else if (idx_arg.is_colon ()) | |
| 1516 { | |
| 1517 // Fast magic colon processing. | |
| 1518 retval = Sparse<T> (nr * nc, 1, nz); | |
| 1519 | |
| 5275 | 1520 for (octave_idx_type i = 0; i < nc; i++) |
| 1521 for (octave_idx_type j = cidx(i); j < cidx(i+1); j++) | |
| 5164 | 1522 { |
| 1523 OCTAVE_QUIT; | |
| 1524 retval.xdata(j) = data(j); | |
| 1525 retval.xridx(j) = ridx(j) + i * nr; | |
| 1526 } | |
| 1527 retval.xcidx(0) = 0; | |
| 1528 retval.xcidx(1) = nz; | |
| 1529 } | |
| 1530 else if (nr == 1 && nc == 1) | |
| 1531 { | |
| 1532 // You have to be pretty sick to get to this bit of code, | |
| 1533 // since you have a scalar stored as a sparse matrix, and | |
| 1534 // then want to make a dense matrix with sparse | |
| 1535 // representation. Ok, we'll do it, but you deserve what | |
| 1536 // you get!! | |
| 5275 | 1537 octave_idx_type n = idx_arg.freeze (length (), "sparse vector", resize_ok); |
| 5164 | 1538 if (n == 0) |
| 1539 if (idx_arg.one_zero_only ()) | |
| 1540 retval = Sparse<T> (dim_vector (0, 0)); | |
| 1541 else | |
| 7299 | 1542 retval = Sparse<T> (idx_orig_dims); |
| 5164 | 1543 else if (nz < 1) |
| 1544 if (n >= idx_orig_dims.numel ()) | |
| 1545 retval = Sparse<T> (idx_orig_dims); | |
| 1546 else | |
| 1547 retval = Sparse<T> (dim_vector (n, 1)); | |
| 1548 else if (n >= idx_orig_dims.numel ()) | |
| 1549 { | |
| 1550 T el = elem (0); | |
| 5275 | 1551 octave_idx_type new_nr = idx_orig_rows; |
| 1552 octave_idx_type new_nc = idx_orig_columns; | |
| 1553 for (octave_idx_type i = 2; i < idx_orig_dims.length (); i++) | |
| 5164 | 1554 new_nc *= idx_orig_dims (i); |
| 1555 | |
| 1556 retval = Sparse<T> (new_nr, new_nc, idx_arg.ones_count ()); | |
| 1557 | |
| 5275 | 1558 octave_idx_type ic = 0; |
| 1559 for (octave_idx_type i = 0; i < n; i++) | |
| 5164 | 1560 { |
| 1561 if (i % new_nr == 0) | |
| 7322 | 1562 retval.xcidx(i / new_nr) = ic; |
| 5164 | 1563 |
| 5275 | 1564 octave_idx_type ii = idx_arg.elem (i); |
| 5164 | 1565 if (ii == 0) |
| 1566 { | |
| 1567 OCTAVE_QUIT; | |
| 1568 retval.xdata(ic) = el; | |
| 1569 retval.xridx(ic++) = i % new_nr; | |
| 1570 } | |
| 1571 } | |
| 1572 retval.xcidx (new_nc) = ic; | |
| 1573 } | |
| 1574 else | |
| 1575 { | |
| 1576 T el = elem (0); | |
| 1577 retval = Sparse<T> (n, 1, nz); | |
| 1578 | |
| 5275 | 1579 for (octave_idx_type i = 0; i < nz; i++) |
| 5164 | 1580 { |
| 1581 OCTAVE_QUIT; | |
| 1582 retval.xdata(i) = el; | |
| 1583 retval.xridx(i) = i; | |
| 1584 } | |
| 1585 retval.xcidx(0) = 0; | |
| 1586 retval.xcidx(1) = n; | |
| 1587 } | |
| 1588 } | |
| 1589 else if (nr == 1 || nc == 1) | |
| 1590 { | |
| 1591 // If indexing a vector with a matrix, return value has same | |
| 1592 // shape as the index. Otherwise, it has same orientation as | |
| 1593 // indexed object. | |
| 5275 | 1594 octave_idx_type len = length (); |
| 1595 octave_idx_type n = idx_arg.freeze (len, "sparse vector", resize_ok); | |
| 5164 | 1596 |
| 1597 if (n == 0) | |
| 1598 if (nr == 1) | |
| 1599 retval = Sparse<T> (dim_vector (1, 0)); | |
| 1600 else | |
| 1601 retval = Sparse<T> (dim_vector (0, 1)); | |
| 1602 else if (nz < 1) | |
| 1603 if ((n != 0 && idx_arg.one_zero_only ()) | |
| 1604 || idx_orig_rows == 1 || idx_orig_columns == 1) | |
| 1605 retval = Sparse<T> ((nr == 1 ? 1 : n), (nr == 1 ? n : 1)); | |
| 1606 else | |
| 1607 retval = Sparse<T> (idx_orig_dims); | |
| 1608 else | |
| 1609 { | |
| 1610 | |
| 5604 | 1611 octave_idx_type new_nzmx = 0; |
| 5164 | 1612 if (nr == 1) |
| 5275 | 1613 for (octave_idx_type i = 0; i < n; i++) |
| 5164 | 1614 { |
| 1615 OCTAVE_QUIT; | |
| 1616 | |
| 5275 | 1617 octave_idx_type ii = idx_arg.elem (i); |
| 5164 | 1618 if (ii < len) |
| 1619 if (cidx(ii) != cidx(ii+1)) | |
| 5604 | 1620 new_nzmx++; |
| 5164 | 1621 } |
| 1622 else | |
| 5275 | 1623 for (octave_idx_type i = 0; i < n; i++) |
| 5164 | 1624 { |
| 5275 | 1625 octave_idx_type ii = idx_arg.elem (i); |
| 5164 | 1626 if (ii < len) |
| 5275 | 1627 for (octave_idx_type j = 0; j < nz; j++) |
| 5164 | 1628 { |
| 1629 OCTAVE_QUIT; | |
| 1630 | |
| 1631 if (ridx(j) == ii) | |
| 5604 | 1632 new_nzmx++; |
| 5164 | 1633 if (ridx(j) >= ii) |
| 1634 break; | |
| 1635 } | |
| 1636 } | |
| 1637 | |
| 1638 if (idx_arg.one_zero_only () || idx_orig_rows == 1 || | |
| 1639 idx_orig_columns == 1) | |
| 1640 { | |
| 1641 if (nr == 1) | |
| 1642 { | |
| 5604 | 1643 retval = Sparse<T> (1, n, new_nzmx); |
| 5275 | 1644 octave_idx_type jj = 0; |
| 5164 | 1645 retval.xcidx(0) = 0; |
| 5275 | 1646 for (octave_idx_type i = 0; i < n; i++) |
| 5164 | 1647 { |
| 1648 OCTAVE_QUIT; | |
| 1649 | |
| 5275 | 1650 octave_idx_type ii = idx_arg.elem (i); |
| 5164 | 1651 if (ii < len) |
| 1652 if (cidx(ii) != cidx(ii+1)) | |
| 1653 { | |
| 1654 retval.xdata(jj) = data(cidx(ii)); | |
| 1655 retval.xridx(jj++) = 0; | |
| 1656 } | |
| 1657 retval.xcidx(i+1) = jj; | |
| 1658 } | |
| 1659 } | |
| 1660 else | |
| 1661 { | |
| 5604 | 1662 retval = Sparse<T> (n, 1, new_nzmx); |
| 5164 | 1663 retval.xcidx(0) = 0; |
| 5604 | 1664 retval.xcidx(1) = new_nzmx; |
| 5275 | 1665 octave_idx_type jj = 0; |
| 1666 for (octave_idx_type i = 0; i < n; i++) | |
| 5164 | 1667 { |
| 5275 | 1668 octave_idx_type ii = idx_arg.elem (i); |
| 5164 | 1669 if (ii < len) |
| 5275 | 1670 for (octave_idx_type j = 0; j < nz; j++) |
| 5164 | 1671 { |
| 1672 OCTAVE_QUIT; | |
| 1673 | |
| 1674 if (ridx(j) == ii) | |
| 1675 { | |
| 1676 retval.xdata(jj) = data(j); | |
| 1677 retval.xridx(jj++) = i; | |
| 1678 } | |
| 1679 if (ridx(j) >= ii) | |
| 1680 break; | |
| 1681 } | |
| 1682 } | |
| 1683 } | |
| 1684 } | |
| 1685 else | |
| 1686 { | |
| 5275 | 1687 octave_idx_type new_nr; |
| 1688 octave_idx_type new_nc; | |
| 5164 | 1689 if (n >= idx_orig_dims.numel ()) |
| 1690 { | |
| 1691 new_nr = idx_orig_rows; | |
| 1692 new_nc = idx_orig_columns; | |
| 1693 } | |
| 1694 else | |
| 1695 { | |
| 1696 new_nr = n; | |
| 1697 new_nc = 1; | |
| 1698 } | |
| 1699 | |
| 5604 | 1700 retval = Sparse<T> (new_nr, new_nc, new_nzmx); |
| 5164 | 1701 |
| 1702 if (nr == 1) | |
| 1703 { | |
| 5275 | 1704 octave_idx_type jj = 0; |
| 5164 | 1705 retval.xcidx(0) = 0; |
| 5275 | 1706 for (octave_idx_type i = 0; i < n; i++) |
| 5164 | 1707 { |
| 1708 OCTAVE_QUIT; | |
| 1709 | |
| 5275 | 1710 octave_idx_type ii = idx_arg.elem (i); |
| 5164 | 1711 if (ii < len) |
| 1712 if (cidx(ii) != cidx(ii+1)) | |
| 1713 { | |
| 1714 retval.xdata(jj) = data(cidx(ii)); | |
| 1715 retval.xridx(jj++) = 0; | |
| 1716 } | |
| 1717 retval.xcidx(i/new_nr+1) = jj; | |
| 1718 } | |
| 1719 } | |
| 1720 else | |
| 1721 { | |
| 5275 | 1722 octave_idx_type jj = 0; |
| 5164 | 1723 retval.xcidx(0) = 0; |
| 5275 | 1724 for (octave_idx_type i = 0; i < n; i++) |
| 5164 | 1725 { |
| 5275 | 1726 octave_idx_type ii = idx_arg.elem (i); |
| 5164 | 1727 if (ii < len) |
| 5275 | 1728 for (octave_idx_type j = 0; j < nz; j++) |
| 5164 | 1729 { |
| 1730 OCTAVE_QUIT; | |
| 1731 | |
| 1732 if (ridx(j) == ii) | |
| 1733 { | |
| 1734 retval.xdata(jj) = data(j); | |
| 1735 retval.xridx(jj++) = i; | |
| 1736 } | |
| 1737 if (ridx(j) >= ii) | |
| 1738 break; | |
| 1739 } | |
| 1740 retval.xcidx(i/new_nr+1) = jj; | |
| 1741 } | |
| 1742 } | |
| 1743 } | |
| 1744 } | |
| 1745 } | |
| 1746 else | |
| 1747 { | |
| 5781 | 1748 if (! (idx_arg.one_zero_only () |
| 1749 && idx_orig_rows == nr | |
| 1750 && idx_orig_columns == nc)) | |
| 1751 (*current_liboctave_warning_with_id_handler) | |
| 1752 ("Octave:fortran-indexing", "single index used for sparse matrix"); | |
| 5164 | 1753 |
| 1754 // This code is only for indexing matrices. The vector | |
| 1755 // cases are handled above. | |
| 1756 | |
| 1757 idx_arg.freeze (nr * nc, "matrix", resize_ok); | |
| 1758 | |
| 1759 if (idx_arg) | |
| 1760 { | |
| 5275 | 1761 octave_idx_type result_nr = idx_orig_rows; |
| 1762 octave_idx_type result_nc = idx_orig_columns; | |
| 5164 | 1763 |
| 1764 if (idx_arg.one_zero_only ()) | |
| 1765 { | |
| 1766 result_nr = idx_arg.ones_count (); | |
| 1767 result_nc = (result_nr > 0 ? 1 : 0); | |
| 1768 } | |
| 1769 | |
| 1770 if (nz < 1) | |
| 1771 retval = Sparse<T> (result_nr, result_nc); | |
| 1772 else | |
| 1773 { | |
| 1774 // Count number of non-zero elements | |
| 5604 | 1775 octave_idx_type new_nzmx = 0; |
| 5275 | 1776 octave_idx_type kk = 0; |
| 1777 for (octave_idx_type j = 0; j < result_nc; j++) | |
| 5164 | 1778 { |
| 5275 | 1779 for (octave_idx_type i = 0; i < result_nr; i++) |
| 5164 | 1780 { |
| 1781 OCTAVE_QUIT; | |
| 1782 | |
| 5275 | 1783 octave_idx_type ii = idx_arg.elem (kk++); |
| 5164 | 1784 if (ii < orig_len) |
| 1785 { | |
| 5275 | 1786 octave_idx_type fr = ii % nr; |
| 1787 octave_idx_type fc = (ii - fr) / nr; | |
| 1788 for (octave_idx_type k = cidx(fc); k < cidx(fc+1); k++) | |
| 5164 | 1789 { |
| 1790 if (ridx(k) == fr) | |
| 5604 | 1791 new_nzmx++; |
| 5164 | 1792 if (ridx(k) >= fr) |
| 1793 break; | |
| 1794 } | |
| 1795 } | |
| 1796 } | |
| 1797 } | |
| 1798 | |
| 5604 | 1799 retval = Sparse<T> (result_nr, result_nc, new_nzmx); |
| 5164 | 1800 |
| 1801 kk = 0; | |
| 5275 | 1802 octave_idx_type jj = 0; |
| 5164 | 1803 retval.xcidx(0) = 0; |
| 5275 | 1804 for (octave_idx_type j = 0; j < result_nc; j++) |
| 5164 | 1805 { |
| 5275 | 1806 for (octave_idx_type i = 0; i < result_nr; i++) |
| 5164 | 1807 { |
| 1808 OCTAVE_QUIT; | |
| 1809 | |
| 5275 | 1810 octave_idx_type ii = idx_arg.elem (kk++); |
| 5164 | 1811 if (ii < orig_len) |
| 1812 { | |
| 5275 | 1813 octave_idx_type fr = ii % nr; |
| 1814 octave_idx_type fc = (ii - fr) / nr; | |
| 1815 for (octave_idx_type k = cidx(fc); k < cidx(fc+1); k++) | |
| 5164 | 1816 { |
| 1817 if (ridx(k) == fr) | |
| 1818 { | |
| 1819 retval.xdata(jj) = data(k); | |
| 1820 retval.xridx(jj++) = i; | |
| 1821 } | |
| 1822 if (ridx(k) >= fr) | |
| 1823 break; | |
| 1824 } | |
| 1825 } | |
| 1826 } | |
| 1827 retval.xcidx(j+1) = jj; | |
| 1828 } | |
| 1829 } | |
| 1830 // idx_vector::freeze() printed an error message for us. | |
| 1831 } | |
| 1832 } | |
| 1833 | |
| 1834 return retval; | |
| 1835 } | |
| 1836 | |
| 6988 | 1837 struct |
| 1838 idx_node | |
| 1839 { | |
| 1840 octave_idx_type i; | |
| 1841 struct idx_node *next; | |
| 1842 }; | |
| 1843 | |
| 5164 | 1844 template <class T> |
| 1845 Sparse<T> | |
| 1846 Sparse<T>::index (idx_vector& idx_i, idx_vector& idx_j, int resize_ok) const | |
| 1847 { | |
| 1848 Sparse<T> retval; | |
| 1849 | |
| 1850 assert (ndims () == 2); | |
| 1851 | |
| 5275 | 1852 octave_idx_type nr = dim1 (); |
| 1853 octave_idx_type nc = dim2 (); | |
| 1854 | |
| 1855 octave_idx_type n = idx_i.freeze (nr, "row", resize_ok); | |
| 1856 octave_idx_type m = idx_j.freeze (nc, "column", resize_ok); | |
| 5164 | 1857 |
| 1858 if (idx_i && idx_j) | |
| 1859 { | |
| 1860 if (idx_i.orig_empty () || idx_j.orig_empty () || n == 0 || m == 0) | |
| 1861 { | |
| 1862 retval.resize_no_fill (n, m); | |
| 1863 } | |
| 5681 | 1864 else |
| 5164 | 1865 { |
| 5681 | 1866 int idx_i_colon = idx_i.is_colon_equiv (nr); |
| 1867 int idx_j_colon = idx_j.is_colon_equiv (nc); | |
| 1868 | |
| 1869 if (idx_i_colon && idx_j_colon) | |
| 1870 { | |
| 1871 retval = *this; | |
| 1872 } | |
| 1873 else | |
| 5164 | 1874 { |
| 5681 | 1875 // Identify if the indices have any repeated values |
| 1876 bool permutation = true; | |
| 1877 | |
| 1878 OCTAVE_LOCAL_BUFFER (octave_idx_type, itmp, | |
| 1879 (nr > nc ? nr : nc)); | |
| 1880 octave_sort<octave_idx_type> sort; | |
| 1881 | |
| 1882 if (n > nr || m > nc) | |
| 1883 permutation = false; | |
| 1884 | |
| 1885 if (permutation && ! idx_i_colon) | |
| 1886 { | |
| 1887 // Can't use something like | |
| 1888 // idx_vector tmp_idx = idx_i; | |
| 1889 // tmp_idx.sort (true); | |
| 1890 // if (tmp_idx.length(nr) != n) | |
| 1891 // permutation = false; | |
| 1892 // here as there is no make_unique function | |
| 1893 // for idx_vector type. | |
| 1894 for (octave_idx_type i = 0; i < n; i++) | |
| 1895 itmp [i] = idx_i.elem (i); | |
| 1896 sort.sort (itmp, n); | |
| 1897 for (octave_idx_type i = 1; i < n; i++) | |
| 1898 if (itmp[i-1] == itmp[i]) | |
| 1899 { | |
| 1900 permutation = false; | |
| 1901 break; | |
| 1902 } | |
| 1903 } | |
| 1904 if (permutation && ! idx_j_colon) | |
| 1905 { | |
| 1906 for (octave_idx_type i = 0; i < m; i++) | |
| 1907 itmp [i] = idx_j.elem (i); | |
| 1908 sort.sort (itmp, m); | |
| 1909 for (octave_idx_type i = 1; i < m; i++) | |
| 1910 if (itmp[i-1] == itmp[i]) | |
| 1911 { | |
| 1912 permutation = false; | |
| 1913 break; | |
| 1914 } | |
| 1915 } | |
| 1916 | |
| 1917 if (permutation) | |
| 5164 | 1918 { |
| 5681 | 1919 // Special case permutation like indexing for speed |
| 1920 retval = Sparse<T> (n, m, nnz ()); | |
| 1921 octave_idx_type *ri = retval.xridx (); | |
| 1922 | |
| 5766 | 1923 std::vector<T> X (n); |
| 5681 | 1924 for (octave_idx_type i = 0; i < nr; i++) |
| 1925 itmp [i] = -1; | |
| 1926 for (octave_idx_type i = 0; i < n; i++) | |
| 1927 itmp[idx_i.elem(i)] = i; | |
| 1928 | |
| 1929 octave_idx_type kk = 0; | |
| 1930 retval.xcidx(0) = 0; | |
| 1931 for (octave_idx_type j = 0; j < m; j++) | |
| 1932 { | |
| 1933 octave_idx_type jj = idx_j.elem (j); | |
| 1934 for (octave_idx_type i = cidx(jj); i < cidx(jj+1); i++) | |
| 1935 { | |
| 6988 | 1936 OCTAVE_QUIT; |
| 1937 | |
| 5681 | 1938 octave_idx_type ii = itmp [ridx(i)]; |
| 1939 if (ii >= 0) | |
| 1940 { | |
| 1941 X [ii] = data (i); | |
| 1942 retval.xridx (kk++) = ii; | |
| 1943 } | |
| 1944 } | |
| 1945 sort.sort (ri + retval.xcidx (j), kk - retval.xcidx (j)); | |
| 1946 for (octave_idx_type p = retval.xcidx (j); p < kk; p++) | |
| 1947 retval.xdata (p) = X [retval.xridx (p)]; | |
| 1948 retval.xcidx(j+1) = kk; | |
| 1949 } | |
| 1950 retval.maybe_compress (); | |
| 1951 } | |
| 1952 else | |
| 1953 { | |
| 6988 | 1954 OCTAVE_LOCAL_BUFFER (struct idx_node, nodes, n); |
| 1955 OCTAVE_LOCAL_BUFFER (octave_idx_type, start_nodes, nr); | |
| 1956 | |
| 1957 for (octave_idx_type i = 0; i < nr; i++) | |
| 1958 start_nodes[i] = -1; | |
| 1959 | |
| 1960 for (octave_idx_type i = 0; i < n; i++) | |
| 1961 { | |
| 1962 octave_idx_type ii = idx_i.elem (i); | |
| 1963 nodes[i].i = i; | |
| 1964 nodes[i].next = 0; | |
| 1965 | |
| 1966 octave_idx_type node = start_nodes[ii]; | |
| 1967 if (node == -1) | |
| 1968 start_nodes[ii] = i; | |
| 1969 else | |
| 1970 { | |
| 7322 | 1971 while (nodes[node].next) |
| 1972 node = nodes[node].next->i; | |
| 1973 nodes[node].next = nodes + i; | |
| 6988 | 1974 } |
| 1975 } | |
| 1976 | |
| 5681 | 1977 // First count the number of non-zero elements |
| 1978 octave_idx_type new_nzmx = 0; | |
| 1979 for (octave_idx_type j = 0; j < m; j++) | |
| 5164 | 1980 { |
| 5681 | 1981 octave_idx_type jj = idx_j.elem (j); |
| 6988 | 1982 |
| 1983 if (jj < nc) | |
| 5164 | 1984 { |
| 6988 | 1985 for (octave_idx_type i = cidx(jj); |
| 1986 i < cidx(jj+1); i++) | |
| 5681 | 1987 { |
| 6988 | 1988 OCTAVE_QUIT; |
| 1989 | |
| 1990 octave_idx_type ii = start_nodes [ridx(i)]; | |
| 1991 | |
| 1992 if (ii >= 0) | |
| 5681 | 1993 { |
| 6988 | 1994 struct idx_node inode = nodes[ii]; |
| 1995 | |
| 1996 while (true) | |
| 1997 { | |
| 7326 | 1998 if (idx_i.elem (inode.i) < nr) |
| 6988 | 1999 new_nzmx ++; |
| 2000 if (inode.next == 0) | |
| 2001 break; | |
| 2002 else | |
| 2003 inode = *inode.next; | |
| 2004 } | |
| 5681 | 2005 } |
| 2006 } | |
| 5164 | 2007 } |
| 2008 } | |
| 5681 | 2009 |
| 6988 | 2010 std::vector<T> X (n); |
| 5681 | 2011 retval = Sparse<T> (n, m, new_nzmx); |
| 6988 | 2012 octave_idx_type *ri = retval.xridx (); |
| 5681 | 2013 |
| 2014 octave_idx_type kk = 0; | |
| 2015 retval.xcidx(0) = 0; | |
| 2016 for (octave_idx_type j = 0; j < m; j++) | |
| 2017 { | |
| 2018 octave_idx_type jj = idx_j.elem (j); | |
| 6988 | 2019 if (jj < nc) |
| 5681 | 2020 { |
| 6988 | 2021 for (octave_idx_type i = cidx(jj); |
| 2022 i < cidx(jj+1); i++) | |
| 5681 | 2023 { |
| 6988 | 2024 OCTAVE_QUIT; |
| 2025 | |
| 2026 octave_idx_type ii = start_nodes [ridx(i)]; | |
| 2027 | |
| 2028 if (ii >= 0) | |
| 5681 | 2029 { |
| 6988 | 2030 struct idx_node inode = nodes[ii]; |
| 2031 | |
| 2032 while (true) | |
| 5681 | 2033 { |
| 7326 | 2034 if (idx_i.elem (inode.i) < nr) |
| 6988 | 2035 { |
| 2036 X [inode.i] = data (i); | |
| 2037 retval.xridx (kk++) = inode.i; | |
| 2038 } | |
| 2039 | |
| 2040 if (inode.next == 0) | |
| 2041 break; | |
| 2042 else | |
| 2043 inode = *inode.next; | |
| 5681 | 2044 } |
| 2045 } | |
| 2046 } | |
| 6988 | 2047 sort.sort (ri + retval.xcidx (j), |
| 2048 kk - retval.xcidx (j)); | |
| 2049 for (octave_idx_type p = retval.xcidx (j); | |
| 2050 p < kk; p++) | |
| 2051 retval.xdata (p) = X [retval.xridx (p)]; | |
| 2052 retval.xcidx(j+1) = kk; | |
| 5681 | 2053 } |
| 2054 } | |
| 5164 | 2055 } |
| 2056 } | |
| 2057 } | |
| 2058 } | |
| 2059 // idx_vector::freeze() printed an error message for us. | |
| 2060 | |
| 2061 return retval; | |
| 2062 } | |
| 2063 | |
| 2064 template <class T> | |
| 2065 Sparse<T> | |
| 2066 Sparse<T>::index (Array<idx_vector>& ra_idx, int resize_ok) const | |
| 2067 { | |
| 2068 | |
| 2069 if (ra_idx.length () != 2) | |
| 2070 { | |
| 2071 (*current_liboctave_error_handler) ("range error for index"); | |
| 2072 return *this; | |
| 2073 } | |
| 2074 | |
| 2075 return index (ra_idx (0), ra_idx (1), resize_ok); | |
| 2076 } | |
| 2077 | |
| 5775 | 2078 // FIXME |
| 5164 | 2079 // Unfortunately numel can overflow for very large but very sparse matrices. |
| 2080 // For now just flag an error when this happens. | |
| 2081 template <class LT, class RT> | |
| 2082 int | |
| 2083 assign1 (Sparse<LT>& lhs, const Sparse<RT>& rhs) | |
| 2084 { | |
| 2085 int retval = 1; | |
| 2086 | |
| 2087 idx_vector *idx_tmp = lhs.get_idx (); | |
| 2088 | |
| 2089 idx_vector lhs_idx = idx_tmp[0]; | |
| 2090 | |
| 5275 | 2091 octave_idx_type lhs_len = lhs.numel (); |
| 2092 octave_idx_type rhs_len = rhs.numel (); | |
| 5164 | 2093 |
| 5828 | 2094 uint64_t long_lhs_len = |
| 2095 static_cast<uint64_t> (lhs.rows ()) * | |
| 2096 static_cast<uint64_t> (lhs.cols ()); | |
| 2097 | |
| 2098 uint64_t long_rhs_len = | |
| 2099 static_cast<uint64_t> (rhs.rows ()) * | |
| 2100 static_cast<uint64_t> (rhs.cols ()); | |
| 2101 | |
| 2102 if (long_rhs_len != static_cast<uint64_t>(rhs_len) || | |
| 2103 long_lhs_len != static_cast<uint64_t>(lhs_len)) | |
| 5164 | 2104 { |
| 2105 (*current_liboctave_error_handler) | |
| 2106 ("A(I) = X: Matrix dimensions too large to ensure correct\n", | |
| 2107 "operation. This is an limitation that should be removed\n", | |
| 2108 "in the future."); | |
| 2109 | |
| 2110 lhs.clear_index (); | |
| 2111 return 0; | |
| 2112 } | |
| 2113 | |
| 5275 | 2114 octave_idx_type nr = lhs.rows (); |
| 2115 octave_idx_type nc = lhs.cols (); | |
| 5681 | 2116 octave_idx_type nz = lhs.nnz (); |
| 5275 | 2117 |
| 5781 | 2118 octave_idx_type n = lhs_idx.freeze (lhs_len, "vector", true); |
| 5164 | 2119 |
| 2120 if (n != 0) | |
| 2121 { | |
| 5275 | 2122 octave_idx_type max_idx = lhs_idx.max () + 1; |
| 5164 | 2123 max_idx = max_idx < lhs_len ? lhs_len : max_idx; |
| 2124 | |
| 2125 // Take a constant copy of lhs. This means that elem won't | |
| 2126 // create missing elements. | |
| 2127 const Sparse<LT> c_lhs (lhs); | |
| 2128 | |
| 2129 if (rhs_len == n) | |
| 2130 { | |
| 5681 | 2131 octave_idx_type new_nzmx = lhs.nnz (); |
| 5164 | 2132 |
| 5603 | 2133 OCTAVE_LOCAL_BUFFER (octave_idx_type, rhs_idx, n); |
| 2134 if (! lhs_idx.is_colon ()) | |
| 2135 { | |
| 2136 // Ok here we have to be careful with the indexing, | |
| 2137 // to treat cases like "a([3,2,1]) = b", and still | |
| 2138 // handle the need for strict sorting of the sparse | |
| 2139 // elements. | |
| 2140 OCTAVE_LOCAL_BUFFER (octave_idx_vector_sort *, sidx, n); | |
| 2141 OCTAVE_LOCAL_BUFFER (octave_idx_vector_sort, sidxX, n); | |
| 2142 | |
| 2143 for (octave_idx_type i = 0; i < n; i++) | |
| 2144 { | |
| 2145 sidx[i] = &sidxX[i]; | |
| 2146 sidx[i]->i = lhs_idx.elem(i); | |
| 2147 sidx[i]->idx = i; | |
| 2148 } | |
| 2149 | |
| 2150 OCTAVE_QUIT; | |
| 2151 octave_sort<octave_idx_vector_sort *> | |
| 2152 sort (octave_idx_vector_comp); | |
| 2153 | |
| 2154 sort.sort (sidx, n); | |
| 2155 | |
| 2156 intNDArray<octave_idx_type> new_idx (dim_vector (n,1)); | |
| 2157 | |
| 2158 for (octave_idx_type i = 0; i < n; i++) | |
| 2159 { | |
| 2160 new_idx.xelem(i) = sidx[i]->i + 1; | |
| 2161 rhs_idx[i] = sidx[i]->idx; | |
| 2162 } | |
| 2163 | |
| 2164 lhs_idx = idx_vector (new_idx); | |
| 2165 } | |
| 2166 else | |
| 2167 for (octave_idx_type i = 0; i < n; i++) | |
| 2168 rhs_idx[i] = i; | |
| 2169 | |
| 5164 | 2170 // First count the number of non-zero elements |
| 5275 | 2171 for (octave_idx_type i = 0; i < n; i++) |
| 5164 | 2172 { |
| 2173 OCTAVE_QUIT; | |
| 2174 | |
| 5275 | 2175 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
|
2176 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
|
2177 continue; |
| 5164 | 2178 if (ii < lhs_len && c_lhs.elem(ii) != LT ()) |
| 5604 | 2179 new_nzmx--; |
| 5603 | 2180 if (rhs.elem(rhs_idx[i]) != RT ()) |
| 5604 | 2181 new_nzmx++; |
| 5164 | 2182 } |
| 2183 | |
| 2184 if (nr > 1) | |
| 2185 { | |
| 7246 | 2186 Sparse<LT> tmp ((max_idx > nr ? max_idx : nr), 1, new_nzmx); |
| 5164 | 2187 tmp.cidx(0) = 0; |
| 5681 | 2188 tmp.cidx(1) = new_nzmx; |
| 5164 | 2189 |
| 5275 | 2190 octave_idx_type i = 0; |
| 2191 octave_idx_type ii = 0; | |
| 5164 | 2192 if (i < nz) |
| 2193 ii = c_lhs.ridx(i); | |
| 2194 | |
| 5275 | 2195 octave_idx_type j = 0; |
| 2196 octave_idx_type jj = lhs_idx.elem(j); | |
| 2197 | |
| 2198 octave_idx_type kk = 0; | |
| 5164 | 2199 |
| 2200 while (j < n || i < nz) | |
| 2201 { | |
| 2202 if (j == n || (i < nz && ii < jj)) | |
| 2203 { | |
| 2204 tmp.xdata (kk) = c_lhs.data (i); | |
| 2205 tmp.xridx (kk++) = ii; | |
| 2206 if (++i < nz) | |
| 2207 ii = c_lhs.ridx(i); | |
| 2208 } | |
| 2209 else | |
| 2210 { | |
| 5603 | 2211 RT rtmp = rhs.elem (rhs_idx[j]); |
| 5164 | 2212 if (rtmp != RT ()) |
| 2213 { | |
| 2214 tmp.xdata (kk) = rtmp; | |
| 2215 tmp.xridx (kk++) = jj; | |
| 2216 } | |
| 2217 | |
| 2218 if (ii == jj && i < nz) | |
| 2219 if (++i < nz) | |
| 2220 ii = c_lhs.ridx(i); | |
| 2221 if (++j < n) | |
| 2222 jj = lhs_idx.elem(j); | |
| 2223 } | |
| 2224 } | |
| 2225 | |
| 2226 lhs = tmp; | |
| 2227 } | |
| 2228 else | |
| 2229 { | |
| 7246 | 2230 Sparse<LT> tmp (1, (max_idx > nc ? max_idx : nc), new_nzmx); |
| 5164 | 2231 |
| 5275 | 2232 octave_idx_type i = 0; |
| 2233 octave_idx_type ii = 0; | |
| 5164 | 2234 while (ii < nc && c_lhs.cidx(ii+1) <= i) |
| 2235 ii++; | |
| 2236 | |
| 5275 | 2237 octave_idx_type j = 0; |
| 2238 octave_idx_type jj = lhs_idx.elem(j); | |
| 2239 | |
| 2240 octave_idx_type kk = 0; | |
| 2241 octave_idx_type ic = 0; | |
| 5164 | 2242 |
| 2243 while (j < n || i < nz) | |
| 2244 { | |
| 2245 if (j == n || (i < nz && ii < jj)) | |
| 2246 { | |
| 2247 while (ic <= ii) | |
| 2248 tmp.xcidx (ic++) = kk; | |
| 2249 tmp.xdata (kk) = c_lhs.data (i); | |
| 5603 | 2250 tmp.xridx (kk++) = 0; |
| 5164 | 2251 i++; |
| 2252 while (ii < nc && c_lhs.cidx(ii+1) <= i) | |
| 2253 ii++; | |
| 2254 } | |
| 2255 else | |
| 2256 { | |
| 2257 while (ic <= jj) | |
| 2258 tmp.xcidx (ic++) = kk; | |
| 2259 | |
| 5603 | 2260 RT rtmp = rhs.elem (rhs_idx[j]); |
| 5164 | 2261 if (rtmp != RT ()) |
| 5603 | 2262 { |
| 2263 tmp.xdata (kk) = rtmp; | |
| 2264 tmp.xridx (kk++) = 0; | |
| 2265 } | |
| 5164 | 2266 if (ii == jj) |
| 2267 { | |
| 2268 i++; | |
| 2269 while (ii < nc && c_lhs.cidx(ii+1) <= i) | |
| 2270 ii++; | |
| 2271 } | |
| 2272 j++; | |
| 2273 if (j < n) | |
| 2274 jj = lhs_idx.elem(j); | |
| 2275 } | |
| 2276 } | |
| 2277 | |
| 5275 | 2278 for (octave_idx_type iidx = ic; iidx < max_idx+1; iidx++) |
| 5164 | 2279 tmp.xcidx(iidx) = kk; |
| 2280 | |
| 2281 lhs = tmp; | |
| 2282 } | |
| 2283 } | |
| 2284 else if (rhs_len == 1) | |
| 2285 { | |
| 5681 | 2286 octave_idx_type new_nzmx = lhs.nnz (); |
| 5164 | 2287 RT scalar = rhs.elem (0); |
| 2288 bool scalar_non_zero = (scalar != RT ()); | |
| 5603 | 2289 lhs_idx.sort (true); |
|
11669
a6e08ecb4050
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7326
diff
changeset
|
2290 n = lhs_idx.length (n); |
| 5164 | 2291 |
| 2292 // First count the number of non-zero elements | |
| 2293 if (scalar != RT ()) | |
| 5604 | 2294 new_nzmx += n; |
| 5275 | 2295 for (octave_idx_type i = 0; i < n; i++) |
| 5164 | 2296 { |
| 2297 OCTAVE_QUIT; | |
| 2298 | |
| 5275 | 2299 octave_idx_type ii = lhs_idx.elem (i); |
| 5164 | 2300 if (ii < lhs_len && c_lhs.elem(ii) != LT ()) |
| 5604 | 2301 new_nzmx--; |
| 5164 | 2302 } |
| 2303 | |
| 2304 if (nr > 1) | |
| 2305 { | |
| 7246 | 2306 Sparse<LT> tmp ((max_idx > nr ? max_idx : nr), 1, new_nzmx); |
| 5164 | 2307 tmp.cidx(0) = 0; |
| 5681 | 2308 tmp.cidx(1) = new_nzmx; |
| 5164 | 2309 |
| 5275 | 2310 octave_idx_type i = 0; |
| 2311 octave_idx_type ii = 0; | |
| 5164 | 2312 if (i < nz) |
| 2313 ii = c_lhs.ridx(i); | |
| 2314 | |
| 5275 | 2315 octave_idx_type j = 0; |
| 2316 octave_idx_type jj = lhs_idx.elem(j); | |
| 2317 | |
| 2318 octave_idx_type kk = 0; | |
| 5164 | 2319 |
| 2320 while (j < n || i < nz) | |
| 2321 { | |
|
11669
a6e08ecb4050
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7326
diff
changeset
|
2322 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
|
2323 { |
|
a6e08ecb4050
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7326
diff
changeset
|
2324 j++; |
|
a6e08ecb4050
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7326
diff
changeset
|
2325 jj = lhs_idx.elem (j); |
|
a6e08ecb4050
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7326
diff
changeset
|
2326 continue; |
|
a6e08ecb4050
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7326
diff
changeset
|
2327 } |
| 5164 | 2328 if (j == n || (i < nz && ii < jj)) |
| 2329 { | |
| 2330 tmp.xdata (kk) = c_lhs.data (i); | |
| 2331 tmp.xridx (kk++) = ii; | |
| 2332 if (++i < nz) | |
| 2333 ii = c_lhs.ridx(i); | |
| 2334 } | |
| 2335 else | |
| 2336 { | |
| 2337 if (scalar_non_zero) | |
| 2338 { | |
| 2339 tmp.xdata (kk) = scalar; | |
| 2340 tmp.xridx (kk++) = jj; | |
| 2341 } | |
| 2342 | |
| 2343 if (ii == jj && i < nz) | |
| 2344 if (++i < nz) | |
| 2345 ii = c_lhs.ridx(i); | |
| 2346 if (++j < n) | |
| 2347 jj = lhs_idx.elem(j); | |
| 2348 } | |
| 2349 } | |
| 2350 | |
| 2351 lhs = tmp; | |
| 2352 } | |
| 2353 else | |
| 2354 { | |
| 7246 | 2355 Sparse<LT> tmp (1, (max_idx > nc ? max_idx : nc), new_nzmx); |
| 5164 | 2356 |
| 5275 | 2357 octave_idx_type i = 0; |
| 2358 octave_idx_type ii = 0; | |
| 5164 | 2359 while (ii < nc && c_lhs.cidx(ii+1) <= i) |
| 2360 ii++; | |
| 2361 | |
| 5275 | 2362 octave_idx_type j = 0; |
| 2363 octave_idx_type jj = lhs_idx.elem(j); | |
| 2364 | |
| 2365 octave_idx_type kk = 0; | |
| 2366 octave_idx_type ic = 0; | |
| 5164 | 2367 |
| 2368 while (j < n || i < nz) | |
| 2369 { | |
|
11669
a6e08ecb4050
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7326
diff
changeset
|
2370 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
|
2371 { |
|
a6e08ecb4050
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7326
diff
changeset
|
2372 j++; |
|
a6e08ecb4050
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7326
diff
changeset
|
2373 jj = lhs_idx.elem (j); |
|
a6e08ecb4050
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7326
diff
changeset
|
2374 continue; |
|
a6e08ecb4050
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7326
diff
changeset
|
2375 } |
| 5164 | 2376 if (j == n || (i < nz && ii < jj)) |
| 2377 { | |
| 2378 while (ic <= ii) | |
| 2379 tmp.xcidx (ic++) = kk; | |
| 2380 tmp.xdata (kk) = c_lhs.data (i); | |
| 2381 i++; | |
| 2382 while (ii < nc && c_lhs.cidx(ii+1) <= i) | |
| 2383 ii++; | |
| 2384 } | |
| 2385 else | |
| 2386 { | |
| 2387 while (ic <= jj) | |
| 2388 tmp.xcidx (ic++) = kk; | |
| 2389 if (scalar_non_zero) | |
| 2390 tmp.xdata (kk) = scalar; | |
| 2391 if (ii == jj) | |
| 2392 { | |
| 2393 i++; | |
| 2394 while (ii < nc && c_lhs.cidx(ii+1) <= i) | |
| 2395 ii++; | |
| 2396 } | |
| 2397 j++; | |
| 2398 if (j < n) | |
| 2399 jj = lhs_idx.elem(j); | |
| 2400 } | |
| 2401 tmp.xridx (kk++) = 0; | |
| 2402 } | |
| 2403 | |
| 5275 | 2404 for (octave_idx_type iidx = ic; iidx < max_idx+1; iidx++) |
| 5164 | 2405 tmp.xcidx(iidx) = kk; |
| 2406 | |
| 2407 lhs = tmp; | |
| 2408 } | |
| 2409 } | |
| 2410 else | |
| 2411 { | |
| 2412 (*current_liboctave_error_handler) | |
| 2413 ("A(I) = X: X must be a scalar or a vector with same length as I"); | |
| 2414 | |
| 2415 retval = 0; | |
| 2416 } | |
| 2417 } | |
| 2418 else if (lhs_idx.is_colon ()) | |
| 2419 { | |
| 2420 if (lhs_len == 0) | |
| 2421 { | |
| 2422 | |
| 5681 | 2423 octave_idx_type new_nzmx = rhs.nnz (); |
| 5604 | 2424 Sparse<LT> tmp (1, rhs_len, new_nzmx); |
| 5164 | 2425 |
| 5275 | 2426 octave_idx_type ii = 0; |
| 2427 octave_idx_type jj = 0; | |
| 2428 for (octave_idx_type i = 0; i < rhs.cols(); i++) | |
| 2429 for (octave_idx_type j = rhs.cidx(i); j < rhs.cidx(i+1); j++) | |
| 5164 | 2430 { |
| 2431 OCTAVE_QUIT; | |
| 5275 | 2432 for (octave_idx_type k = jj; k <= i * rhs.rows() + rhs.ridx(j); k++) |
| 5164 | 2433 tmp.cidx(jj++) = ii; |
| 2434 | |
| 2435 tmp.data(ii) = rhs.data(j); | |
| 2436 tmp.ridx(ii++) = 0; | |
| 2437 } | |
| 2438 | |
| 5275 | 2439 for (octave_idx_type i = jj; i < rhs_len + 1; i++) |
| 5164 | 2440 tmp.cidx(i) = ii; |
| 2441 | |
| 2442 lhs = tmp; | |
| 2443 } | |
| 2444 else | |
| 2445 (*current_liboctave_error_handler) | |
| 2446 ("A(:) = X: A must be the same size as X"); | |
| 2447 } | |
| 2448 else if (! (rhs_len == 1 || rhs_len == 0)) | |
| 2449 { | |
| 2450 (*current_liboctave_error_handler) | |
| 2451 ("A([]) = X: X must also be an empty matrix or a scalar"); | |
| 2452 | |
| 2453 retval = 0; | |
| 2454 } | |
| 2455 | |
| 2456 lhs.clear_index (); | |
| 2457 | |
| 2458 return retval; | |
| 2459 } | |
| 2460 | |
| 2461 template <class LT, class RT> | |
| 2462 int | |
| 2463 assign (Sparse<LT>& lhs, const Sparse<RT>& rhs) | |
| 2464 { | |
| 2465 int retval = 1; | |
| 2466 | |
| 2467 int n_idx = lhs.index_count (); | |
| 2468 | |
| 5275 | 2469 octave_idx_type lhs_nr = lhs.rows (); |
| 2470 octave_idx_type lhs_nc = lhs.cols (); | |
| 5681 | 2471 octave_idx_type lhs_nz = lhs.nnz (); |
| 5275 | 2472 |
| 2473 octave_idx_type rhs_nr = rhs.rows (); | |
| 2474 octave_idx_type rhs_nc = rhs.cols (); | |
| 5164 | 2475 |
| 2476 idx_vector *tmp = lhs.get_idx (); | |
| 2477 | |
| 2478 idx_vector idx_i; | |
| 2479 idx_vector idx_j; | |
| 2480 | |
| 2481 if (n_idx > 2) | |
| 2482 { | |
| 2483 (*current_liboctave_error_handler) | |
| 2484 ("A(I, J) = X: can only have 1 or 2 indexes for sparse matrices"); | |
| 6092 | 2485 |
| 2486 lhs.clear_index (); | |
| 5164 | 2487 return 0; |
| 2488 } | |
| 2489 | |
| 2490 if (n_idx > 1) | |
| 2491 idx_j = tmp[1]; | |
| 2492 | |
| 2493 if (n_idx > 0) | |
| 2494 idx_i = tmp[0]; | |
| 2495 | |
| 6988 | 2496 // Take a constant copy of lhs. This means that ridx and family won't |
| 2497 // call make_unique. | |
| 2498 const Sparse<LT> c_lhs (lhs); | |
| 2499 | |
| 5164 | 2500 if (n_idx == 2) |
| 2501 { | |
| 5781 | 2502 octave_idx_type n = idx_i.freeze (lhs_nr, "row", true); |
| 2503 octave_idx_type m = idx_j.freeze (lhs_nc, "column", true); | |
| 5164 | 2504 |
| 2505 int idx_i_is_colon = idx_i.is_colon (); | |
| 2506 int idx_j_is_colon = idx_j.is_colon (); | |
| 2507 | |
| 7238 | 2508 if (lhs_nr == 0 && lhs_nc == 0) |
| 2509 { | |
| 2510 if (idx_i_is_colon) | |
| 2511 n = rhs_nr; | |
| 2512 | |
| 2513 if (idx_j_is_colon) | |
| 2514 m = rhs_nc; | |
| 2515 } | |
| 5164 | 2516 |
| 2517 if (idx_i && idx_j) | |
| 2518 { | |
| 2519 if (rhs_nr == 0 && rhs_nc == 0) | |
| 2520 { | |
| 2521 lhs.maybe_delete_elements (idx_i, idx_j); | |
| 2522 } | |
| 2523 else | |
| 2524 { | |
| 2525 if (rhs_nr == 1 && rhs_nc == 1 && n >= 0 && m >= 0) | |
| 2526 { | |
| 2527 if (n > 0 && m > 0) | |
| 2528 { | |
| 5603 | 2529 idx_i.sort (true); |
|
11669
a6e08ecb4050
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7326
diff
changeset
|
2530 n = idx_i.length (n); |
| 5603 | 2531 idx_j.sort (true); |
|
11669
a6e08ecb4050
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7326
diff
changeset
|
2532 m = idx_j.length (m); |
| 5603 | 2533 |
| 5275 | 2534 octave_idx_type max_row_idx = idx_i_is_colon ? rhs_nr : |
| 5164 | 2535 idx_i.max () + 1; |
| 5275 | 2536 octave_idx_type max_col_idx = idx_j_is_colon ? rhs_nc : |
| 5164 | 2537 idx_j.max () + 1; |
| 5603 | 2538 octave_idx_type new_nr = max_row_idx > lhs_nr ? |
| 2539 max_row_idx : lhs_nr; | |
| 2540 octave_idx_type new_nc = max_col_idx > lhs_nc ? | |
| 2541 max_col_idx : lhs_nc; | |
| 5164 | 2542 RT scalar = rhs.elem (0, 0); |
| 2543 | |
| 2544 // Count the number of non-zero terms | |
| 5681 | 2545 octave_idx_type new_nzmx = lhs.nnz (); |
| 5275 | 2546 for (octave_idx_type j = 0; j < m; j++) |
| 5164 | 2547 { |
| 5275 | 2548 octave_idx_type jj = idx_j.elem (j); |
| 5164 | 2549 if (jj < lhs_nc) |
| 2550 { | |
| 5275 | 2551 for (octave_idx_type i = 0; i < n; i++) |
| 5164 | 2552 { |
| 2553 OCTAVE_QUIT; | |
| 2554 | |
| 5275 | 2555 octave_idx_type ii = idx_i.elem (i); |
| 5164 | 2556 |
| 2557 if (ii < lhs_nr) | |
| 2558 { | |
| 6988 | 2559 for (octave_idx_type k = c_lhs.cidx(jj); |
| 2560 k < c_lhs.cidx(jj+1); k++) | |
| 5164 | 2561 { |
| 6988 | 2562 if (c_lhs.ridx(k) == ii) |
| 5604 | 2563 new_nzmx--; |
| 6988 | 2564 if (c_lhs.ridx(k) >= ii) |
| 5164 | 2565 break; |
| 2566 } | |
| 2567 } | |
| 2568 } | |
| 2569 } | |
| 2570 } | |
| 2571 | |
| 2572 if (scalar != RT()) | |
| 5604 | 2573 new_nzmx += m * n; |
| 2574 | |
| 2575 Sparse<LT> stmp (new_nr, new_nc, new_nzmx); | |
| 5164 | 2576 |
| 5275 | 2577 octave_idx_type jji = 0; |
| 2578 octave_idx_type jj = idx_j.elem (jji); | |
| 2579 octave_idx_type kk = 0; | |
| 5164 | 2580 stmp.cidx(0) = 0; |
| 5275 | 2581 for (octave_idx_type j = 0; j < new_nc; j++) |
| 5164 | 2582 { |
| 2583 if (jji < m && jj == j) | |
| 2584 { | |
| 5275 | 2585 octave_idx_type iii = 0; |
| 2586 octave_idx_type ii = idx_i.elem (iii); | |
| 5760 | 2587 octave_idx_type ppp = 0; |
| 6092 | 2588 octave_idx_type ppi = (j >= lhs_nc ? 0 : |
| 6988 | 2589 c_lhs.cidx(j+1) - |
| 2590 c_lhs.cidx(j)); | |
| 5760 | 2591 octave_idx_type pp = (ppp < ppi ? |
| 6988 | 2592 c_lhs.ridx(c_lhs.cidx(j)+ppp) : |
| 2593 new_nr); | |
| 5760 | 2594 while (ppp < ppi || iii < n) |
| 5164 | 2595 { |
| 5760 | 2596 if (iii < n && ii <= pp) |
| 5164 | 2597 { |
| 2598 if (scalar != RT ()) | |
| 2599 { | |
| 2600 stmp.data(kk) = scalar; | |
| 5760 | 2601 stmp.ridx(kk++) = ii; |
| 5164 | 2602 } |
| 5760 | 2603 if (ii == pp) |
| 6988 | 2604 pp = (++ppp < ppi ? c_lhs.ridx(c_lhs.cidx(j)+ppp) : new_nr); |
| 5164 | 2605 if (++iii < n) |
| 2606 ii = idx_i.elem(iii); | |
| 2607 } | |
| 5760 | 2608 else |
| 5164 | 2609 { |
| 5760 | 2610 stmp.data(kk) = |
| 6988 | 2611 c_lhs.data(c_lhs.cidx(j)+ppp); |
| 5760 | 2612 stmp.ridx(kk++) = pp; |
| 6988 | 2613 pp = (++ppp < ppi ? c_lhs.ridx(c_lhs.cidx(j)+ppp) : new_nr); |
| 5164 | 2614 } |
| 2615 } | |
| 2616 if (++jji < m) | |
| 2617 jj = idx_j.elem(jji); | |
| 2618 } | |
| 6988 | 2619 else if (j < lhs_nc) |
| 5164 | 2620 { |
| 6988 | 2621 for (octave_idx_type i = c_lhs.cidx(j); |
| 2622 i < c_lhs.cidx(j+1); i++) | |
| 5164 | 2623 { |
| 6988 | 2624 stmp.data(kk) = c_lhs.data(i); |
| 2625 stmp.ridx(kk++) = c_lhs.ridx(i); | |
| 5164 | 2626 } |
| 2627 } | |
| 2628 stmp.cidx(j+1) = kk; | |
| 2629 } | |
| 2630 | |
| 2631 lhs = stmp; | |
| 2632 } | |
| 7246 | 2633 else |
| 2634 { | |
| 7253 | 2635 #if 0 |
| 2636 // FIXME -- the following code will make this | |
| 2637 // function behave the same as the full matrix | |
| 2638 // case for things like | |
| 2639 // | |
| 2640 // x = sparse (ones (2)); | |
| 2641 // x([],3) = 2; | |
| 2642 // | |
| 2643 // x = | |
| 2644 // | |
| 2645 // Compressed Column Sparse (rows = 2, cols = 3, nnz = 4) | |
| 2646 // | |
| 2647 // (1, 1) -> 1 | |
| 2648 // (2, 1) -> 1 | |
| 2649 // (1, 2) -> 1 | |
| 2650 // (2, 2) -> 1 | |
| 2651 // | |
| 2652 // However, Matlab doesn't resize in this case | |
| 2653 // even though it does in the full matrix case. | |
| 2654 | |
| 7246 | 2655 if (n > 0) |
| 2656 { | |
| 2657 octave_idx_type max_row_idx = idx_i_is_colon ? | |
| 2658 rhs_nr : idx_i.max () + 1; | |
| 2659 octave_idx_type new_nr = max_row_idx > lhs_nr ? | |
| 2660 max_row_idx : lhs_nr; | |
| 2661 octave_idx_type new_nc = lhs_nc; | |
| 2662 | |
| 7253 | 2663 lhs.resize (new_nr, new_nc); |
| 7246 | 2664 } |
| 2665 else if (m > 0) | |
| 2666 { | |
| 2667 octave_idx_type max_col_idx = idx_j_is_colon ? | |
| 2668 rhs_nc : idx_j.max () + 1; | |
| 2669 octave_idx_type new_nr = lhs_nr; | |
| 2670 octave_idx_type new_nc = max_col_idx > lhs_nc ? | |
| 2671 max_col_idx : lhs_nc; | |
| 2672 | |
| 7253 | 2673 lhs.resize (new_nr, new_nc); |
| 7246 | 2674 } |
| 7253 | 2675 #endif |
| 7246 | 2676 } |
| 5164 | 2677 } |
| 2678 else if (n == rhs_nr && m == rhs_nc) | |
| 2679 { | |
| 2680 if (n > 0 && m > 0) | |
| 2681 { | |
| 5275 | 2682 octave_idx_type max_row_idx = idx_i_is_colon ? rhs_nr : |
| 5164 | 2683 idx_i.max () + 1; |
| 5275 | 2684 octave_idx_type max_col_idx = idx_j_is_colon ? rhs_nc : |
| 5164 | 2685 idx_j.max () + 1; |
| 5603 | 2686 octave_idx_type new_nr = max_row_idx > lhs_nr ? |
| 2687 max_row_idx : lhs_nr; | |
| 2688 octave_idx_type new_nc = max_col_idx > lhs_nc ? | |
| 2689 max_col_idx : lhs_nc; | |
| 2690 | |
| 2691 OCTAVE_LOCAL_BUFFER (octave_idx_type, rhs_idx_i, n); | |
| 2692 if (! idx_i.is_colon ()) | |
| 2693 { | |
| 2694 // Ok here we have to be careful with the indexing, | |
| 2695 // to treat cases like "a([3,2,1],:) = b", and still | |
| 2696 // handle the need for strict sorting of the sparse | |
| 2697 // elements. | |
| 2698 OCTAVE_LOCAL_BUFFER (octave_idx_vector_sort *, | |
| 2699 sidx, n); | |
| 2700 OCTAVE_LOCAL_BUFFER (octave_idx_vector_sort, | |
| 2701 sidxX, n); | |
| 2702 | |
| 2703 for (octave_idx_type i = 0; i < n; i++) | |
| 2704 { | |
| 2705 sidx[i] = &sidxX[i]; | |
| 2706 sidx[i]->i = idx_i.elem(i); | |
| 2707 sidx[i]->idx = i; | |
| 2708 } | |
| 2709 | |
| 2710 OCTAVE_QUIT; | |
| 2711 octave_sort<octave_idx_vector_sort *> | |
| 2712 sort (octave_idx_vector_comp); | |
| 2713 | |
| 2714 sort.sort (sidx, n); | |
| 2715 | |
| 2716 intNDArray<octave_idx_type> new_idx (dim_vector (n,1)); | |
| 2717 | |
| 2718 for (octave_idx_type i = 0; i < n; i++) | |
| 2719 { | |
| 2720 new_idx.xelem(i) = sidx[i]->i + 1; | |
| 2721 rhs_idx_i[i] = sidx[i]->idx; | |
| 2722 } | |
| 2723 | |
| 2724 idx_i = idx_vector (new_idx); | |
| 2725 } | |
| 2726 else | |
| 2727 for (octave_idx_type i = 0; i < n; i++) | |
| 2728 rhs_idx_i[i] = i; | |
| 2729 | |
| 2730 OCTAVE_LOCAL_BUFFER (octave_idx_type, rhs_idx_j, m); | |
| 2731 if (! idx_j.is_colon ()) | |
| 2732 { | |
| 2733 // Ok here we have to be careful with the indexing, | |
| 2734 // to treat cases like "a([3,2,1],:) = b", and still | |
| 2735 // handle the need for strict sorting of the sparse | |
| 2736 // elements. | |
| 2737 OCTAVE_LOCAL_BUFFER (octave_idx_vector_sort *, | |
| 2738 sidx, m); | |
| 2739 OCTAVE_LOCAL_BUFFER (octave_idx_vector_sort, | |
| 2740 sidxX, m); | |
| 2741 | |
| 2742 for (octave_idx_type i = 0; i < m; i++) | |
| 2743 { | |
| 2744 sidx[i] = &sidxX[i]; | |
| 2745 sidx[i]->i = idx_j.elem(i); | |
| 2746 sidx[i]->idx = i; | |
| 2747 } | |
| 2748 | |
| 2749 OCTAVE_QUIT; | |
| 2750 octave_sort<octave_idx_vector_sort *> | |
| 2751 sort (octave_idx_vector_comp); | |
| 2752 | |
| 2753 sort.sort (sidx, m); | |
| 2754 | |
| 2755 intNDArray<octave_idx_type> new_idx (dim_vector (m,1)); | |
| 2756 | |
| 2757 for (octave_idx_type i = 0; i < m; i++) | |
| 2758 { | |
| 2759 new_idx.xelem(i) = sidx[i]->i + 1; | |
| 2760 rhs_idx_j[i] = sidx[i]->idx; | |
| 2761 } | |
| 2762 | |
| 2763 idx_j = idx_vector (new_idx); | |
| 2764 } | |
| 2765 else | |
| 2766 for (octave_idx_type i = 0; i < m; i++) | |
| 2767 rhs_idx_j[i] = i; | |
| 5164 | 2768 |
| 5760 | 2769 // Maximum number of non-zero elements |
| 2770 octave_idx_type new_nzmx = lhs.nnz() + rhs.nnz(); | |
| 5164 | 2771 |
| 5604 | 2772 Sparse<LT> stmp (new_nr, new_nc, new_nzmx); |
| 5164 | 2773 |
| 5275 | 2774 octave_idx_type jji = 0; |
| 2775 octave_idx_type jj = idx_j.elem (jji); | |
| 2776 octave_idx_type kk = 0; | |
| 5164 | 2777 stmp.cidx(0) = 0; |
| 5275 | 2778 for (octave_idx_type j = 0; j < new_nc; j++) |
| 5164 | 2779 { |
| 2780 if (jji < m && jj == j) | |
| 2781 { | |
| 5275 | 2782 octave_idx_type iii = 0; |
| 2783 octave_idx_type ii = idx_i.elem (iii); | |
| 5760 | 2784 octave_idx_type ppp = 0; |
| 6092 | 2785 octave_idx_type ppi = (j >= lhs_nc ? 0 : |
| 6988 | 2786 c_lhs.cidx(j+1) - |
| 2787 c_lhs.cidx(j)); | |
| 5760 | 2788 octave_idx_type pp = (ppp < ppi ? |
| 6988 | 2789 c_lhs.ridx(c_lhs.cidx(j)+ppp) : |
| 2790 new_nr); | |
| 5760 | 2791 while (ppp < ppi || iii < n) |
| 5164 | 2792 { |
| 5760 | 2793 if (iii < n && ii <= pp) |
| 5164 | 2794 { |
|
11669
a6e08ecb4050
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7326
diff
changeset
|
2795 if (iii < n - 1 && |
|
a6e08ecb4050
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7326
diff
changeset
|
2796 idx_i.elem (iii + 1) == ii) |
|
a6e08ecb4050
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7326
diff
changeset
|
2797 { |
|
a6e08ecb4050
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7326
diff
changeset
|
2798 iii++; |
|
a6e08ecb4050
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7326
diff
changeset
|
2799 ii = idx_i.elem(iii); |
|
a6e08ecb4050
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7326
diff
changeset
|
2800 continue; |
|
a6e08ecb4050
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7326
diff
changeset
|
2801 } |
|
a6e08ecb4050
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7326
diff
changeset
|
2802 |
| 5603 | 2803 RT rtmp = rhs.elem (rhs_idx_i[iii], |
| 2804 rhs_idx_j[jji]); | |
| 5164 | 2805 if (rtmp != RT ()) |
| 2806 { | |
| 2807 stmp.data(kk) = rtmp; | |
| 5760 | 2808 stmp.ridx(kk++) = ii; |
| 5164 | 2809 } |
| 5760 | 2810 if (ii == pp) |
| 6988 | 2811 pp = (++ppp < ppi ? c_lhs.ridx(c_lhs.cidx(j)+ppp) : new_nr); |
| 5164 | 2812 if (++iii < n) |
| 2813 ii = idx_i.elem(iii); | |
| 2814 } | |
| 5760 | 2815 else |
| 5164 | 2816 { |
| 5760 | 2817 stmp.data(kk) = |
| 6988 | 2818 c_lhs.data(c_lhs.cidx(j)+ppp); |
| 5760 | 2819 stmp.ridx(kk++) = pp; |
| 6988 | 2820 pp = (++ppp < ppi ? c_lhs.ridx(c_lhs.cidx(j)+ppp) : new_nr); |
| 5164 | 2821 } |
| 2822 } | |
| 2823 if (++jji < m) | |
| 2824 jj = idx_j.elem(jji); | |
| 2825 } | |
| 6988 | 2826 else if (j < lhs_nc) |
| 5164 | 2827 { |
| 6988 | 2828 for (octave_idx_type i = c_lhs.cidx(j); |
| 2829 i < c_lhs.cidx(j+1); i++) | |
| 5164 | 2830 { |
| 6988 | 2831 stmp.data(kk) = c_lhs.data(i); |
| 2832 stmp.ridx(kk++) = c_lhs.ridx(i); | |
| 5164 | 2833 } |
| 2834 } | |
| 2835 stmp.cidx(j+1) = kk; | |
| 2836 } | |
| 2837 | |
| 5760 | 2838 stmp.maybe_compress(); |
| 5164 | 2839 lhs = stmp; |
| 2840 } | |
| 2841 } | |
| 2842 else if (n == 0 && m == 0) | |
| 2843 { | |
| 2844 if (! ((rhs_nr == 1 && rhs_nc == 1) | |
| 2845 || (rhs_nr == 0 || rhs_nc == 0))) | |
| 2846 { | |
| 2847 (*current_liboctave_error_handler) | |
| 2848 ("A([], []) = X: X must be an empty matrix or a scalar"); | |
| 2849 | |
| 2850 retval = 0; | |
| 2851 } | |
| 2852 } | |
| 2853 else | |
| 2854 { | |
| 2855 (*current_liboctave_error_handler) | |
| 2856 ("A(I, J) = X: X must be a scalar or the number of elements in I must"); | |
| 2857 (*current_liboctave_error_handler) | |
| 2858 ("match the number of rows in X and the number of elements in J must"); | |
| 2859 (*current_liboctave_error_handler) | |
| 2860 ("match the number of columns in X"); | |
| 2861 | |
| 2862 retval = 0; | |
| 2863 } | |
| 2864 } | |
| 2865 } | |
| 2866 // idx_vector::freeze() printed an error message for us. | |
| 2867 } | |
| 2868 else if (n_idx == 1) | |
| 2869 { | |
| 2870 int lhs_is_empty = lhs_nr == 0 || lhs_nc == 0; | |
| 2871 | |
| 2872 if (lhs_is_empty || (lhs_nr == 1 && lhs_nc == 1)) | |
| 2873 { | |
| 5275 | 2874 octave_idx_type lhs_len = lhs.length (); |
| 2875 | |
| 5781 | 2876 octave_idx_type n = idx_i.freeze (lhs_len, 0, true); |
| 5164 | 2877 |
| 2878 if (idx_i) | |
| 2879 { | |
| 2880 if (rhs_nr == 0 && rhs_nc == 0) | |
| 2881 { | |
| 2882 if (n != 0 && (lhs_nr != 0 || lhs_nc != 0)) | |
| 2883 lhs.maybe_delete_elements (idx_i); | |
| 2884 } | |
| 2885 else | |
| 2886 { | |
| 5781 | 2887 if (lhs_is_empty |
| 2888 && idx_i.is_colon () | |
| 2889 && ! (rhs_nr == 1 || rhs_nc == 1)) | |
| 5164 | 2890 { |
| 5781 | 2891 (*current_liboctave_warning_with_id_handler) |
| 2892 ("Octave:fortran-indexing", | |
| 2893 "A(:) = X: X is not a vector or scalar"); | |
| 2894 } | |
| 2895 else | |
| 2896 { | |
| 2897 octave_idx_type idx_nr = idx_i.orig_rows (); | |
| 2898 octave_idx_type idx_nc = idx_i.orig_columns (); | |
| 2899 | |
| 2900 if (! (rhs_nr == idx_nr && rhs_nc == idx_nc)) | |
| 2901 (*current_liboctave_warning_with_id_handler) | |
| 2902 ("Octave:fortran-indexing", | |
| 2903 "A(I) = X: X does not have same shape as I"); | |
| 5164 | 2904 } |
| 2905 | |
| 5760 | 2906 if (! assign1 (lhs, rhs)) |
| 5164 | 2907 retval = 0; |
| 2908 } | |
| 2909 } | |
| 2910 // idx_vector::freeze() printed an error message for us. | |
| 2911 } | |
| 2912 else if (lhs_nr == 1) | |
| 2913 { | |
| 5781 | 2914 idx_i.freeze (lhs_nc, "vector", true); |
| 5164 | 2915 |
| 2916 if (idx_i) | |
| 2917 { | |
| 2918 if (rhs_nr == 0 && rhs_nc == 0) | |
| 2919 lhs.maybe_delete_elements (idx_i); | |
| 5760 | 2920 else if (! assign1 (lhs, rhs)) |
| 5164 | 2921 retval = 0; |
| 2922 } | |
| 2923 // idx_vector::freeze() printed an error message for us. | |
| 2924 } | |
| 2925 else if (lhs_nc == 1) | |
| 2926 { | |
| 5781 | 2927 idx_i.freeze (lhs_nr, "vector", true); |
| 5164 | 2928 |
| 2929 if (idx_i) | |
| 2930 { | |
| 2931 if (rhs_nr == 0 && rhs_nc == 0) | |
| 2932 lhs.maybe_delete_elements (idx_i); | |
| 5760 | 2933 else if (! assign1 (lhs, rhs)) |
| 5164 | 2934 retval = 0; |
| 2935 } | |
| 2936 // idx_vector::freeze() printed an error message for us. | |
| 2937 } | |
| 2938 else | |
| 2939 { | |
| 5781 | 2940 if (! (idx_i.is_colon () |
| 2941 || (idx_i.one_zero_only () | |
| 2942 && idx_i.orig_rows () == lhs_nr | |
| 2943 && idx_i.orig_columns () == lhs_nc))) | |
| 2944 (*current_liboctave_warning_with_id_handler) | |
| 2945 ("Octave:fortran-indexing", "single index used for matrix"); | |
| 5164 | 2946 |
| 5275 | 2947 octave_idx_type lhs_len = lhs.length (); |
| 2948 | |
| 2949 octave_idx_type len = idx_i.freeze (lhs_nr * lhs_nc, "matrix"); | |
| 5164 | 2950 |
| 2951 if (idx_i) | |
| 2952 { | |
| 2953 if (rhs_nr == 0 && rhs_nc == 0) | |
| 2954 lhs.maybe_delete_elements (idx_i); | |
| 2955 else if (len == 0) | |
| 2956 { | |
| 2957 if (! ((rhs_nr == 1 && rhs_nc == 1) | |
| 2958 || (rhs_nr == 0 || rhs_nc == 0))) | |
| 2959 (*current_liboctave_error_handler) | |
| 2960 ("A([]) = X: X must be an empty matrix or scalar"); | |
| 2961 } | |
| 2962 else if (len == rhs_nr * rhs_nc) | |
| 2963 { | |
| 5604 | 2964 octave_idx_type new_nzmx = lhs_nz; |
| 5603 | 2965 OCTAVE_LOCAL_BUFFER (octave_idx_type, rhs_idx, len); |
| 2966 | |
| 2967 if (! idx_i.is_colon ()) | |
| 2968 { | |
| 2969 // Ok here we have to be careful with the indexing, to | |
| 2970 // treat cases like "a([3,2,1]) = b", and still handle | |
| 2971 // the need for strict sorting of the sparse elements. | |
| 2972 | |
| 2973 OCTAVE_LOCAL_BUFFER (octave_idx_vector_sort *, sidx, | |
| 2974 len); | |
| 2975 OCTAVE_LOCAL_BUFFER (octave_idx_vector_sort, sidxX, | |
| 2976 len); | |
| 2977 | |
| 2978 for (octave_idx_type i = 0; i < len; i++) | |
| 2979 { | |
| 2980 sidx[i] = &sidxX[i]; | |
| 2981 sidx[i]->i = idx_i.elem(i); | |
| 2982 sidx[i]->idx = i; | |
| 2983 } | |
| 2984 | |
| 2985 OCTAVE_QUIT; | |
| 2986 octave_sort<octave_idx_vector_sort *> | |
| 2987 sort (octave_idx_vector_comp); | |
| 2988 | |
| 2989 sort.sort (sidx, len); | |
| 2990 | |
| 2991 intNDArray<octave_idx_type> new_idx (dim_vector (len,1)); | |
| 2992 | |
| 2993 for (octave_idx_type i = 0; i < len; i++) | |
| 2994 { | |
| 2995 new_idx.xelem(i) = sidx[i]->i + 1; | |
| 2996 rhs_idx[i] = sidx[i]->idx; | |
| 2997 } | |
| 2998 | |
| 2999 idx_i = idx_vector (new_idx); | |
| 3000 } | |
| 3001 else | |
| 3002 for (octave_idx_type i = 0; i < len; i++) | |
| 3003 rhs_idx[i] = i; | |
| 5164 | 3004 |
| 3005 // First count the number of non-zero elements | |
| 5275 | 3006 for (octave_idx_type i = 0; i < len; i++) |
| 5164 | 3007 { |
| 3008 OCTAVE_QUIT; | |
| 3009 | |
| 5275 | 3010 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
|
3011 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
|
3012 continue; |
| 5164 | 3013 if (ii < lhs_len && c_lhs.elem(ii) != LT ()) |
| 5604 | 3014 new_nzmx--; |
| 5603 | 3015 if (rhs.elem(rhs_idx[i]) != RT ()) |
| 5604 | 3016 new_nzmx++; |
| 5164 | 3017 } |
| 3018 | |
| 5604 | 3019 Sparse<LT> stmp (lhs_nr, lhs_nc, new_nzmx); |
| 5164 | 3020 |
| 5275 | 3021 octave_idx_type i = 0; |
| 3022 octave_idx_type ii = 0; | |
| 3023 octave_idx_type ic = 0; | |
| 5164 | 3024 if (i < lhs_nz) |
| 3025 { | |
| 3026 while (ic < lhs_nc && i >= c_lhs.cidx(ic+1)) | |
| 3027 ic++; | |
| 3028 ii = ic * lhs_nr + c_lhs.ridx(i); | |
| 3029 } | |
| 3030 | |
| 5275 | 3031 octave_idx_type j = 0; |
| 3032 octave_idx_type jj = idx_i.elem (j); | |
| 3033 octave_idx_type jr = jj % lhs_nr; | |
| 3034 octave_idx_type jc = (jj - jr) / lhs_nr; | |
| 3035 | |
| 3036 octave_idx_type kk = 0; | |
| 3037 octave_idx_type kc = 0; | |
| 5164 | 3038 |
| 3039 while (j < len || i < lhs_nz) | |
| 3040 { | |
| 3041 if (j == len || (i < lhs_nz && ii < jj)) | |
| 3042 { | |
| 3043 while (kc <= ic) | |
| 3044 stmp.xcidx (kc++) = kk; | |
| 3045 stmp.xdata (kk) = c_lhs.data (i); | |
| 3046 stmp.xridx (kk++) = c_lhs.ridx (i); | |
| 3047 i++; | |
| 3048 while (ic < lhs_nc && i >= c_lhs.cidx(ic+1)) | |
| 3049 ic++; | |
| 3050 if (i < lhs_nz) | |
| 3051 ii = ic * lhs_nr + c_lhs.ridx(i); | |
| 3052 } | |
| 3053 else | |
| 3054 { | |
| 3055 while (kc <= jc) | |
| 3056 stmp.xcidx (kc++) = kk; | |
| 5603 | 3057 RT rtmp = rhs.elem (rhs_idx[j]); |
| 5164 | 3058 if (rtmp != RT ()) |
| 3059 { | |
| 3060 stmp.xdata (kk) = rtmp; | |
| 3061 stmp.xridx (kk++) = jr; | |
| 3062 } | |
| 3063 if (ii == jj) | |
| 3064 { | |
| 3065 i++; | |
| 3066 while (ic < lhs_nc && i >= c_lhs.cidx(ic+1)) | |
| 3067 ic++; | |
| 3068 if (i < lhs_nz) | |
| 3069 ii = ic * lhs_nr + c_lhs.ridx(i); | |
| 3070 } | |
| 3071 j++; | |
| 3072 if (j < len) | |
| 3073 { | |
| 3074 jj = idx_i.elem (j); | |
| 3075 jr = jj % lhs_nr; | |
| 3076 jc = (jj - jr) / lhs_nr; | |
| 3077 } | |
| 3078 } | |
| 3079 } | |
| 3080 | |
| 5275 | 3081 for (octave_idx_type iidx = kc; iidx < lhs_nc+1; iidx++) |
| 5603 | 3082 stmp.xcidx(iidx) = kk; |
| 5164 | 3083 |
| 3084 lhs = stmp; | |
| 3085 } | |
| 3086 else if (rhs_nr == 1 && rhs_nc == 1) | |
| 3087 { | |
| 3088 RT scalar = rhs.elem (0, 0); | |
| 5604 | 3089 octave_idx_type new_nzmx = lhs_nz; |
| 5603 | 3090 idx_i.sort (true); |
|
11669
a6e08ecb4050
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7326
diff
changeset
|
3091 len = idx_i.length (len); |
| 5164 | 3092 |
| 3093 // First count the number of non-zero elements | |
| 3094 if (scalar != RT ()) | |
| 5604 | 3095 new_nzmx += len; |
| 5275 | 3096 for (octave_idx_type i = 0; i < len; i++) |
| 5164 | 3097 { |
| 3098 OCTAVE_QUIT; | |
| 5275 | 3099 octave_idx_type ii = idx_i.elem (i); |
| 5164 | 3100 if (ii < lhs_len && c_lhs.elem(ii) != LT ()) |
| 5604 | 3101 new_nzmx--; |
| 5164 | 3102 } |
| 3103 | |
| 5604 | 3104 Sparse<LT> stmp (lhs_nr, lhs_nc, new_nzmx); |
| 5164 | 3105 |
| 5275 | 3106 octave_idx_type i = 0; |
| 3107 octave_idx_type ii = 0; | |
| 3108 octave_idx_type ic = 0; | |
| 5164 | 3109 if (i < lhs_nz) |
| 3110 { | |
| 3111 while (ic < lhs_nc && i >= c_lhs.cidx(ic+1)) | |
| 3112 ic++; | |
| 3113 ii = ic * lhs_nr + c_lhs.ridx(i); | |
| 3114 } | |
| 3115 | |
| 5275 | 3116 octave_idx_type j = 0; |
| 3117 octave_idx_type jj = idx_i.elem (j); | |
| 3118 octave_idx_type jr = jj % lhs_nr; | |
| 3119 octave_idx_type jc = (jj - jr) / lhs_nr; | |
| 3120 | |
| 3121 octave_idx_type kk = 0; | |
| 3122 octave_idx_type kc = 0; | |
| 5164 | 3123 |
| 3124 while (j < len || i < lhs_nz) | |
| 3125 { | |
|
11669
a6e08ecb4050
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7326
diff
changeset
|
3126 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
|
3127 { |
|
a6e08ecb4050
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7326
diff
changeset
|
3128 j++; |
|
a6e08ecb4050
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7326
diff
changeset
|
3129 jj = idx_i.elem (j); |
|
a6e08ecb4050
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7326
diff
changeset
|
3130 jr = jj % lhs_nr; |
|
a6e08ecb4050
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7326
diff
changeset
|
3131 jc = (jj - jr) / lhs_nr; |
|
a6e08ecb4050
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7326
diff
changeset
|
3132 continue; |
|
a6e08ecb4050
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7326
diff
changeset
|
3133 } |
|
a6e08ecb4050
Treat repeated indices in the sparse assignments
David Bateman <dbateman@free.fr>
parents:
7326
diff
changeset
|
3134 |
| 5164 | 3135 if (j == len || (i < lhs_nz && ii < jj)) |
| 3136 { | |
| 3137 while (kc <= ic) | |
| 3138 stmp.xcidx (kc++) = kk; | |
| 3139 stmp.xdata (kk) = c_lhs.data (i); | |
| 3140 stmp.xridx (kk++) = c_lhs.ridx (i); | |
| 3141 i++; | |
| 3142 while (ic < lhs_nc && i >= c_lhs.cidx(ic+1)) | |
| 3143 ic++; | |
| 3144 if (i < lhs_nz) | |
| 3145 ii = ic * lhs_nr + c_lhs.ridx(i); | |
| 3146 } | |
| 3147 else | |
| 3148 { | |
| 3149 while (kc <= jc) | |
| 3150 stmp.xcidx (kc++) = kk; | |
| 3151 if (scalar != RT ()) | |
| 3152 { | |
| 3153 stmp.xdata (kk) = scalar; | |
| 3154 stmp.xridx (kk++) = jr; | |
| 3155 } | |
| 3156 if (ii == jj) | |
| 3157 { | |
| 3158 i++; | |
| 3159 while (ic < lhs_nc && i >= c_lhs.cidx(ic+1)) | |
| 3160 ic++; | |
| 3161 if (i < lhs_nz) | |
| 3162 ii = ic * lhs_nr + c_lhs.ridx(i); | |
| 3163 } | |
| 3164 j++; | |
| 3165 if (j < len) | |
| 3166 { | |
| 3167 jj = idx_i.elem (j); | |
| 3168 jr = jj % lhs_nr; | |
| 3169 jc = (jj - jr) / lhs_nr; | |
| 3170 } | |
| 3171 } | |
| 3172 } | |
| 3173 | |
| 5275 | 3174 for (octave_idx_type iidx = kc; iidx < lhs_nc+1; iidx++) |
| 5164 | 3175 stmp.xcidx(iidx) = kk; |
| 3176 | |
| 3177 lhs = stmp; | |
| 3178 } | |
| 3179 else | |
| 3180 { | |
| 3181 (*current_liboctave_error_handler) | |
| 3182 ("A(I) = X: X must be a scalar or a matrix with the same size as I"); | |
| 3183 | |
| 3184 retval = 0; | |
| 3185 } | |
| 3186 } | |
| 3187 // idx_vector::freeze() printed an error message for us. | |
| 3188 } | |
| 3189 } | |
| 3190 else | |
| 3191 { | |
| 3192 (*current_liboctave_error_handler) | |
| 3193 ("invalid number of indices for matrix expression"); | |
| 3194 | |
| 3195 retval = 0; | |
| 3196 } | |
| 3197 | |
| 3198 lhs.clear_index (); | |
| 3199 | |
| 3200 return retval; | |
| 3201 } | |
| 3202 | |
| 3203 template <class T> | |
| 3204 void | |
| 3205 Sparse<T>::print_info (std::ostream& os, const std::string& prefix) const | |
| 3206 { | |
| 3207 os << prefix << "rep address: " << rep << "\n" | |
| 5604 | 3208 << prefix << "rep->nzmx: " << rep->nzmx << "\n" |
| 5164 | 3209 << prefix << "rep->nrows: " << rep->nrows << "\n" |
| 3210 << prefix << "rep->ncols: " << rep->ncols << "\n" | |
| 3211 << prefix << "rep->data: " << static_cast<void *> (rep->d) << "\n" | |
| 3212 << prefix << "rep->ridx: " << static_cast<void *> (rep->r) << "\n" | |
| 3213 << prefix << "rep->cidx: " << static_cast<void *> (rep->c) << "\n" | |
| 3214 << prefix << "rep->count: " << rep->count << "\n"; | |
| 3215 } | |
| 3216 | |
| 3217 /* | |
| 3218 ;;; Local Variables: *** | |
| 3219 ;;; mode: C++ *** | |
| 3220 ;;; End: *** | |
| 3221 */ |