Mercurial > octave

comparison liboctave/array/Sparse.cc @ 17769:49a5a4be04a1

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
maint: Use GNU style coding conventions for code in liboctave/ * liboctave/array/Array-C.cc, liboctave/array/Array-b.cc, liboctave/array/Array-ch.cc, liboctave/array/Array-d.cc, liboctave/array/Array-f.cc, liboctave/array/Array-fC.cc, liboctave/array/Array-util.cc, liboctave/array/Array-util.h, liboctave/array/Array.cc, liboctave/array/Array.h, liboctave/array/Array3.h, liboctave/array/CColVector.cc, liboctave/array/CColVector.h, liboctave/array/CDiagMatrix.cc, liboctave/array/CDiagMatrix.h, liboctave/array/CMatrix.cc, liboctave/array/CMatrix.h, liboctave/array/CNDArray.cc, liboctave/array/CNDArray.h, liboctave/array/CRowVector.cc, liboctave/array/CRowVector.h, liboctave/array/CSparse.cc, liboctave/array/CSparse.h, liboctave/array/DiagArray2.h, liboctave/array/MArray.cc, liboctave/array/MArray.h, liboctave/array/MDiagArray2.cc, liboctave/array/MDiagArray2.h, liboctave/array/MSparse.cc, liboctave/array/MSparse.h, liboctave/array/MatrixType.cc, liboctave/array/MatrixType.h, liboctave/array/PermMatrix.h, liboctave/array/Range.cc, liboctave/array/Range.h, liboctave/array/Sparse.cc, liboctave/array/Sparse.h, liboctave/array/boolMatrix.cc, liboctave/array/boolMatrix.h, liboctave/array/boolNDArray.cc, liboctave/array/boolNDArray.h, liboctave/array/boolSparse.cc, liboctave/array/boolSparse.h, liboctave/array/chMatrix.cc, liboctave/array/chMatrix.h, liboctave/array/chNDArray.cc, liboctave/array/chNDArray.h, liboctave/array/dColVector.h, liboctave/array/dDiagMatrix.cc, liboctave/array/dDiagMatrix.h, liboctave/array/dMatrix.cc, liboctave/array/dMatrix.h, liboctave/array/dNDArray.cc, liboctave/array/dNDArray.h, liboctave/array/dRowVector.h, liboctave/array/dSparse.cc, liboctave/array/dSparse.h, liboctave/array/dim-vector.cc, liboctave/array/dim-vector.h, liboctave/array/fCColVector.cc, liboctave/array/fCColVector.h, liboctave/array/fCDiagMatrix.cc, liboctave/array/fCDiagMatrix.h, liboctave/array/fCMatrix.cc, liboctave/array/fCMatrix.h, liboctave/array/fCNDArray.cc, liboctave/array/fCNDArray.h, liboctave/array/fCRowVector.cc, liboctave/array/fCRowVector.h, liboctave/array/fColVector.h, liboctave/array/fDiagMatrix.cc, liboctave/array/fDiagMatrix.h, liboctave/array/fMatrix.cc, liboctave/array/fMatrix.h, liboctave/array/fNDArray.cc, liboctave/array/fNDArray.h, liboctave/array/fRowVector.h, liboctave/array/idx-vector.cc, liboctave/array/idx-vector.h, liboctave/array/intNDArray.cc, liboctave/array/intNDArray.h, liboctave/cruft/misc/blaswrap.c, liboctave/cruft/misc/quit.cc, liboctave/numeric/CmplxCHOL.cc, liboctave/numeric/CmplxCHOL.h, liboctave/numeric/CmplxGEPBAL.cc, liboctave/numeric/CmplxGEPBAL.h, liboctave/numeric/CmplxHESS.h, liboctave/numeric/CmplxLU.cc, liboctave/numeric/CmplxLU.h, liboctave/numeric/CmplxQR.cc, liboctave/numeric/CmplxQRP.cc, liboctave/numeric/CmplxQRP.h, liboctave/numeric/CmplxSCHUR.h, liboctave/numeric/CmplxSVD.cc, liboctave/numeric/CmplxSVD.h, liboctave/numeric/CollocWt.h, liboctave/numeric/DAE.h, liboctave/numeric/DAEFunc.h, liboctave/numeric/DAERT.h, liboctave/numeric/DAERTFunc.h, liboctave/numeric/DASPK.cc, liboctave/numeric/DASRT.cc, liboctave/numeric/DASRT.h, liboctave/numeric/DASSL.cc, liboctave/numeric/DET.h, liboctave/numeric/EIG.cc, liboctave/numeric/EIG.h, liboctave/numeric/LSODE.cc, liboctave/numeric/ODE.h, liboctave/numeric/ODEFunc.h, liboctave/numeric/ODES.h, liboctave/numeric/ODESFunc.h, liboctave/numeric/Quad.cc, liboctave/numeric/Quad.h, liboctave/numeric/SparseCmplxCHOL.h, liboctave/numeric/SparseCmplxLU.cc, liboctave/numeric/SparseCmplxLU.h, liboctave/numeric/SparseCmplxQR.cc, liboctave/numeric/SparseCmplxQR.h, liboctave/numeric/SparseQR.cc, liboctave/numeric/SparseQR.h, liboctave/numeric/SparsedbleCHOL.h, liboctave/numeric/SparsedbleLU.cc, liboctave/numeric/SparsedbleLU.h, liboctave/numeric/base-aepbal.h, liboctave/numeric/base-dae.h, liboctave/numeric/base-de.h, liboctave/numeric/base-lu.cc, liboctave/numeric/base-lu.h, liboctave/numeric/base-min.h, liboctave/numeric/base-qr.h, liboctave/numeric/bsxfun.h, liboctave/numeric/dbleCHOL.cc, liboctave/numeric/dbleCHOL.h, liboctave/numeric/dbleGEPBAL.h, liboctave/numeric/dbleHESS.h, liboctave/numeric/dbleLU.cc, liboctave/numeric/dbleLU.h, liboctave/numeric/dbleQR.cc, liboctave/numeric/dbleQRP.cc, liboctave/numeric/dbleQRP.h, liboctave/numeric/dbleSCHUR.cc, liboctave/numeric/dbleSCHUR.h, liboctave/numeric/dbleSVD.cc, liboctave/numeric/dbleSVD.h, liboctave/numeric/eigs-base.cc, liboctave/numeric/fCmplxAEPBAL.cc, liboctave/numeric/fCmplxAEPBAL.h, liboctave/numeric/fCmplxCHOL.cc, liboctave/numeric/fCmplxCHOL.h, liboctave/numeric/fCmplxGEPBAL.cc, liboctave/numeric/fCmplxGEPBAL.h, liboctave/numeric/fCmplxHESS.h, liboctave/numeric/fCmplxLU.cc, liboctave/numeric/fCmplxLU.h, liboctave/numeric/fCmplxQR.cc, liboctave/numeric/fCmplxQR.h, liboctave/numeric/fCmplxQRP.cc, liboctave/numeric/fCmplxQRP.h, liboctave/numeric/fCmplxSCHUR.cc, liboctave/numeric/fCmplxSCHUR.h, liboctave/numeric/fCmplxSVD.h, liboctave/numeric/fEIG.cc, liboctave/numeric/fEIG.h, liboctave/numeric/floatCHOL.cc, liboctave/numeric/floatCHOL.h, liboctave/numeric/floatGEPBAL.cc, liboctave/numeric/floatGEPBAL.h, liboctave/numeric/floatHESS.h, liboctave/numeric/floatLU.cc, liboctave/numeric/floatLU.h, liboctave/numeric/floatQR.cc, liboctave/numeric/floatQRP.cc, liboctave/numeric/floatQRP.h, liboctave/numeric/floatSCHUR.cc, liboctave/numeric/floatSCHUR.h, liboctave/numeric/floatSVD.cc, liboctave/numeric/floatSVD.h, liboctave/numeric/lo-mappers.cc, liboctave/numeric/lo-mappers.h, liboctave/numeric/lo-specfun.cc, liboctave/numeric/lo-specfun.h, liboctave/numeric/oct-convn.cc, liboctave/numeric/oct-fftw.cc, liboctave/numeric/oct-fftw.h, liboctave/numeric/oct-norm.cc, liboctave/numeric/oct-rand.cc, liboctave/numeric/oct-rand.h, liboctave/numeric/randgamma.c, liboctave/numeric/randgamma.h, liboctave/numeric/randmtzig.c, liboctave/numeric/randpoisson.c, liboctave/numeric/randpoisson.h, liboctave/numeric/sparse-base-chol.h, liboctave/numeric/sparse-base-lu.h, liboctave/numeric/sparse-dmsolve.cc, liboctave/operators/Sparse-diag-op-defs.h, liboctave/operators/Sparse-op-defs.h, liboctave/operators/mx-inlines.cc, liboctave/system/dir-ops.h, liboctave/system/file-ops.cc, liboctave/system/file-stat.cc, liboctave/system/file-stat.h, liboctave/system/lo-sysdep.cc, liboctave/system/lo-sysdep.h, liboctave/system/mach-info.cc, liboctave/system/mach-info.h, liboctave/system/oct-env.cc, liboctave/system/oct-group.cc, liboctave/system/oct-syscalls.cc, liboctave/system/oct-syscalls.h, liboctave/system/oct-time.h, liboctave/system/tempname.c, liboctave/util/action-container.h, liboctave/util/base-list.h, liboctave/util/cmd-edit.cc, liboctave/util/cmd-edit.h, liboctave/util/cmd-hist.cc, liboctave/util/cmd-hist.h, liboctave/util/data-conv.cc, liboctave/util/data-conv.h, liboctave/util/kpse.cc, liboctave/util/lo-array-gripes.cc, liboctave/util/lo-cieee.c, liboctave/util/lo-regexp.cc, liboctave/util/lo-utils.cc, liboctave/util/oct-alloc.cc, liboctave/util/oct-base64.cc, liboctave/util/oct-binmap.h, liboctave/util/oct-cmplx.h, liboctave/util/oct-glob.cc, liboctave/util/oct-inttypes.cc, liboctave/util/oct-inttypes.h, liboctave/util/oct-locbuf.cc, liboctave/util/oct-locbuf.h, liboctave/util/oct-mem.h, liboctave/util/oct-mutex.cc, liboctave/util/oct-refcount.h, liboctave/util/oct-shlib.cc, liboctave/util/oct-shlib.h, liboctave/util/oct-sort.cc, liboctave/util/oct-sort.h, liboctave/util/pathsearch.cc, liboctave/util/pathsearch.h, liboctave/util/sparse-util.cc, liboctave/util/str-vec.cc, liboctave/util/str-vec.h, liboctave/util/unwind-prot.h, liboctave/util/url-transfer.cc, liboctave/util/url-transfer.h: Use GNU style coding conventions.
author Rik <rik@octave.org>
date Sat, 26 Oct 2013 18:57:05 -0700
parents d63878346099
children eaf5c3ef3e8d
comparison
equal deleted inserted replaced
17768:271c0cce0f64 17769:49a5a4be04a1
243 : rep (0), dimensions (dv) 243 : rep (0), dimensions (dv)
244 { 244 {
245 245
246 // Work in unsigned long long to avoid overflow issues with numel 246 // Work in unsigned long long to avoid overflow issues with numel
247 unsigned long long a_nel = static_cast<unsigned long long>(a.rows ()) * 247 unsigned long long a_nel = static_cast<unsigned long long>(a.rows ()) *
248 static_cast<unsigned long long>(a.cols ()); 248 static_cast<unsigned long long>(a.cols ());
249 unsigned long long dv_nel = static_cast<unsigned long long>(dv (0)) * 249 unsigned long long dv_nel = static_cast<unsigned long long>(dv (0)) *
250 static_cast<unsigned long long>(dv (1)); 250 static_cast<unsigned long long>(dv (1));
251 251
252 if (a_nel != dv_nel) 252 if (a_nel != dv_nel)
253 (*current_liboctave_error_handler) 253 (*current_liboctave_error_handler)
254 ("Sparse::Sparse (const Sparse&, const dim_vector&): dimension mismatch"); 254 ("Sparse::Sparse (const Sparse&, const dim_vector&): dimension mismatch");
255 else 255 else
735 return foo; 735 return foo;
736 } 736 }
737 737
738 template <class T> 738 template <class T>
739 T 739 T
740 Sparse<T>::range_error (const char *fcn, octave_idx_type i, octave_idx_type j) const 740 Sparse<T>::range_error (const char *fcn, octave_idx_type i,
741 octave_idx_type j) const
741 { 742 {
742 (*current_liboctave_error_handler) 743 (*current_liboctave_error_handler)
743 ("%s (%d, %d): range error", fcn, i, j); 744 ("%s (%d, %d): range error", fcn, i, j);
744 return T (); 745 return T ();
745 } 746 }
754 return foo; 755 return foo;
755 } 756 }
756 757
757 template <class T> 758 template <class T>
758 T 759 T
759 Sparse<T>::range_error (const char *fcn, const Array<octave_idx_type>& ra_idx) const 760 Sparse<T>::range_error (const char *fcn,
761 const Array<octave_idx_type>& ra_idx) const
760 { 762 {
761 std::ostringstream buf; 763 std::ostringstream buf;
762 764
763 buf << fcn << " ("; 765 buf << fcn << " (";
764 766
970 copy_or_memcpy (std::min (c, rep->ncols)+1, rep->c, new_cidx); 972 copy_or_memcpy (std::min (c, rep->ncols)+1, rep->c, new_cidx);
971 delete [] rep->c; 973 delete [] rep->c;
972 rep->c = new_cidx; 974 rep->c = new_cidx;
973 975
974 if (c > rep->ncols) 976 if (c > rep->ncols)
975 fill_or_memset (c - rep->ncols, rep->c[rep->ncols], rep->c + rep->ncols + 1); 977 fill_or_memset (c - rep->ncols, rep->c[rep->ncols],
978 rep->c + rep->ncols + 1);
976 } 979 }
977 980
978 rep->ncols = dimensions(1) = c; 981 rep->ncols = dimensions(1) = c;
979 982
980 rep->change_length (rep->nnz ()); 983 rep->change_length (rep->nnz ());
1182 const octave_idx_type *sj = sidx.raw (); 1185 const octave_idx_type *sj = sidx.raw ();
1183 octave_idx_type sl = sidx.length (nel), nz_new = 0, j = 0; 1186 octave_idx_type sl = sidx.length (nel), nz_new = 0, j = 0;
1184 for (octave_idx_type i = 0; i < nz; i++) 1187 for (octave_idx_type i = 0; i < nz; i++)
1185 { 1188 {
1186 octave_idx_type r = tmp.ridx (i); 1189 octave_idx_type r = tmp.ridx (i);
1187 for (;j < sl && sj[j] < r; j++) ; 1190 for (; j < sl && sj[j] < r; j++) ;
1188 if (j == sl || sj[j] > r) 1191 if (j == sl || sj[j] > r)
1189 { 1192 {
1190 data_new[nz_new] = tmp.data (i); 1193 data_new[nz_new] = tmp.data (i);
1191 ridx_new[nz_new++] = r - j; 1194 ridx_new[nz_new++] = r - j;
1192 } 1195 }
1203 // Sparse row vector. 1206 // Sparse row vector.
1204 octave_idx_type lb, ub; 1207 octave_idx_type lb, ub;
1205 if (idx.is_cont_range (nc, lb, ub)) 1208 if (idx.is_cont_range (nc, lb, ub))
1206 { 1209 {
1207 const Sparse<T> tmp = *this; 1210 const Sparse<T> tmp = *this;
1208 octave_idx_type lbi = tmp.cidx (lb), ubi = tmp.cidx (ub), new_nz = nz - (ubi - lbi); 1211 octave_idx_type lbi = tmp.cidx (lb);
1212 octave_idx_type ubi = tmp.cidx (ub);
1213 octave_idx_type new_nz = nz - (ubi - lbi);
1209 *this = Sparse<T> (1, nc - (ub - lb), new_nz); 1214 *this = Sparse<T> (1, nc - (ub - lb), new_nz);
1210 copy_or_memcpy (lbi, tmp.data (), data ()); 1215 copy_or_memcpy (lbi, tmp.data (), data ());
1211 copy_or_memcpy (nz - ubi, tmp.data () + ubi, xdata () + lbi); 1216 copy_or_memcpy (nz - ubi, tmp.data () + ubi, xdata () + lbi);
1212 fill_or_memset (new_nz, static_cast<octave_idx_type> (0), ridx ()); 1217 fill_or_memset (new_nz, static_cast<octave_idx_type> (0), ridx ());
1213 copy_or_memcpy (lb, tmp.cidx () + 1, cidx () + 1); 1218 copy_or_memcpy (lb, tmp.cidx () + 1, cidx () + 1);
1214 mx_inline_sub (nc - ub, xcidx () + 1, tmp.cidx () + ub + 1, ubi - lbi); 1219 mx_inline_sub (nc - ub, xcidx () + 1, tmp.cidx () + ub + 1,
1220 ubi - lbi);
1215 } 1221 }
1216 else 1222 else
1217 *this = index (idx.complement (nc)); 1223 *this = index (idx.complement (nc));
1218 } 1224 }
1219 else if (idx.length (nel) != 0) 1225 else if (idx.length (nel) != 0)
1258 *this = Sparse<T> (nr, nc - (ub - lb)); 1264 *this = Sparse<T> (nr, nc - (ub - lb));
1259 } 1265 }
1260 else 1266 else
1261 { 1267 {
1262 const Sparse<T> tmp = *this; 1268 const Sparse<T> tmp = *this;
1263 octave_idx_type lbi = tmp.cidx(lb), ubi = tmp.cidx(ub), 1269 octave_idx_type lbi = tmp.cidx(lb);
1264 new_nz = nz - (ubi - lbi); 1270 octave_idx_type ubi = tmp.cidx(ub);
1271 octave_idx_type new_nz = nz - (ubi - lbi);
1265 1272
1266 *this = Sparse<T> (nr, nc - (ub - lb), new_nz); 1273 *this = Sparse<T> (nr, nc - (ub - lb), new_nz);
1267 copy_or_memcpy (lbi, tmp.data (), data ()); 1274 copy_or_memcpy (lbi, tmp.data (), data ());
1268 copy_or_memcpy (lbi, tmp.ridx (), ridx ()); 1275 copy_or_memcpy (lbi, tmp.ridx (), ridx ());
1269 copy_or_memcpy (nz - ubi, tmp.data () + ubi, xdata () + lbi); 1276 copy_or_memcpy (nz - ubi, tmp.data () + ubi, xdata () + lbi);
1530 { 1537 {
1531 // Special-case a contiguous range. 1538 // Special-case a contiguous range.
1532 octave_idx_type lbi = cidx (lb), ubi = cidx (ub), new_nz = ubi - lbi; 1539 octave_idx_type lbi = cidx (lb), ubi = cidx (ub), new_nz = ubi - lbi;
1533 retval = Sparse<T> (1, ub - lb, new_nz); 1540 retval = Sparse<T> (1, ub - lb, new_nz);
1534 copy_or_memcpy (new_nz, data () + lbi, retval.data ()); 1541 copy_or_memcpy (new_nz, data () + lbi, retval.data ());
1535 fill_or_memset (new_nz, static_cast<octave_idx_type> (0), retval.ridx ()); 1542 fill_or_memset (new_nz, static_cast<octave_idx_type> (0),
1543 retval.ridx ());
1536 mx_inline_sub (ub - lb + 1, retval.cidx (), cidx () + lb, lbi); 1544 mx_inline_sub (ub - lb + 1, retval.cidx (), cidx () + lb, lbi);
1537 } 1545 }
1538 else 1546 else
1539 { 1547 {
1540 // Sparse row vectors occupy O(nr) storage anyway, so let's just 1548 // Sparse row vectors occupy O(nr) storage anyway, so let's just
1547 if (nr != 0 && idx.is_scalar ()) 1555 if (nr != 0 && idx.is_scalar ())
1548 retval = Sparse<T> (1, 1, elem (idx(0) % nr, idx(0) / nr)); 1556 retval = Sparse<T> (1, 1, elem (idx(0) % nr, idx(0) / nr));
1549 else 1557 else
1550 { 1558 {
1551 // Indexing a non-vector sparse matrix by linear indexing. 1559 // Indexing a non-vector sparse matrix by linear indexing.
1552 // I suppose this is rare (and it may easily overflow), so let's take the easy way, 1560 // I suppose this is rare (and it may easily overflow), so let's take
1553 // and reshape first to column vector, which is already handled above. 1561 // the easy way, and reshape first to column vector, which is already
1562 // handled above.
1554 retval = index (idx_vector::colon).index (idx); 1563 retval = index (idx_vector::colon).index (idx);
1555 // In this case we're supposed to always inherit the shape, but column(row) doesn't do 1564 // In this case we're supposed to always inherit the shape, but
1556 // it, so we'll do it instead. 1565 // column(row) doesn't do it, so we'll do it instead.
1557 if (idx_dims(0) == 1 && idx_dims(1) != 1) 1566 if (idx_dims(0) == 1 && idx_dims(1) != 1)
1558 retval = retval.transpose (); 1567 retval = retval.transpose ();
1559 } 1568 }
1560 } 1569 }
1561 1570
1562 return retval; 1571 return retval;
1563 } 1572 }
1564 1573
1565 template <class T> 1574 template <class T>
1566 Sparse<T> 1575 Sparse<T>
1567 Sparse<T>::index (const idx_vector& idx_i, const idx_vector& idx_j, bool resize_ok) const 1576 Sparse<T>::index (const idx_vector& idx_i, const idx_vector& idx_j,
1577 bool resize_ok) const
1568 { 1578 {
1569 Sparse<T> retval; 1579 Sparse<T> retval;
1570 1580
1571 assert (ndims () == 2); 1581 assert (ndims () == 2);
1572 1582
1594 else 1604 else
1595 gripe_index_out_of_range (2, 2, idx_j.extent (nc), nc); 1605 gripe_index_out_of_range (2, 2, idx_j.extent (nc), nc);
1596 } 1606 }
1597 else if (nr == 1 && nc == 1) 1607 else if (nr == 1 && nc == 1)
1598 { 1608 {
1599 // Scalars stored as sparse matrices occupy more memory than 1609 // Scalars stored as sparse matrices occupy more memory than
1600 // a scalar, so let's just convert the matrix to full, index, 1610 // a scalar, so let's just convert the matrix to full, index,
1601 // and sparsify the result. 1611 // and sparsify the result.
1602 1612
1603 retval = Sparse<T> (array_value ().index (idx_i, idx_j)); 1613 retval = Sparse<T> (array_value ().index (idx_i, idx_j));
1604 } 1614 }
1605 else if (idx_i.is_colon ()) 1615 else if (idx_i.is_colon ())
1631 1641
1632 // Copy data & indices. 1642 // Copy data & indices.
1633 for (octave_idx_type j = 0; j < m; j++) 1643 for (octave_idx_type j = 0; j < m; j++)
1634 { 1644 {
1635 octave_idx_type ljj = cidx (idx_j(j)); 1645 octave_idx_type ljj = cidx (idx_j(j));
1636 octave_idx_type lj = retval.xcidx (j), nzj = retval.xcidx (j+1) - lj; 1646 octave_idx_type lj = retval.xcidx (j);
1647 octave_idx_type nzj = retval.xcidx (j+1) - lj;
1648
1637 copy_or_memcpy (nzj, data () + ljj, retval.data () + lj); 1649 copy_or_memcpy (nzj, data () + ljj, retval.data () + lj);
1638 copy_or_memcpy (nzj, ridx () + ljj, retval.ridx () + lj); 1650 copy_or_memcpy (nzj, ridx () + ljj, retval.ridx () + lj);
1639 } 1651 }
1640 } 1652 }
1641 } 1653 }
1654 retval = Sparse<T> (1, m); 1666 retval = Sparse<T> (1, m);
1655 OCTAVE_LOCAL_BUFFER (octave_idx_type, ij, m); 1667 OCTAVE_LOCAL_BUFFER (octave_idx_type, ij, m);
1656 for (octave_idx_type j = 0; j < m; j++) 1668 for (octave_idx_type j = 0; j < m; j++)
1657 { 1669 {
1658 octave_quit (); 1670 octave_quit ();
1659 octave_idx_type jj = idx_j(j), lj = cidx (jj), nzj = cidx (jj+1) - cidx (jj); 1671 octave_idx_type jj = idx_j(j);
1672 octave_idx_type lj = cidx (jj);
1673 octave_idx_type nzj = cidx (jj+1) - cidx (jj);
1674
1660 // Scalar index - just a binary lookup. 1675 // Scalar index - just a binary lookup.
1661 octave_idx_type i = lblookup (ridx () + lj, nzj, ii); 1676 octave_idx_type i = lblookup (ridx () + lj, nzj, ii);
1662 if (i < nzj && ridx (i+lj) == ii) 1677 if (i < nzj && ridx (i+lj) == ii)
1663 { 1678 {
1664 ij[j] = i + lj; 1679 ij[j] = i + lj;
1687 OCTAVE_LOCAL_BUFFER (octave_idx_type, li, m); 1702 OCTAVE_LOCAL_BUFFER (octave_idx_type, li, m);
1688 OCTAVE_LOCAL_BUFFER (octave_idx_type, ui, m); 1703 OCTAVE_LOCAL_BUFFER (octave_idx_type, ui, m);
1689 for (octave_idx_type j = 0; j < m; j++) 1704 for (octave_idx_type j = 0; j < m; j++)
1690 { 1705 {
1691 octave_quit (); 1706 octave_quit ();
1692 octave_idx_type jj = idx_j(j), lj = cidx (jj), nzj = cidx (jj+1) - cidx (jj); 1707 octave_idx_type jj = idx_j(j);
1708 octave_idx_type lj = cidx (jj);
1709 octave_idx_type nzj = cidx (jj+1) - cidx (jj);
1693 octave_idx_type lij, uij; 1710 octave_idx_type lij, uij;
1711
1694 // Lookup indices. 1712 // Lookup indices.
1695 li[j] = lij = lblookup (ridx () + lj, nzj, lb) + lj; 1713 li[j] = lij = lblookup (ridx () + lj, nzj, lb) + lj;
1696 ui[j] = uij = lblookup (ridx () + lj, nzj, ub) + lj; 1714 ui[j] = uij = lblookup (ridx () + lj, nzj, ub) + lj;
1697 retval.xcidx (j+1) = retval.xcidx (j) + ui[j] - li[j]; 1715 retval.xcidx (j+1) = retval.xcidx (j) + ui[j] - li[j];
1698 } 1716 }
1710 } 1728 }
1711 } 1729 }
1712 } 1730 }
1713 else if (idx_i.is_permutation (nr)) 1731 else if (idx_i.is_permutation (nr))
1714 { 1732 {
1715 // The columns preserve their length, we just need to renumber and sort them. 1733 // The columns preserve their length, just need to renumber and sort them.
1716 // Count new nonzero elements. 1734 // Count new nonzero elements.
1717 retval = Sparse<T> (nr, m); 1735 retval = Sparse<T> (nr, m);
1718 for (octave_idx_type j = 0; j < m; j++) 1736 for (octave_idx_type j = 0; j < m; j++)
1719 { 1737 {
1720 octave_idx_type jj = idx_j(j); 1738 octave_idx_type jj = idx_j(j);
1729 { 1747 {
1730 // It's nr:-1:1. Just flip all columns. 1748 // It's nr:-1:1. Just flip all columns.
1731 for (octave_idx_type j = 0; j < m; j++) 1749 for (octave_idx_type j = 0; j < m; j++)
1732 { 1750 {
1733 octave_quit (); 1751 octave_quit ();
1734 octave_idx_type jj = idx_j(j), lj = cidx (jj), nzj = cidx (jj+1) - cidx (jj); 1752 octave_idx_type jj = idx_j(j);
1753 octave_idx_type lj = cidx (jj);
1754 octave_idx_type nzj = cidx (jj+1) - cidx (jj);
1735 octave_idx_type li = retval.xcidx (j), uj = lj + nzj - 1; 1755 octave_idx_type li = retval.xcidx (j), uj = lj + nzj - 1;
1736 for (octave_idx_type i = 0; i < nzj; i++) 1756 for (octave_idx_type i = 0; i < nzj; i++)
1737 { 1757 {
1738 retval.xdata (li + i) = data (uj - i); // Copy in reverse order. 1758 retval.xdata (li + i) = data (uj - i); // Copy in reverse order.
1739 retval.xridx (li + i) = nr - 1 - ridx (uj - i); // Ditto with transform. 1759 retval.xridx (li + i) = nr - 1 - ridx (uj - i); // Ditto with transform.
1751 octave_idx_type *rri = retval.ridx (); 1771 octave_idx_type *rri = retval.ridx ();
1752 1772
1753 for (octave_idx_type j = 0; j < m; j++) 1773 for (octave_idx_type j = 0; j < m; j++)
1754 { 1774 {
1755 octave_quit (); 1775 octave_quit ();
1756 octave_idx_type jj = idx_j(j), lj = cidx (jj), nzj = cidx (jj+1) - cidx (jj); 1776 octave_idx_type jj = idx_j(j);
1777 octave_idx_type lj = cidx (jj);
1778 octave_idx_type nzj = cidx (jj+1) - cidx (jj);
1757 octave_idx_type li = retval.xcidx (j); 1779 octave_idx_type li = retval.xcidx (j);
1758 // Scatter the column, transform indices. 1780 // Scatter the column, transform indices.
1759 for (octave_idx_type i = 0; i < nzj; i++) 1781 for (octave_idx_type i = 0; i < nzj; i++)
1760 scb[rri[li + i] = iinv[ridx (lj + i)]] = data (lj + i); 1782 scb[rri[li + i] = iinv[ridx (lj + i)]] = data (lj + i);
1761 1783
1844 // Adding/overwriting elements, enough capacity allocated. 1866 // Adding/overwriting elements, enough capacity allocated.
1845 1867
1846 if (new_nz > nz) 1868 if (new_nz > nz)
1847 { 1869 {
1848 // Make room first. 1870 // Make room first.
1849 std::copy_backward (data () + ui, data () + nz, data () + nz + rnz); 1871 std::copy_backward (data () + ui, data () + nz,
1850 std::copy_backward (ridx () + ui, ridx () + nz, ridx () + nz + rnz); 1872 data () + nz + rnz);
1873 std::copy_backward (ridx () + ui, ridx () + nz,
1874 ridx () + nz + rnz);
1851 } 1875 }
1852 1876
1853 // Copy data and adjust indices from rhs. 1877 // Copy data and adjust indices from rhs.
1854 copy_or_memcpy (rnz, rhs.data (), data () + li); 1878 copy_or_memcpy (rnz, rhs.data (), data () + li);
1855 mx_inline_add (rnz, ridx () + li, rhs.ridx (), lb); 1879 mx_inline_add (rnz, ridx () + li, rhs.ridx (), lb);
1868 // new stuff ... 1892 // new stuff ...
1869 copy_or_memcpy (rnz, rhs.data (), data () + li); 1893 copy_or_memcpy (rnz, rhs.data (), data () + li);
1870 mx_inline_add (rnz, ridx () + li, rhs.ridx (), lb); 1894 mx_inline_add (rnz, ridx () + li, rhs.ridx (), lb);
1871 1895
1872 // ...tail 1896 // ...tail
1873 copy_or_memcpy (nz - ui, tmp.data () + ui, data () + li + rnz); 1897 copy_or_memcpy (nz - ui, tmp.data () + ui,
1874 copy_or_memcpy (nz - ui, tmp.ridx () + ui, ridx () + li + rnz); 1898 data () + li + rnz);
1899 copy_or_memcpy (nz - ui, tmp.ridx () + ui,
1900 ridx () + li + rnz);
1875 } 1901 }
1876 1902
1877 cidx (1) = new_nz; 1903 cidx (1) = new_nz;
1878 } 1904 }
1879 else if (idx.is_range () && idx.increment () == -1) 1905 else if (idx.is_range () && idx.increment () == -1)
1951 octave_idx_type nz = nnz (); 1977 octave_idx_type nz = nnz ();
1952 1978
1953 octave_idx_type n = rhs.rows (); 1979 octave_idx_type n = rhs.rows ();
1954 octave_idx_type m = rhs.columns (); 1980 octave_idx_type m = rhs.columns ();
1955 1981
1956 // FIXME -- this should probably be written more like the 1982 // FIXME: this should probably be written more like the
1957 // Array<T>::assign function... 1983 // Array<T>::assign function...
1958 1984
1959 bool orig_zero_by_zero = (nr == 0 && nc == 0); 1985 bool orig_zero_by_zero = (nr == 0 && nc == 0);
1960 1986
1961 if (orig_zero_by_zero || (idx_i.length (nr) == n && idx_j.length (nc) == m)) 1987 if (orig_zero_by_zero || (idx_i.length (nr) == n && idx_j.length (nc) == m))
2031 } 2057 }
2032 2058
2033 // Copy data and indices from rhs. 2059 // Copy data and indices from rhs.
2034 copy_or_memcpy (rnz, rhs.data (), data () + li); 2060 copy_or_memcpy (rnz, rhs.data (), data () + li);
2035 copy_or_memcpy (rnz, rhs.ridx (), ridx () + li); 2061 copy_or_memcpy (rnz, rhs.ridx (), ridx () + li);
2036 mx_inline_add (ub - lb, cidx () + lb + 1, rhs.cidx () + 1, li); 2062 mx_inline_add (ub - lb, cidx () + lb + 1, rhs.cidx () + 1,
2063 li);
2037 2064
2038 assert (nnz () == new_nz); 2065 assert (nnz () == new_nz);
2039 } 2066 }
2040 else 2067 else
2041 { 2068 {
2050 copy_or_memcpy (lb, tmp.cidx () + 1, cidx () + 1); 2077 copy_or_memcpy (lb, tmp.cidx () + 1, cidx () + 1);
2051 2078
2052 // new stuff... 2079 // new stuff...
2053 copy_or_memcpy (rnz, rhs.data (), data () + li); 2080 copy_or_memcpy (rnz, rhs.data (), data () + li);
2054 copy_or_memcpy (rnz, rhs.ridx (), ridx () + li); 2081 copy_or_memcpy (rnz, rhs.ridx (), ridx () + li);
2055 mx_inline_add (ub - lb, cidx () + lb + 1, rhs.cidx () + 1, li); 2082 mx_inline_add (ub - lb, cidx () + lb + 1, rhs.cidx () + 1,
2083 li);
2056 2084
2057 // ...tail. 2085 // ...tail.
2058 copy_or_memcpy (nz - ui, tmp.data () + ui, data () + li + rnz); 2086 copy_or_memcpy (nz - ui, tmp.data () + ui,
2059 copy_or_memcpy (nz - ui, tmp.ridx () + ui, ridx () + li + rnz); 2087 data () + li + rnz);
2060 mx_inline_add (nc - ub, cidx () + ub + 1, tmp.cidx () + ub + 1, new_nz - nz); 2088 copy_or_memcpy (nz - ui, tmp.ridx () + ui,
2089 ridx () + li + rnz);
2090 mx_inline_add (nc - ub, cidx () + ub + 1,
2091 tmp.cidx () + ub + 1, new_nz - nz);
2061 2092
2062 assert (nnz () == new_nz); 2093 assert (nnz () == new_nz);
2063 } 2094 }
2064 } 2095 }
2065 else if (idx_j.is_range () && idx_j.increment () == -1) 2096 else if (idx_j.is_range () && idx_j.increment () == -1)
2066 { 2097 {
2067 // It's s(:,u:-1:l) = r. Reverse the assignment. 2098 // It's s(:,u:-1:l) = r. Reverse the assignment.
2068 assign (idx_i, idx_j.sorted (), rhs.index (idx_i, idx_vector (m - 1, 0, -1))); 2099 assign (idx_i, idx_j.sorted (),
2100 rhs.index (idx_i, idx_vector (m - 1, 0, -1)));
2069 } 2101 }
2070 else if (idx_j.is_permutation (nc)) 2102 else if (idx_j.is_permutation (nc))
2071 { 2103 {
2072 *this = rhs.index (idx_i, idx_j.inverse_permutation (nc)); 2104 *this = rhs.index (idx_i, idx_j.inverse_permutation (nc));
2073 } 2105 }
2116 2148
2117 } 2149 }
2118 } 2150 }
2119 else if (nc == 1 && idx_j.is_colon_equiv (nc) && idx_i.is_vector ()) 2151 else if (nc == 1 && idx_j.is_colon_equiv (nc) && idx_i.is_vector ())
2120 { 2152 {
2121 // It's actually vector indexing. The 1D assign is specialized for that. 2153 // It's just vector indexing. The 1D assign is specialized for that.
2122 assign (idx_i, rhs); 2154 assign (idx_i, rhs);
2123 } 2155 }
2124 else if (idx_j.is_colon ()) 2156 else if (idx_j.is_colon ())
2125 { 2157 {
2126 if (idx_i.is_permutation (nr)) 2158 if (idx_i.is_permutation (nr))
2235 v += ns; 2267 v += ns;
2236 mridx += ns; 2268 mridx += ns;
2237 } 2269 }
2238 2270
2239 if (dim > 0) 2271 if (dim > 0)
2240 m = m.transpose (); 2272 m = m.transpose ();
2241 2273
2242 return m; 2274 return m;
2243 } 2275 }
2244 2276
2245 template <class T> 2277 template <class T>
2556 2588
2557 switch (dim) 2589 switch (dim)
2558 { 2590 {
2559 case 0: 2591 case 0:
2560 { 2592 {
2561 // sparse vertcat. This is not efficiently handled by assignment, so 2593 // sparse vertcat. This is not efficiently handled by assignment,
2562 // we'll do it directly. 2594 // so we'll do it directly.
2563 octave_idx_type l = 0; 2595 octave_idx_type l = 0;
2564 for (octave_idx_type j = 0; j < dv(1); j++) 2596 for (octave_idx_type j = 0; j < dv(1); j++)
2565 { 2597 {
2566 octave_quit (); 2598 octave_quit ();
2567 2599
2598 // Skipping empty matrices. See the comment in Array.cc. 2630 // Skipping empty matrices. See the comment in Array.cc.
2599 if (sparse_list[i].is_empty ()) 2631 if (sparse_list[i].is_empty ())
2600 continue; 2632 continue;
2601 2633
2602 octave_idx_type u = l + sparse_list[i].columns (); 2634 octave_idx_type u = l + sparse_list[i].columns ();
2603 retval.assign (idx_vector::colon, idx_vector (l, u), sparse_list[i]); 2635 retval.assign (idx_vector::colon, idx_vector (l, u),
2636 sparse_list[i]);
2604 l = u; 2637 l = u;
2605 } 2638 }
2606 2639
2607 break; 2640 break;
2608 } 2641 }