Mercurial > octave-nkf

comparison liboctave/array/Sparse.cc @ 18084:8e056300994b

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Follow coding convention of defining and initializing only 1 variable per line in liboctave. * liboctave/array/Array-b.cc, liboctave/array/Array-util.cc, liboctave/array/Array.cc, liboctave/array/CDiagMatrix.cc, liboctave/array/CMatrix.cc, liboctave/array/CSparse.cc, liboctave/array/MDiagArray2.cc, liboctave/array/MatrixType.cc, liboctave/array/PermMatrix.cc, liboctave/array/Sparse.cc, liboctave/array/Sparse.h, liboctave/array/boolSparse.cc, liboctave/array/dDiagMatrix.cc, liboctave/array/dMatrix.cc, liboctave/array/dSparse.cc, liboctave/array/dim-vector.cc, liboctave/array/fCDiagMatrix.cc, liboctave/array/fCMatrix.cc, liboctave/array/fDiagMatrix.cc, liboctave/array/fMatrix.cc, liboctave/array/idx-vector.cc, liboctave/array/idx-vector.h, liboctave/numeric/CmplxLU.cc, liboctave/numeric/CmplxQR.cc, liboctave/numeric/base-qr.cc, liboctave/numeric/bsxfun-defs.cc, liboctave/numeric/bsxfun.h, liboctave/numeric/dbleLU.cc, liboctave/numeric/dbleQR.cc, liboctave/numeric/fCmplxLU.cc, liboctave/numeric/fCmplxQR.cc, liboctave/numeric/floatLU.cc, liboctave/numeric/floatQR.cc, liboctave/numeric/lo-specfun.cc, liboctave/numeric/oct-convn.cc, liboctave/numeric/oct-norm.cc, liboctave/numeric/sparse-dmsolve.cc, liboctave/operators/mx-inlines.cc, liboctave/operators/mx-op-defs.h, liboctave/util/caseless-str.h, liboctave/util/kpse.cc, liboctave/util/lo-utils.cc, liboctave/util/oct-binmap.h, liboctave/util/oct-cmplx.h, liboctave/util/oct-inttypes.cc, liboctave/util/oct-inttypes.h, liboctave/util/oct-sort.cc: Follow coding convention of defining and initializing only 1 variable per line in liboctave.
author Rik <rik@octave.org>
date Wed, 04 Dec 2013 22:13:18 -0800
parents 938f01339043
children a86d608c413c 9254ff4036b2
comparison
equal deleted inserted replaced
18083:938f01339043 18084:8e056300994b
142 void 142 void
143 Sparse<T>::SparseRep::maybe_compress (bool remove_zeros) 143 Sparse<T>::SparseRep::maybe_compress (bool remove_zeros)
144 { 144 {
145 if (remove_zeros) 145 if (remove_zeros)
146 { 146 {
147 octave_idx_type i = 0, k = 0; 147 octave_idx_type i = 0;
148 octave_idx_type k = 0;
148 for (octave_idx_type j = 1; j <= ncols; j++) 149 for (octave_idx_type j = 1; j <= ncols; j++)
149 { 150 {
150 octave_idx_type u = c[j]; 151 octave_idx_type u = c[j];
151 for (i = i; i < u; i++) 152 for (i = i; i < u; i++)
152 if (d[i] != T ()) 153 if (d[i] != T ())
301 (*current_liboctave_error_handler) 302 (*current_liboctave_error_handler)
302 ("sparse: column index %d out of bound %d", r.extent (nc), nc); 303 ("sparse: column index %d out of bound %d", r.extent (nc), nc);
303 304
304 dimensions = dim_vector (nr, nc); 305 dimensions = dim_vector (nr, nc);
305 306
306 octave_idx_type n = a.numel (), rl = r.length (nr), cl = c.length (nc); 307 octave_idx_type n = a.numel ();
308 octave_idx_type rl = r.length (nr);
309 octave_idx_type cl = c.length (nc);
307 bool a_scalar = n == 1; 310 bool a_scalar = n == 1;
308 if (a_scalar) 311 if (a_scalar)
309 { 312 {
310 if (rl != 1) 313 if (rl != 1)
311 n = rl; 314 n = rl;
359 octave_idx_type *rri = ridx (); 362 octave_idx_type *rri = ridx ();
360 T *rrd = data (); 363 T *rrd = data ();
361 364
362 octave_quit (); 365 octave_quit ();
363 366
364 octave_idx_type k = -1, l = -1; 367 octave_idx_type k = -1;
368 octave_idx_type l = -1;
365 369
366 if (sum_terms) 370 if (sum_terms)
367 { 371 {
368 // Sum repeated indices. 372 // Sum repeated indices.
369 for (octave_idx_type i = 0; i < n; i++) 373 for (octave_idx_type i = 0; i < n; i++)
393 } 397 }
394 398
395 } 399 }
396 else 400 else
397 { 401 {
398 idx_vector rr = r, cc = c; 402 idx_vector rr = r;
399 const octave_idx_type *rd = rr.raw (), *cd = cc.raw (); 403 idx_vector cc = c;
404 const octave_idx_type *rd = rr.raw ();
405 const octave_idx_type *cd = cc.raw ();
400 OCTAVE_LOCAL_BUFFER_INIT (octave_idx_type, ci, nc+1, 0); 406 OCTAVE_LOCAL_BUFFER_INIT (octave_idx_type, ci, nc+1, 0);
401 ci[0] = 0; 407 ci[0] = 0;
402 // Bin counts of column indices. 408 // Bin counts of column indices.
403 for (octave_idx_type i = 0; i < n; i++) 409 for (octave_idx_type i = 0; i < n; i++)
404 ci[cd[i]+1]++; 410 ci[cd[i]+1]++;
423 // Subsorts. We don't need a stable sort, all values are equal. 429 // Subsorts. We don't need a stable sort, all values are equal.
424 xcidx (0) = 0; 430 xcidx (0) = 0;
425 for (octave_idx_type j = 0; j < nc; j++) 431 for (octave_idx_type j = 0; j < nc; j++)
426 { 432 {
427 std::sort (sidx + ci[j], sidx + ci[j+1]); 433 std::sort (sidx + ci[j], sidx + ci[j+1]);
428 octave_idx_type l = -1, nzj = 0; 434 octave_idx_type l = -1;
435 octave_idx_type nzj = 0;
429 // Count. 436 // Count.
430 for (octave_idx_type i = ci[j]; i < ci[j+1]; i++) 437 for (octave_idx_type i = ci[j]; i < ci[j+1]; i++)
431 { 438 {
432 octave_idx_type k = sidx[i]; 439 octave_idx_type k = sidx[i];
433 if (k != l) 440 if (k != l)
488 Array<octave_idx_type> rsi; 495 Array<octave_idx_type> rsi;
489 idx_vector rs = r.sorted (rsi); 496 idx_vector rs = r.sorted (rsi);
490 497
491 octave_quit (); 498 octave_quit ();
492 499
493 const octave_idx_type *rd = rs.raw (), *rdi = rsi.data (); 500 const octave_idx_type *rd = rs.raw ();
501 const octave_idx_type *rdi = rsi.data ();
494 // Count unique indices. 502 // Count unique indices.
495 octave_idx_type new_nz = 1; 503 octave_idx_type new_nz = 1;
496 for (octave_idx_type i = 1; i < n; i++) 504 for (octave_idx_type i = 1; i < n; i++)
497 new_nz += rd[i-1] != rd[i]; 505 new_nz += rd[i-1] != rd[i];
498 // Allocate result. 506 // Allocate result.
535 543
536 maybe_compress (true); 544 maybe_compress (true);
537 } 545 }
538 else 546 else
539 { 547 {
540 idx_vector rr = r, cc = c; 548 idx_vector rr = r;
541 const octave_idx_type *rd = rr.raw (), *cd = cc.raw (); 549 idx_vector cc = c;
550 const octave_idx_type *rd = rr.raw ();
551 const octave_idx_type *cd = cc.raw ();
542 OCTAVE_LOCAL_BUFFER_INIT (octave_idx_type, ci, nc+1, 0); 552 OCTAVE_LOCAL_BUFFER_INIT (octave_idx_type, ci, nc+1, 0);
543 ci[0] = 0; 553 ci[0] = 0;
544 // Bin counts of column indices. 554 // Bin counts of column indices.
545 for (octave_idx_type i = 0; i < n; i++) 555 for (octave_idx_type i = 0; i < n; i++)
546 ci[cd[i]+1]++; 556 ci[cd[i]+1]++;
570 // Subsorts. We don't need a stable sort, the second index stabilizes it. 580 // Subsorts. We don't need a stable sort, the second index stabilizes it.
571 xcidx (0) = 0; 581 xcidx (0) = 0;
572 for (octave_idx_type j = 0; j < nc; j++) 582 for (octave_idx_type j = 0; j < nc; j++)
573 { 583 {
574 std::sort (spairs + ci[j], spairs + ci[j+1]); 584 std::sort (spairs + ci[j], spairs + ci[j+1]);
575 octave_idx_type l = -1, nzj = 0; 585 octave_idx_type l = -1;
586 octave_idx_type nzj = 0;
576 // Count. 587 // Count.
577 for (octave_idx_type i = ci[j]; i < ci[j+1]; i++) 588 for (octave_idx_type i = ci[j]; i < ci[j+1]; i++)
578 { 589 {
579 octave_idx_type k = spairs[i].first; 590 octave_idx_type k = spairs[i].first;
580 if (k != l) 591 if (k != l)
899 910
900 template <class T> 911 template <class T>
901 void 912 void
902 Sparse<T>::resize1 (octave_idx_type n) 913 Sparse<T>::resize1 (octave_idx_type n)
903 { 914 {
904 octave_idx_type nr = rows (), nc = cols (); 915 octave_idx_type nr = rows ();
916 octave_idx_type nc = cols ();
905 917
906 if (nr == 0) 918 if (nr == 0)
907 resize (1, std::max (nc, n)); 919 resize (1, std::max (nc, n));
908 else if (nc == 0) 920 else if (nc == 0)
909 resize (nr, (n + nr - 1) / nr); // Ain't it wicked? 921 resize (nr, (n + nr - 1) / nr); // Ain't it wicked?
948 // Sparse rep. It is not good for anything in there. 960 // Sparse rep. It is not good for anything in there.
949 make_unique (); 961 make_unique ();
950 962
951 if (r < rows ()) 963 if (r < rows ())
952 { 964 {
953 octave_idx_type i = 0, k = 0; 965 octave_idx_type i = 0;
966 octave_idx_type k = 0;
954 for (octave_idx_type j = 1; j <= rep->ncols; j++) 967 for (octave_idx_type j = 1; j <= rep->ncols; j++)
955 { 968 {
956 octave_idx_type u = xcidx (j); 969 octave_idx_type u = xcidx (j);
957 for (i = i; i < u; i++) 970 for (i = i; i < u; i++)
958 if (xridx (i) < r) 971 if (xridx (i) < r)
1181 { 1194 {
1182 OCTAVE_LOCAL_BUFFER (octave_idx_type, ridx_new, nz); 1195 OCTAVE_LOCAL_BUFFER (octave_idx_type, ridx_new, nz);
1183 OCTAVE_LOCAL_BUFFER (T, data_new, nz); 1196 OCTAVE_LOCAL_BUFFER (T, data_new, nz);
1184 idx_vector sidx = idx.sorted (true); 1197 idx_vector sidx = idx.sorted (true);
1185 const octave_idx_type *sj = sidx.raw (); 1198 const octave_idx_type *sj = sidx.raw ();
1186 octave_idx_type sl = sidx.length (nel), nz_new = 0, j = 0; 1199 octave_idx_type sl = sidx.length (nel);
1200 octave_idx_type nz_new = 0;
1201 octave_idx_type j = 0;
1187 for (octave_idx_type i = 0; i < nz; i++) 1202 for (octave_idx_type i = 0; i < nz; i++)
1188 { 1203 {
1189 octave_idx_type r = tmp.ridx (i); 1204 octave_idx_type r = tmp.ridx (i);
1190 for (; j < sl && sj[j] < r; j++) ; 1205 for (; j < sl && sj[j] < r; j++) ;
1191 if (j == sl || sj[j] > r) 1206 if (j == sl || sj[j] > r)
1486 // Vector transpose is cheap, so do it right here. 1501 // Vector transpose is cheap, so do it right here.
1487 const Array<octave_idx_type> idxa = (idx_dims(0) == 1 1502 const Array<octave_idx_type> idxa = (idx_dims(0) == 1
1488 ? idx.as_array ().transpose () 1503 ? idx.as_array ().transpose ()
1489 : idx.as_array ()); 1504 : idx.as_array ());
1490 1505
1491 octave_idx_type new_nr = idxa.rows (), new_nc = idxa.cols (); 1506 octave_idx_type new_nr = idxa.rows ();
1507 octave_idx_type new_nc = idxa.cols ();
1492 1508
1493 // Lookup. 1509 // Lookup.
1494 // FIXME: Could specialize for sorted idx? 1510 // FIXME: Could specialize for sorted idx?
1495 NoAlias< Array<octave_idx_type> > lidx (dim_vector (new_nr, new_nc)); 1511 NoAlias< Array<octave_idx_type> > lidx (dim_vector (new_nr, new_nc));
1496 for (octave_idx_type i = 0; i < new_nr*new_nc; i++) 1512 for (octave_idx_type i = 0; i < new_nr*new_nc; i++)
1534 if (idx.is_scalar ()) 1550 if (idx.is_scalar ())
1535 retval = Sparse<T> (1, 1, elem (0, idx(0))); 1551 retval = Sparse<T> (1, 1, elem (0, idx(0)));
1536 else if (idx.is_cont_range (nel, lb, ub)) 1552 else if (idx.is_cont_range (nel, lb, ub))
1537 { 1553 {
1538 // Special-case a contiguous range. 1554 // Special-case a contiguous range.
1539 octave_idx_type lbi = cidx (lb), ubi = cidx (ub), new_nz = ubi - lbi; 1555 octave_idx_type lbi = cidx (lb);
1556 octave_idx_type ubi = cidx (ub);
1557 octave_idx_type new_nz = ubi - lbi;
1540 retval = Sparse<T> (1, ub - lb, new_nz); 1558 retval = Sparse<T> (1, ub - lb, new_nz);
1541 copy_or_memcpy (new_nz, data () + lbi, retval.data ()); 1559 copy_or_memcpy (new_nz, data () + lbi, retval.data ());
1542 fill_or_memset (new_nz, static_cast<octave_idx_type> (0), 1560 fill_or_memset (new_nz, static_cast<octave_idx_type> (0),
1543 retval.ridx ()); 1561 retval.ridx ());
1544 mx_inline_sub (ub - lb + 1, retval.cidx (), cidx () + lb, lbi); 1562 mx_inline_sub (ub - lb + 1, retval.cidx (), cidx () + lb, lbi);
1619 if (idx_j.is_colon ()) 1637 if (idx_j.is_colon ())
1620 retval = *this; // Shallow copy. 1638 retval = *this; // Shallow copy.
1621 else if (idx_j.is_cont_range (nc, lb, ub)) 1639 else if (idx_j.is_cont_range (nc, lb, ub))
1622 { 1640 {
1623 // Special-case a contiguous range. 1641 // Special-case a contiguous range.
1624 octave_idx_type lbi = cidx (lb), ubi = cidx (ub), new_nz = ubi - lbi; 1642 octave_idx_type lbi = cidx (lb);
1643 octave_idx_type ubi = cidx (ub);
1644 octave_idx_type new_nz = ubi - lbi;
1625 retval = Sparse<T> (nr, ub - lb, new_nz); 1645 retval = Sparse<T> (nr, ub - lb, new_nz);
1626 copy_or_memcpy (new_nz, data () + lbi, retval.data ()); 1646 copy_or_memcpy (new_nz, data () + lbi, retval.data ());
1627 copy_or_memcpy (new_nz, ridx () + lbi, retval.ridx ()); 1647 copy_or_memcpy (new_nz, ridx () + lbi, retval.ridx ());
1628 mx_inline_sub (ub - lb + 1, retval.cidx (), cidx () + lb, lbi); 1648 mx_inline_sub (ub - lb + 1, retval.cidx (), cidx () + lb, lbi);
1629 } 1649 }
1750 { 1770 {
1751 octave_quit (); 1771 octave_quit ();
1752 octave_idx_type jj = idx_j(j); 1772 octave_idx_type jj = idx_j(j);
1753 octave_idx_type lj = cidx (jj); 1773 octave_idx_type lj = cidx (jj);
1754 octave_idx_type nzj = cidx (jj+1) - cidx (jj); 1774 octave_idx_type nzj = cidx (jj+1) - cidx (jj);
1755 octave_idx_type li = retval.xcidx (j), uj = lj + nzj - 1; 1775 octave_idx_type li = retval.xcidx (j);
1776 octave_idx_type uj = lj + nzj - 1;
1756 for (octave_idx_type i = 0; i < nzj; i++) 1777 for (octave_idx_type i = 0; i < nzj; i++)
1757 { 1778 {
1758 retval.xdata (li + i) = data (uj - i); // Copy in reverse order. 1779 retval.xdata (li + i) = data (uj - i); // Copy in reverse order.
1759 retval.xridx (li + i) = nr - 1 - ridx (uj - i); // Ditto with transform. 1780 retval.xridx (li + i) = nr - 1 - ridx (uj - i); // Ditto with transform.
1760 } 1781 }
1857 { 1878 {
1858 // Special-case a contiguous range. 1879 // Special-case a contiguous range.
1859 // Look-up indices first. 1880 // Look-up indices first.
1860 octave_idx_type li = lblookup (ridx (), nz, lb); 1881 octave_idx_type li = lblookup (ridx (), nz, lb);
1861 octave_idx_type ui = lblookup (ridx (), nz, ub); 1882 octave_idx_type ui = lblookup (ridx (), nz, ub);
1862 octave_idx_type rnz = rhs.nnz (), new_nz = nz - (ui - li) + rnz; 1883 octave_idx_type rnz = rhs.nnz ();
1884 octave_idx_type new_nz = nz - (ui - li) + rnz;
1863 1885
1864 if (new_nz >= nz && new_nz <= capacity ()) 1886 if (new_nz >= nz && new_nz <= capacity ())
1865 { 1887 {
1866 // Adding/overwriting elements, enough capacity allocated. 1888 // Adding/overwriting elements, enough capacity allocated.
1867 1889
2037 if (idx_j.is_colon ()) 2059 if (idx_j.is_colon ())
2038 *this = rhs; // Shallow copy. 2060 *this = rhs; // Shallow copy.
2039 else if (idx_j.is_cont_range (nc, lb, ub)) 2061 else if (idx_j.is_cont_range (nc, lb, ub))
2040 { 2062 {
2041 // Special-case a contiguous range. 2063 // Special-case a contiguous range.
2042 octave_idx_type li = cidx (lb), ui = cidx (ub); 2064 octave_idx_type li = cidx (lb);
2043 octave_idx_type rnz = rhs.nnz (), new_nz = nz - (ui - li) + rnz; 2065 octave_idx_type ui = cidx (ub);
2066 octave_idx_type rnz = rhs.nnz ();
2067 octave_idx_type new_nz = nz - (ui - li) + rnz;
2044 2068
2045 if (new_nz >= nz && new_nz <= capacity ()) 2069 if (new_nz >= nz && new_nz <= capacity ())
2046 { 2070 {
2047 // Adding/overwriting elements, enough capacity allocated. 2071 // Adding/overwriting elements, enough capacity allocated.
2048 2072
2127 change_capacity (nnz ()); 2151 change_capacity (nnz ());
2128 2152
2129 // Merge columns. 2153 // Merge columns.
2130 for (octave_idx_type i = 0; i < nc; i++) 2154 for (octave_idx_type i = 0; i < nc; i++)
2131 { 2155 {
2132 octave_idx_type l = xcidx (i), u = xcidx (i+1), j = jsav[i]; 2156 octave_idx_type l = xcidx (i);
2157 octave_idx_type u = xcidx (i+1);
2158 octave_idx_type j = jsav[i];
2133 if (j >= 0) 2159 if (j >= 0)
2134 { 2160 {
2135 // from rhs 2161 // from rhs
2136 octave_idx_type k = rhs.cidx (j); 2162 octave_idx_type k = rhs.cidx (j);
2137 copy_or_memcpy (u - l, rhs.data () + k, xdata () + l); 2163 copy_or_memcpy (u - l, rhs.data () + k, xdata () + l);
2603 const Sparse<T>& spi = sparse_list[i]; 2629 const Sparse<T>& spi = sparse_list[i];
2604 // Skipping empty matrices. See the comment in Array.cc. 2630 // Skipping empty matrices. See the comment in Array.cc.
2605 if (spi.is_empty ()) 2631 if (spi.is_empty ())
2606 continue; 2632 continue;
2607 2633
2608 octave_idx_type kl = spi.cidx (j), ku = spi.cidx (j+1); 2634 octave_idx_type kl = spi.cidx (j);
2635 octave_idx_type ku = spi.cidx (j+1);
2609 for (octave_idx_type k = kl; k < ku; k++, l++) 2636 for (octave_idx_type k = kl; k < ku; k++, l++)
2610 { 2637 {
2611 retval.xridx (l) = spi.ridx (k) + rcum; 2638 retval.xridx (l) = spi.ridx (k) + rcum;
2612 retval.xdata (l) = spi.data (k); 2639 retval.xdata (l) = spi.data (k);
2613 } 2640 }