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