//////////////////////////////////////////////////////////////////////// // // Copyright (C) 1998-2021 The Octave Project Developers // // See the file COPYRIGHT.md in the top-level directory of this // distribution or . // // This file is part of Octave. // // Octave is free software: you can redistribute it and/or modify it // under the terms of the GNU General Public License as published by // the Free Software Foundation, either version 3 of the License, or // (at your option) any later version. // // Octave is distributed in the hope that it will be useful, but // WITHOUT ANY WARRANTY; without even the implied warranty of // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the // GNU General Public License for more details. // // You should have received a copy of the GNU General Public License // along with Octave; see the file COPYING. If not, see // . // //////////////////////////////////////////////////////////////////////// #if ! defined (octave_Sparse_h) #define octave_Sparse_h 1 #include "octave-config.h" #include #include #include #include #include class PermMatrix; template class Array; template void delete_array (void* p); using deleter_t = void (*)(void *); // Two dimensional sparse class. Handles the reference counting for // all the derived classes. // forward declare template with visibility attribute template class OCTAVE_API Sparse; template class Sparse { public: typedef T element_type; protected: //-------------------------------------------------------------------- // The real representation of all Sparse arrays. //-------------------------------------------------------------------- class SparseRep { public: T *m_data; octave_idx_type *m_ridx; octave_idx_type *m_cidx; octave_idx_type m_nzmax; octave_idx_type m_nrows; octave_idx_type m_ncols; octave::refcount m_count; deleter_t m_data_deleter; deleter_t m_ridx_deleter; deleter_t m_cidx_deleter; SparseRep (void) : m_data (new T [1]), m_ridx (new octave_idx_type [1] {}), m_cidx (new octave_idx_type [1] {}), m_nzmax (1), m_nrows (0), m_ncols (0), m_count (1), m_data_deleter (delete_array), m_ridx_deleter (delete_array), m_cidx_deleter (delete_array) { } SparseRep (octave_idx_type n) : m_data (new T [1]), m_ridx (new octave_idx_type [1] {}), m_cidx (new octave_idx_type [n+1] {}), m_nzmax (1), m_nrows (n), m_ncols (n), m_count (1), m_data_deleter (delete_array), m_ridx_deleter (delete_array), m_cidx_deleter (delete_array) { } SparseRep (octave_idx_type nr, octave_idx_type nc, octave_idx_type nz = 1) : m_data (nz > 0 ? new T [nz] : new T [1]), m_ridx (nz > 0 ? new octave_idx_type [nz] {} : new octave_idx_type [1] {}), m_cidx (new octave_idx_type [nc+1] {}), m_nzmax (nz > 0 ? nz : 1), m_nrows (nr), m_ncols (nc), m_count (1), m_data_deleter (delete_array), m_ridx_deleter (delete_array), m_cidx_deleter (delete_array) { } SparseRep (octave_idx_type nr, octave_idx_type nc, octave_idx_type nz, const T *d, const octave_idx_type *r, const octave_idx_type *c) : m_data (new T [nz]), m_ridx (new octave_idx_type [nz] {}), m_cidx (new octave_idx_type [nc+1] {}), m_nzmax (nz), m_nrows (nr), m_ncols (nc), m_count (1), m_data_deleter (delete_array), m_ridx_deleter (delete_array), m_cidx_deleter (delete_array) { std::copy_n (d, nz, m_data); std::copy_n (r, nz, m_ridx); std::copy_n (c, m_ncols + 1, m_cidx); } template SparseRep (octave_idx_type nr, octave_idx_type nc, octave_idx_type nz, const U *d, const octave_idx_type *r, const octave_idx_type *c) : m_data (new T [nz]), m_ridx (new octave_idx_type [nz] {}), m_cidx (new octave_idx_type [nc+1] {}), m_nzmax (nz), m_nrows (nr), m_ncols (nc), m_count (1), m_data_deleter (delete_array), m_ridx_deleter (delete_array), m_cidx_deleter (delete_array) { std::copy_n (d, nz, m_data); std::copy_n (r, nz, m_ridx); std::copy_n (c, nc + 1, m_cidx); } template SparseRep (const dim_vector& dv, octave_idx_type nz, U *d, octave_idx_type *r, octave_idx_type *c, deleter_t d_del, deleter_t r_del, deleter_t c_del) : m_data (d), m_ridx (r), m_cidx (c), m_nzmax (nz), m_nrows (dv(0)), m_ncols (dv(1)), m_count (1), m_data_deleter (d_del), m_ridx_deleter (r_del), m_cidx_deleter (c_del) { } SparseRep (const SparseRep& a) : m_data (new T [a.m_nzmax]), m_ridx (new octave_idx_type [a.m_nzmax]), m_cidx (new octave_idx_type [a.m_ncols + 1]), m_nzmax (a.m_nzmax), m_nrows (a.m_nrows), m_ncols (a.m_ncols), m_count (1), m_data_deleter (delete_array), m_ridx_deleter (delete_array), m_cidx_deleter (delete_array) { octave_idx_type nz = a.nnz (); std::copy_n (a.m_data, nz, m_data); std::copy_n (a.m_ridx, nz, m_ridx); std::copy_n (a.m_cidx, m_ncols + 1, m_cidx); } ~SparseRep (void) { m_data_deleter (m_data); m_ridx_deleter (m_ridx); m_cidx_deleter (m_cidx); } octave_idx_type nzmax (void) const { return m_nzmax; } octave_idx_type nnz (void) const { return m_cidx[m_ncols]; } octave_idx_type rows (void) const { return m_nrows; } octave_idx_type cols (void) const { return m_ncols; } octave_idx_type columns (void) const { return m_ncols; } OCTAVE_API T& elem (octave_idx_type r, octave_idx_type c); OCTAVE_API T celem (octave_idx_type r, octave_idx_type c) const; T& data (octave_idx_type i) { return m_data[i]; } T cdata (octave_idx_type i) const { return m_data[i]; } octave_idx_type& ridx (octave_idx_type i) { return m_ridx[i]; } octave_idx_type cridx (octave_idx_type i) const { return m_ridx[i]; } octave_idx_type& cidx (octave_idx_type i) { return m_cidx[i]; } octave_idx_type ccidx (octave_idx_type i) const { return m_cidx[i]; } OCTAVE_API void maybe_compress (bool remove_zeros); OCTAVE_API void change_length (octave_idx_type nz); OCTAVE_API bool indices_ok (void) const; OCTAVE_API bool any_element_is_nan (void) const; // Prefer nzmax. octave_idx_type length (void) const { return m_nzmax; } template friend class Sparse; // No assignment! SparseRep& operator = (const SparseRep&) = delete; }; //-------------------------------------------------------------------- void make_unique (void) { if (m_rep->m_count > 1) { SparseRep *r = new SparseRep (*m_rep); if (--m_rep->m_count == 0) delete m_rep; m_rep = r; } } protected: typename Sparse::SparseRep *m_rep; dim_vector m_dimensions; private: static OCTAVE_API typename Sparse::SparseRep * nil_rep (void); public: Sparse (void) : m_rep (nil_rep ()), m_dimensions (dim_vector (0, 0)) { m_rep->m_count++; } explicit Sparse (octave_idx_type n) : m_rep (new typename Sparse::SparseRep (n)), m_dimensions (dim_vector (n, n)) { } explicit Sparse (octave_idx_type nr, octave_idx_type nc) : m_rep (new typename Sparse::SparseRep (nr, nc)), m_dimensions (dim_vector (nr, nc)) { } explicit OCTAVE_API Sparse (octave_idx_type nr, octave_idx_type nc, T val); Sparse (const dim_vector& dv, octave_idx_type nz) : m_rep (new typename Sparse::SparseRep (dv(0), dv(1), nz)), m_dimensions (dv) { } Sparse (octave_idx_type nr, octave_idx_type nc, octave_idx_type nz) : m_rep (new typename Sparse::SparseRep (nr, nc, nz)), m_dimensions (dim_vector (nr, nc)) { } // Construct a Sparse array from pointers to externally allocated // arrays of values and indices. PTR, RIDX, and CIDX must be // allocated with operator new. The Sparse object takes ownership of // these arrays and will delete them when the Sparse object is // deleted. The dimension vector DV must be consistent with the sizes // of the allocated PTR, CIDX, and RIDX arrays. Sparse (const dim_vector& dv, octave_idx_type nz, T *ptr, octave_idx_type *ridx, octave_idx_type *cidx, deleter_t d_del, deleter_t r_del, deleter_t c_del) : m_rep (new typename Sparse::SparseRep (dv, nz, ptr, ridx, cidx, d_del, r_del, c_del)), m_dimensions (dv) { } // Both SparseMatrix and SparseBoolMatrix need this ctor, and this // is their only common ancestor. explicit OCTAVE_API Sparse (const PermMatrix& a); // Type conversion case. Preserves nzmax. template Sparse (const Sparse& a) : m_rep (new typename Sparse::SparseRep (a.rows (), a.cols (), a.nzmax (), a.data (), a.ridx (), a.cidx ())), m_dimensions (a.dims ()) { } // No type conversion case. Sparse (const Sparse& a) : m_rep (a.m_rep), m_dimensions (a.m_dimensions) { m_rep->m_count++; } public: OCTAVE_API Sparse (const dim_vector& dv); OCTAVE_API Sparse (const Sparse& a, const dim_vector& dv); OCTAVE_API Sparse (const Array& a, const octave::idx_vector& r, const octave::idx_vector& c, octave_idx_type nr = -1, octave_idx_type nc = -1, bool sum_terms = true, octave_idx_type nzm = -1); // Sparsify a normal matrix OCTAVE_API Sparse (const Array& a); virtual ~Sparse (void); OCTAVE_API Sparse& operator = (const Sparse& a); //! Amount of storage for nonzero elements. //! This may differ from the actual number of elements, see nnz(). octave_idx_type nzmax (void) const { return m_rep->nzmax (); } //! Actual number of nonzero terms. octave_idx_type nnz (void) const { return m_rep->nnz (); } // Querying the number of elements (incl. zeros) may overflow the index type, // so don't do it unless you really need it. octave_idx_type numel (void) const { return m_dimensions.safe_numel (); } octave_idx_type dim1 (void) const { return m_dimensions(0); } octave_idx_type dim2 (void) const { return m_dimensions(1); } octave_idx_type rows (void) const { return dim1 (); } octave_idx_type cols (void) const { return dim2 (); } octave_idx_type columns (void) const { return dim2 (); } octave_idx_type get_row_index (octave_idx_type k) { return ridx (k); } octave_idx_type get_col_index (octave_idx_type k) { octave_idx_type ret = 0; while (cidx (ret+1) < k) ret++; return ret; } std::size_t byte_size (void) const { return (static_cast (cols () + 1) * sizeof (octave_idx_type) + static_cast (nzmax ()) * (sizeof (T) + sizeof (octave_idx_type))); } dim_vector dims (void) const { return m_dimensions; } Sparse squeeze (void) const { return *this; } OCTAVE_API octave_idx_type compute_index (const Array& ra_idx) const; OCTAVE_NORETURN OCTAVE_API T range_error (const char *fcn, octave_idx_type n) const; OCTAVE_NORETURN OCTAVE_API T& range_error (const char *fcn, octave_idx_type n); OCTAVE_NORETURN OCTAVE_API T range_error (const char *fcn, octave_idx_type i, octave_idx_type j) const; OCTAVE_NORETURN OCTAVE_API T& range_error (const char *fcn, octave_idx_type i, octave_idx_type j); OCTAVE_NORETURN OCTAVE_API T range_error (const char *fcn, const Array& ra_idx) const; OCTAVE_NORETURN OCTAVE_API T& range_error (const char *fcn, const Array& ra_idx); // No checking, even for multiple references, ever. T& xelem (octave_idx_type n) { octave_idx_type i = n % rows (); octave_idx_type j = n / rows (); return xelem (i, j); } T xelem (octave_idx_type n) const { octave_idx_type i = n % rows (); octave_idx_type j = n / rows (); return xelem (i, j); } T& xelem (octave_idx_type i, octave_idx_type j) { return m_rep->elem (i, j); } T xelem (octave_idx_type i, octave_idx_type j) const { return m_rep->celem (i, j); } T& xelem (const Array& ra_idx) { return xelem (compute_index (ra_idx)); } T xelem (const Array& ra_idx) const { return xelem (compute_index (ra_idx)); } // FIXME: would be nice to fix this so that we don't unnecessarily force a // copy, but that is not so easy, and I see no clean way to do it. T& checkelem (octave_idx_type n) { if (n < 0 || n >= numel ()) range_error ("T& Sparse::checkelem", n); else { make_unique (); return xelem (n); } } T& checkelem (octave_idx_type i, octave_idx_type j) { if (i < 0 || j < 0 || i >= dim1 () || j >= dim2 ()) range_error ("T& Sparse::checkelem", i, j); else { make_unique (); return xelem (i, j); } } T& checkelem (const Array& ra_idx) { octave_idx_type i = compute_index (ra_idx); if (i < 0) range_error ("T& Sparse::checkelem", ra_idx); else return elem (i); } T& elem (octave_idx_type n) { make_unique (); return xelem (n); } T& elem (octave_idx_type i, octave_idx_type j) { make_unique (); return xelem (i, j); } T& elem (const Array& ra_idx) { return Sparse::elem (compute_index (ra_idx)); } T& operator () (octave_idx_type n) { return elem (n); } T& operator () (octave_idx_type i, octave_idx_type j) { return elem (i, j); } T& operator () (const Array& ra_idx) { return elem (ra_idx); } T checkelem (octave_idx_type n) const { if (n < 0 || n >= numel ()) range_error ("T Sparse::checkelem", n); else return xelem (n); } T checkelem (octave_idx_type i, octave_idx_type j) const { if (i < 0 || j < 0 || i >= dim1 () || j >= dim2 ()) range_error ("T Sparse::checkelem", i, j); else return xelem (i, j); } T checkelem (const Array& ra_idx) const { octave_idx_type i = compute_index (ra_idx); if (i < 0) range_error ("T Sparse::checkelem", ra_idx); else return Sparse::elem (i); } T elem (octave_idx_type n) const { return xelem (n); } T elem (octave_idx_type i, octave_idx_type j) const { return xelem (i, j); } T elem (const Array& ra_idx) const { return Sparse::elem (compute_index (ra_idx)); } T operator () (octave_idx_type n) const { return elem (n); } T operator () (octave_idx_type i, octave_idx_type j) const { return elem (i, j); } T operator () (const Array& ra_idx) const { return elem (ra_idx); } Sparse maybe_compress (bool remove_zeros = false) { if (remove_zeros) make_unique (); // Need to unshare because elements are removed. m_rep->maybe_compress (remove_zeros); return (*this); } OCTAVE_API Sparse reshape (const dim_vector& new_dims) const; OCTAVE_API Sparse permute (const Array& vec, bool inv = false) const; Sparse ipermute (const Array& vec) const { return permute (vec, true); } OCTAVE_API void resize1 (octave_idx_type n); OCTAVE_API void resize (octave_idx_type r, octave_idx_type c); OCTAVE_API void resize (const dim_vector& dv); void change_capacity (octave_idx_type nz) { if (nz < nnz ()) make_unique (); // Unshare now because elements will be truncated. m_rep->change_length (nz); } OCTAVE_API Sparse& insert (const Sparse& a, octave_idx_type r, octave_idx_type c); OCTAVE_API Sparse& insert (const Sparse& a, const Array& idx); bool issquare (void) const { return (dim1 () == dim2 ()); } bool isempty (void) const { return (rows () < 1 || cols () < 1); } OCTAVE_API Sparse transpose (void) const; T * data (void) { make_unique (); return m_rep->m_data; } T& data (octave_idx_type i) { make_unique (); return m_rep->data (i); } T * xdata (void) { return m_rep->m_data; } T& xdata (octave_idx_type i) { return m_rep->data (i); } T data (octave_idx_type i) const { return m_rep->data (i); } // FIXME: shouldn't this be returning const T*? T * data (void) const { return m_rep->m_data; } octave_idx_type * ridx (void) { make_unique (); return m_rep->m_ridx; } octave_idx_type& ridx (octave_idx_type i) { make_unique (); return m_rep->ridx (i); } octave_idx_type * xridx (void) { return m_rep->m_ridx; } octave_idx_type& xridx (octave_idx_type i) { return m_rep->ridx (i); } octave_idx_type ridx (octave_idx_type i) const { return m_rep->cridx (i); } // FIXME: shouldn't this be returning const octave_idx_type*? octave_idx_type * ridx (void) const { return m_rep->m_ridx; } octave_idx_type * cidx (void) { make_unique (); return m_rep->m_cidx; } octave_idx_type& cidx (octave_idx_type i) { make_unique (); return m_rep->cidx (i); } octave_idx_type * xcidx (void) { return m_rep->m_cidx; } octave_idx_type& xcidx (octave_idx_type i) { return m_rep->cidx (i); } octave_idx_type cidx (octave_idx_type i) const { return m_rep->ccidx (i); } // FIXME: shouldn't this be returning const octave_idx_type*? octave_idx_type * cidx (void) const { return m_rep->m_cidx; } octave_idx_type ndims (void) const { return m_dimensions.ndims (); } OCTAVE_API void delete_elements (const octave::idx_vector& i); OCTAVE_API void delete_elements (int dim, const octave::idx_vector& i); OCTAVE_API void delete_elements (const octave::idx_vector& i, const octave::idx_vector& j); OCTAVE_API Sparse index (const octave::idx_vector& i, bool resize_ok = false) const; OCTAVE_API Sparse index (const octave::idx_vector& i, const octave::idx_vector& j, bool resize_ok = false) const; OCTAVE_API void assign (const octave::idx_vector& i, const Sparse& rhs); OCTAVE_API void assign (const octave::idx_vector& i, const T& rhs); OCTAVE_API void assign (const octave::idx_vector& i, const octave::idx_vector& j, const Sparse& rhs); OCTAVE_API void assign (const octave::idx_vector& i, const octave::idx_vector& j, const T& rhs); OCTAVE_API void print_info (std::ostream& os, const std::string& prefix) const; OCTAVE_API Sparse sort (octave_idx_type dim = 0, sortmode mode = ASCENDING) const; OCTAVE_API Sparse sort (Array& sidx, octave_idx_type dim = 0, sortmode mode = ASCENDING) const; OCTAVE_API Sparse diag (octave_idx_type k = 0) const; // dim = -1 and dim = -2 are special; see Array::cat description. static OCTAVE_API Sparse cat (int dim, octave_idx_type n, const Sparse *sparse_list); OCTAVE_API Array array_value (void) const; // Generic any/all test functionality with arbitrary predicate. template bool test (F fcn) const { return octave::any_all_test (fcn, data (), nnz ()); } // Simpler calls. template bool test_any (F fcn) const { return test (fcn); } template bool test_all (F fcn) const { return test (fcn); } // Overloads for function references. bool test_any (bool (&fcn) (T)) const { return test (fcn); } bool test_any (bool (&fcn) (const T&)) const { return test (fcn); } bool test_all (bool (&fcn) (T)) const { return test (fcn); } bool test_all (bool (&fcn) (const T&)) const { return test (fcn); } template Sparse map (F fcn) const { Sparse result; U f_zero = fcn (0.0); if (f_zero != 0.0) { octave_idx_type nr = rows (); octave_idx_type nc = cols (); result = Sparse (nr, nc, f_zero); for (octave_idx_type j = 0; j < nc; j++) for (octave_idx_type i = cidx (j); i < cidx (j+1); i++) { octave_quit (); /* Use data instead of elem for better performance. */ result.data (ridx (i) + j * nr) = fcn (data (i)); } result.maybe_compress (true); } else { octave_idx_type nz = nnz (); octave_idx_type nr = rows (); octave_idx_type nc = cols (); result = Sparse (nr, nc, nz); octave_idx_type ii = 0; result.cidx (ii) = 0; for (octave_idx_type j = 0; j < nc; j++) { for (octave_idx_type i = cidx (j); i < cidx (j+1); i++) { U val = fcn (data (i)); if (val != 0.0) { result.data (ii) = val; result.ridx (ii++) = ridx (i); } octave_quit (); } result.cidx (j+1) = ii; } result.maybe_compress (false); } return result; } // Overloads for function references. template Sparse map (U (&fcn) (T)) const { return map (fcn); } template Sparse map (U (&fcn) (const T&)) const { return map (fcn); } bool indices_ok (void) const { return m_rep->indices_ok (); } bool any_element_is_nan (void) const { return m_rep->any_element_is_nan (); } }; template OCTAVE_API std::istream& read_sparse_matrix (std::istream& is, Sparse& a, T (*read_fcn) (std::istream&)); #endif