////////////////////////////////////////////////////////////////////////
//
// 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