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