Mercurial > octave-nkf

comparison liboctave/Sparse.cc @ 10421:99e9bae2d81e

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
improve sparse indexing interface
author Jaroslav Hajek <highegg@gmail.com>
date Thu, 18 Mar 2010 14:55:52 +0100
parents 1766c133674c
children 0677c5d80b77
comparison
equal deleted inserted replaced
10420:3373fdc0b14a 10421:99e9bae2d81e
1609 return retval; 1609 return retval;
1610 } 1610 }
1611 1611
1612 template <class T> 1612 template <class T>
1613 Sparse<T> 1613 Sparse<T>
1614 Sparse<T>::index (idx_vector& idx_arg, int resize_ok) const 1614 Sparse<T>::index (const idx_vector& idx_arg, bool resize_ok) const
1615 { 1615 {
1616 Sparse<T> retval; 1616 Sparse<T> retval;
1617 1617
1618 assert (ndims () == 2); 1618 assert (ndims () == 2);
1619 1619
1620 octave_idx_type nr = dim1 (); 1620 octave_idx_type nr = dim1 ();
1621 octave_idx_type nc = dim2 (); 1621 octave_idx_type nc = dim2 ();
1622 octave_idx_type nz = nnz (); 1622 octave_idx_type nz = nnz ();
1623 1623
1624 // Use safe_numel so that we get an error if the matrix is too big to be indexed.
1624 octave_idx_type orig_len = nr * nc; 1625 octave_idx_type orig_len = nr * nc;
1625 1626
1626 dim_vector idx_orig_dims = idx_arg.orig_dimensions (); 1627 dim_vector idx_orig_dims = idx_arg.orig_dimensions ();
1627 1628
1628 octave_idx_type idx_orig_rows = idx_arg.orig_rows (); 1629 octave_idx_type idx_orig_rows = idx_arg.orig_rows ();
1629 octave_idx_type idx_orig_columns = idx_arg.orig_columns (); 1630 octave_idx_type idx_orig_columns = idx_arg.orig_columns ();
1630 1631
1631 if (idx_orig_dims.length () > 2) 1632 if (idx_orig_dims.length () > 2)
1632 (*current_liboctave_error_handler) 1633 (*current_liboctave_error_handler)
1633 ("Sparse<T>::index: Can not index Sparse<T> with an N-D Array"); 1634 ("cannot index sparse matrix with an N-D Array");
1634 else if (idx_arg.is_colon ()) 1635 else if (idx_arg.is_colon ())
1635 { 1636 {
1636 // Fast magic colon processing. 1637 // Fast magic colon processing.
1637 retval = Sparse<T> (nr * nc, 1, nz); 1638 retval = Sparse<T> (nr * nc, 1, nz);
1638 1639
1644 retval.xridx(j) = ridx(j) + i * nr; 1645 retval.xridx(j) = ridx(j) + i * nr;
1645 } 1646 }
1646 retval.xcidx(0) = 0; 1647 retval.xcidx(0) = 0;
1647 retval.xcidx(1) = nz; 1648 retval.xcidx(1) = nz;
1648 } 1649 }
1650 else if (! resize_ok && idx_arg.extent (length ()) > length ())
1651 {
1652 gripe_index_out_of_range (1, 1, idx_arg.extent (orig_len), orig_len);
1653 }
1649 else if (nr == 1 && nc == 1) 1654 else if (nr == 1 && nc == 1)
1650 { 1655 {
1651 // You have to be pretty sick to get to this bit of code, 1656 // You have to be pretty sick to get to this bit of code,
1652 // since you have a scalar stored as a sparse matrix, and 1657 // since you have a scalar stored as a sparse matrix, and
1653 // then want to make a dense matrix with sparse 1658 // then want to make a dense matrix with sparse
1654 // representation. Ok, we'll do it, but you deserve what 1659 // representation. Ok, we'll do it, but you deserve what
1655 // you get!! 1660 // you get!!
1656 octave_idx_type n = idx_arg.freeze (length (), "sparse vector", resize_ok); 1661 octave_idx_type n = idx_arg.length (length ());
1657 if (n == 0) 1662 if (n == 0)
1658 1663
1659 retval = Sparse<T> (idx_orig_dims); 1664 retval = Sparse<T> (idx_orig_dims);
1660 else if (nz < 1) 1665 else if (nz < 1)
1661 if (n >= idx_orig_dims.numel ()) 1666 if (n >= idx_orig_dims.numel ())
1707 { 1712 {
1708 // If indexing a vector with a matrix, return value has same 1713 // If indexing a vector with a matrix, return value has same
1709 // shape as the index. Otherwise, it has same orientation as 1714 // shape as the index. Otherwise, it has same orientation as
1710 // indexed object. 1715 // indexed object.
1711 octave_idx_type len = length (); 1716 octave_idx_type len = length ();
1712 octave_idx_type n = idx_arg.freeze (len, "sparse vector", resize_ok); 1717 octave_idx_type n = idx_arg.length (len);
1713 octave_idx_type l, u; 1718 octave_idx_type l, u;
1714 if (n == 0) 1719 if (n == 0)
1715 if (nr == 1) 1720 if (nr == 1)
1716 retval = Sparse<T> (dim_vector (1, 0)); 1721 retval = Sparse<T> (dim_vector (1, 0));
1717 else 1722 else
1915 ("Octave:fortran-indexing", "single index used for sparse matrix"); 1920 ("Octave:fortran-indexing", "single index used for sparse matrix");
1916 1921
1917 // This code is only for indexing matrices. The vector 1922 // This code is only for indexing matrices. The vector
1918 // cases are handled above. 1923 // cases are handled above.
1919 1924
1920 idx_arg.freeze (nr * nc, "matrix", resize_ok); 1925 octave_idx_type result_nr = idx_orig_rows;
1921 1926 octave_idx_type result_nc = idx_orig_columns;
1922 if (idx_arg) 1927
1923 { 1928 if (nz < 1)
1924 octave_idx_type result_nr = idx_orig_rows; 1929 retval = Sparse<T> (result_nr, result_nc);
1925 octave_idx_type result_nc = idx_orig_columns; 1930 else
1926 1931 {
1927 if (nz < 1) 1932 // Count number of non-zero elements
1928 retval = Sparse<T> (result_nr, result_nc); 1933 octave_idx_type new_nzmx = 0;
1929 else 1934 octave_idx_type kk = 0;
1930 { 1935 for (octave_idx_type j = 0; j < result_nc; j++)
1931 // Count number of non-zero elements 1936 {
1932 octave_idx_type new_nzmx = 0; 1937 for (octave_idx_type i = 0; i < result_nr; i++)
1933 octave_idx_type kk = 0; 1938 {
1934 for (octave_idx_type j = 0; j < result_nc; j++) 1939 octave_quit ();
1935 { 1940
1936 for (octave_idx_type i = 0; i < result_nr; i++) 1941 octave_idx_type ii = idx_arg.elem (kk++);
1937 { 1942 if (ii < orig_len)
1938 octave_quit (); 1943 {
1939 1944 octave_idx_type fr = ii % nr;
1940 octave_idx_type ii = idx_arg.elem (kk++); 1945 octave_idx_type fc = (ii - fr) / nr;
1941 if (ii < orig_len) 1946 for (octave_idx_type k = cidx(fc); k < cidx(fc+1); k++)
1942 { 1947 {
1943 octave_idx_type fr = ii % nr; 1948 if (ridx(k) == fr)
1944 octave_idx_type fc = (ii - fr) / nr; 1949 new_nzmx++;
1945 for (octave_idx_type k = cidx(fc); k < cidx(fc+1); k++) 1950 if (ridx(k) >= fr)
1951 break;
1952 }
1953 }
1954 }
1955 }
1956
1957 retval = Sparse<T> (result_nr, result_nc, new_nzmx);
1958
1959 kk = 0;
1960 octave_idx_type jj = 0;
1961 retval.xcidx(0) = 0;
1962 for (octave_idx_type j = 0; j < result_nc; j++)
1963 {
1964 for (octave_idx_type i = 0; i < result_nr; i++)
1965 {
1966 octave_quit ();
1967
1968 octave_idx_type ii = idx_arg.elem (kk++);
1969 if (ii < orig_len)
1970 {
1971 octave_idx_type fr = ii % nr;
1972 octave_idx_type fc = (ii - fr) / nr;
1973 for (octave_idx_type k = cidx(fc); k < cidx(fc+1); k++)
1974 {
1975 if (ridx(k) == fr)
1946 { 1976 {
1947 if (ridx(k) == fr) 1977 retval.xdata(jj) = data(k);
1948 new_nzmx++; 1978 retval.xridx(jj++) = i;
1949 if (ridx(k) >= fr)
1950 break;
1951 } 1979 }
1980 if (ridx(k) >= fr)
1981 break;
1952 } 1982 }
1953 } 1983 }
1954 } 1984 }
1955 1985 retval.xcidx(j+1) = jj;
1956 retval = Sparse<T> (result_nr, result_nc, new_nzmx); 1986 }
1957
1958 kk = 0;
1959 octave_idx_type jj = 0;
1960 retval.xcidx(0) = 0;
1961 for (octave_idx_type j = 0; j < result_nc; j++)
1962 {
1963 for (octave_idx_type i = 0; i < result_nr; i++)
1964 {
1965 octave_quit ();
1966
1967 octave_idx_type ii = idx_arg.elem (kk++);
1968 if (ii < orig_len)
1969 {
1970 octave_idx_type fr = ii % nr;
1971 octave_idx_type fc = (ii - fr) / nr;
1972 for (octave_idx_type k = cidx(fc); k < cidx(fc+1); k++)
1973 {
1974 if (ridx(k) == fr)
1975 {
1976 retval.xdata(jj) = data(k);
1977 retval.xridx(jj++) = i;
1978 }
1979 if (ridx(k) >= fr)
1980 break;
1981 }
1982 }
1983 }
1984 retval.xcidx(j+1) = jj;
1985 }
1986 }
1987 // idx_vector::freeze() printed an error message for us.
1988 } 1987 }
1989 } 1988 }
1990 1989
1991 return retval; 1990 return retval;
1992 } 1991 }
1998 struct idx_node *next; 1997 struct idx_node *next;
1999 }; 1998 };
2000 1999
2001 template <class T> 2000 template <class T>
2002 Sparse<T> 2001 Sparse<T>
2003 Sparse<T>::index (idx_vector& idx_i, idx_vector& idx_j, int resize_ok) const 2002 Sparse<T>::index (const idx_vector& idx_i, const idx_vector& idx_j, bool resize_ok) const
2004 { 2003 {
2005 Sparse<T> retval; 2004 Sparse<T> retval;
2006 2005
2007 assert (ndims () == 2); 2006 assert (ndims () == 2);
2008 2007
2009 octave_idx_type nr = dim1 (); 2008 octave_idx_type nr = dim1 ();
2010 octave_idx_type nc = dim2 (); 2009 octave_idx_type nc = dim2 ();
2011 2010
2012 octave_idx_type n = idx_i.freeze (nr, "row", resize_ok); 2011 octave_idx_type n = idx_i.length (nr);
2013 octave_idx_type m = idx_j.freeze (nc, "column", resize_ok); 2012 octave_idx_type m = idx_j.length (nc);
2014 2013
2015 if (idx_i && idx_j) 2014 if (! resize_ok && idx_i.extent (nr) > nr)
2016 { 2015 gripe_index_out_of_range (2, 1, idx_i.extent (nr), nr);
2017 if (idx_i.orig_empty () || idx_j.orig_empty () || n == 0 || m == 0) 2016 else if (! resize_ok && idx_j.extent (nc) > nc)
2017 gripe_index_out_of_range (2, 2, idx_j.extent (nc), nc);
2018 else
2019 {
2020 if (n == 0 || m == 0)
2018 { 2021 {
2019 retval.resize_no_fill (n, m); 2022 retval.resize_no_fill (n, m);
2020 } 2023 }
2021 else 2024 else
2022 { 2025 {
2023 int idx_i_colon = idx_i.is_colon_equiv (nr); 2026 bool idx_i_colon = idx_i.is_colon_equiv (nr);
2024 int idx_j_colon = idx_j.is_colon_equiv (nc); 2027 bool idx_j_colon = idx_j.is_colon_equiv (nc);
2025 2028
2026 if (idx_i_colon && idx_j_colon) 2029 if (idx_i_colon && idx_j_colon)
2027 { 2030 {
2028 retval = *this; 2031 retval = *this;
2029 } 2032 }
2211 } 2214 }
2212 } 2215 }
2213 } 2216 }
2214 } 2217 }
2215 } 2218 }
2216 // idx_vector::freeze() printed an error message for us.
2217 2219
2218 return retval; 2220 return retval;
2219 }
2220
2221 template <class T>
2222 Sparse<T>
2223 Sparse<T>::index (Array<idx_vector>& ra_idx, int resize_ok) const
2224 {
2225
2226 if (ra_idx.length () != 2)
2227 {
2228 (*current_liboctave_error_handler) ("range error for index");
2229 return *this;
2230 }
2231
2232 return index (ra_idx (0), ra_idx (1), resize_ok);
2233 } 2221 }
2234 2222
2235 // Can't use versions of these in Array.cc due to duplication of the 2223 // Can't use versions of these in Array.cc due to duplication of the
2236 // instantiations for Array<double and Sparse<double>, etc 2224 // instantiations for Array<double and Sparse<double>, etc
2237 template <class T> 2225 template <class T>