////////////////////////////////////////////////////////////////////////
//
// 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
// .
//
////////////////////////////////////////////////////////////////////////
// This file should not include config.h. It is only included in other
// C++ source files that should have included config.h before including
// this file.
#include
#include
#include
#include
#include
#include
#include
#include "Array.h"
#include "MArray.h"
#include "Array-util.h"
#include "Range.h"
#include "idx-vector.h"
#include "lo-error.h"
#include "quit.h"
#include "oct-locbuf.h"
#include "Sparse.h"
#include "sparse-util.h"
#include "oct-spparms.h"
#include "mx-inlines.cc"
#include "PermMatrix.h"
template
OCTAVE_API typename Sparse::SparseRep *
Sparse::nil_rep (void)
{
static typename Sparse::SparseRep nr;
return &nr;
}
template
OCTAVE_API
T&
Sparse::SparseRep::elem (octave_idx_type r, octave_idx_type c)
{
octave_idx_type i;
if (m_nzmax <= 0)
(*current_liboctave_error_handler)
("Sparse::SparseRep::elem (octave_idx_type, octave_idx_type): sparse matrix filled");
for (i = m_cidx[c]; i < m_cidx[c + 1]; i++)
if (m_ridx[i] == r)
return m_data[i];
else if (m_ridx[i] > r)
break;
// Ok, If we've gotten here, we're in trouble. Have to create a
// new element in the sparse array. This' gonna be slow!!!
if (m_cidx[m_ncols] == m_nzmax)
(*current_liboctave_error_handler)
("Sparse::SparseRep::elem (octave_idx_type, octave_idx_type): sparse matrix filled");
octave_idx_type to_move = m_cidx[m_ncols] - i;
if (to_move != 0)
{
for (octave_idx_type j = m_cidx[m_ncols]; j > i; j--)
{
m_data[j] = m_data[j-1];
m_ridx[j] = m_ridx[j-1];
}
}
for (octave_idx_type j = c + 1; j < m_ncols + 1; j++)
m_cidx[j] = m_cidx[j] + 1;
m_data[i] = 0.;
m_ridx[i] = r;
return m_data[i];
}
template
OCTAVE_API
T
Sparse::SparseRep::celem (octave_idx_type r, octave_idx_type c) const
{
if (m_nzmax > 0)
for (octave_idx_type i = m_cidx[c]; i < m_cidx[c + 1]; i++)
if (m_ridx[i] == r)
return m_data[i];
return T ();
}
template
OCTAVE_API
void
Sparse::SparseRep::maybe_compress (bool remove_zeros)
{
if (remove_zeros)
{
octave_idx_type i = 0;
octave_idx_type k = 0;
for (octave_idx_type j = 1; j <= m_ncols; j++)
{
octave_idx_type u = m_cidx[j];
for (; i < u; i++)
if (m_data[i] != T ())
{
m_data[k] = m_data[i];
m_ridx[k++] = m_ridx[i];
}
m_cidx[j] = k;
}
}
change_length (m_cidx[m_ncols]);
}
template
OCTAVE_API
void
Sparse::SparseRep::change_length (octave_idx_type nz)
{
for (octave_idx_type j = m_ncols; j > 0 && m_cidx[j] > nz; j--)
m_cidx[j] = nz;
// Always preserve space for 1 element.
nz = (nz > 0 ? nz : 1);
// Skip reallocation if we have less than 1/frac extra elements to discard.
static const int frac = 5;
if (nz > m_nzmax || nz < m_nzmax - m_nzmax/frac)
{
// Reallocate.
octave_idx_type min_nzmax = std::min (nz, m_nzmax);
octave_idx_type *new_ridx = new octave_idx_type [nz];
std::copy_n (m_ridx, min_nzmax, new_ridx);
m_ridx_deleter (m_ridx);
m_ridx = new_ridx;
m_ridx_deleter = delete_array;
T *new_data = new T [nz];
std::copy_n (m_data, min_nzmax, new_data);
m_data_deleter (m_data);
m_data = new_data;
m_data_deleter = delete_array;
m_nzmax = nz;
}
}
template
OCTAVE_API
bool
Sparse::SparseRep::indices_ok (void) const
{
return sparse_indices_ok (m_ridx, m_cidx, m_nrows, m_ncols, nnz ());
}
template
OCTAVE_API
bool
Sparse::SparseRep::any_element_is_nan (void) const
{
octave_idx_type nz = nnz ();
for (octave_idx_type i = 0; i < nz; i++)
if (octave::math::isnan (m_data[i]))
return true;
return false;
}
template
OCTAVE_API
Sparse::Sparse (octave_idx_type nr, octave_idx_type nc, T val)
: m_rep (nullptr), m_dimensions (dim_vector (nr, nc))
{
if (val != T ())
{
m_rep = new typename Sparse::SparseRep (nr, nc, m_dimensions.safe_numel ());
octave_idx_type ii = 0;
xcidx (0) = 0;
for (octave_idx_type j = 0; j < nc; j++)
{
for (octave_idx_type i = 0; i < nr; i++)
{
xdata (ii) = val;
xridx (ii++) = i;
}
xcidx (j+1) = ii;
}
}
else
{
m_rep = new typename Sparse::SparseRep (nr, nc, 0);
for (octave_idx_type j = 0; j < nc+1; j++)
xcidx (j) = 0;
}
}
template
OCTAVE_API
Sparse::Sparse (const PermMatrix& a)
: m_rep (new typename Sparse::SparseRep (a.rows (), a.cols (), a.rows ())),
m_dimensions (dim_vector (a.rows (), a.cols ()))
{
octave_idx_type n = a.rows ();
for (octave_idx_type i = 0; i <= n; i++)
cidx (i) = i;
const Array pv = a.col_perm_vec ();
for (octave_idx_type i = 0; i < n; i++)
ridx (i) = pv(i);
for (octave_idx_type i = 0; i < n; i++)
data (i) = 1.0;
}
template
OCTAVE_API
Sparse::Sparse (const dim_vector& dv)
: m_rep (nullptr), m_dimensions (dv)
{
if (dv.ndims () != 2)
(*current_liboctave_error_handler)
("Sparse::Sparse (const dim_vector&): dimension mismatch");
m_rep = new typename Sparse::SparseRep (dv(0), dv(1), 0);
}
template
OCTAVE_API
Sparse::Sparse (const Sparse& a, const dim_vector& dv)
: m_rep (nullptr), m_dimensions (dv)
{
// Work in unsigned long long to avoid overflow issues with numel
unsigned long long a_nel = static_cast(a.rows ()) *
static_cast(a.cols ());
unsigned long long dv_nel = static_cast(dv(0)) *
static_cast(dv(1));
if (a_nel != dv_nel)
(*current_liboctave_error_handler)
("Sparse::Sparse (const Sparse&, const dim_vector&): dimension mismatch");
dim_vector old_dims = a.dims ();
octave_idx_type new_nzmax = a.nnz ();
octave_idx_type new_nr = dv(0);
octave_idx_type new_nc = dv(1);
octave_idx_type old_nr = old_dims(0);
octave_idx_type old_nc = old_dims(1);
m_rep = new typename Sparse::SparseRep (new_nr, new_nc, new_nzmax);
octave_idx_type kk = 0;
xcidx (0) = 0;
for (octave_idx_type i = 0; i < old_nc; i++)
for (octave_idx_type j = a.cidx (i); j < a.cidx (i+1); j++)
{
octave_idx_type tmp = i * old_nr + a.ridx (j);
octave_idx_type ii = tmp % new_nr;
octave_idx_type jj = (tmp - ii) / new_nr;
for (octave_idx_type k = kk; k < jj; k++)
xcidx (k+1) = j;
kk = jj;
xdata (j) = a.data (j);
xridx (j) = ii;
}
for (octave_idx_type k = kk; k < new_nc; k++)
xcidx (k+1) = new_nzmax;
}
template
OCTAVE_API
Sparse::Sparse (const Array& a, const octave::idx_vector& r,
const octave::idx_vector& c, octave_idx_type nr,
octave_idx_type nc, bool sum_terms,
octave_idx_type nzm)
: m_rep (nullptr), m_dimensions ()
{
if (nr < 0)
nr = r.extent (0);
else if (r.extent (nr) > nr)
(*current_liboctave_error_handler)
("sparse: row index %" OCTAVE_IDX_TYPE_FORMAT "out of bound "
"%" OCTAVE_IDX_TYPE_FORMAT, r.extent (nr), nr);
if (nc < 0)
nc = c.extent (0);
else if (c.extent (nc) > nc)
(*current_liboctave_error_handler)
("sparse: column index %" OCTAVE_IDX_TYPE_FORMAT " out of bound "
"%" OCTAVE_IDX_TYPE_FORMAT, r.extent (nc), nc);
m_dimensions = dim_vector (nr, nc);
octave_idx_type n = a.numel ();
octave_idx_type rl = r.length (nr);
octave_idx_type cl = c.length (nc);
bool a_scalar = n == 1;
if (a_scalar)
{
if (rl != 1)
n = rl;
else if (cl != 1)
n = cl;
}
if ((rl != 1 && rl != n) || (cl != 1 && cl != n))
(*current_liboctave_error_handler) ("sparse: dimension mismatch");
// Only create m_rep after input validation to avoid memory leak.
m_rep = new typename Sparse::SparseRep (nr, nc, (nzm > 0 ? nzm : 0));
if (rl <= 1 && cl <= 1)
{
if (n == 1 && a(0) != T ())
{
change_capacity (nzm > 1 ? nzm : 1);
xridx (0) = r(0);
xdata (0) = a(0);
std::fill_n (xcidx () + c(0) + 1, nc - c(0), 1);
}
}
else if (a_scalar)
{
// This is completely specialized, because the sorts can be simplified.
T a0 = a(0);
if (a0 == T ())
{
// Do nothing, it's an empty matrix.
}
else if (cl == 1)
{
// Sparse column vector. Sort row indices.
octave::idx_vector rs = r.sorted ();
octave_quit ();
const octave_idx_type *rd = rs.raw ();
// Count unique indices.
octave_idx_type new_nz = 1;
for (octave_idx_type i = 1; i < n; i++)
new_nz += rd[i-1] != rd[i];
// Allocate result.
change_capacity (nzm > new_nz ? nzm : new_nz);
std::fill_n (xcidx () + c(0) + 1, nc - c(0), new_nz);
octave_idx_type *rri = ridx ();
T *rrd = data ();
octave_quit ();
octave_idx_type k = -1;
octave_idx_type l = -1;
if (sum_terms)
{
// Sum repeated indices.
for (octave_idx_type i = 0; i < n; i++)
{
if (rd[i] != l)
{
l = rd[i];
rri[++k] = rd[i];
rrd[k] = a0;
}
else
rrd[k] += a0;
}
}
else
{
// Pick the last one.
for (octave_idx_type i = 0; i < n; i++)
{
if (rd[i] != l)
{
l = rd[i];
rri[++k] = rd[i];
rrd[k] = a0;
}
}
}
}
else
{
octave::idx_vector rr = r;
octave::idx_vector cc = c;
const octave_idx_type *rd = rr.raw ();
const octave_idx_type *cd = cc.raw ();
OCTAVE_LOCAL_BUFFER_INIT (octave_idx_type, ci, nc+1, 0);
ci[0] = 0;
// Bin counts of column indices.
for (octave_idx_type i = 0; i < n; i++)
ci[cd[i]+1]++;
// Make them cumulative, shifted one to right.
for (octave_idx_type i = 1, s = 0; i <= nc; i++)
{
octave_idx_type s1 = s + ci[i];
ci[i] = s;
s = s1;
}
octave_quit ();
// Bucket sort.
OCTAVE_LOCAL_BUFFER (octave_idx_type, sidx, n);
for (octave_idx_type i = 0; i < n; i++)
if (rl == 1)
sidx[ci[cd[i]+1]++] = rd[0];
else
sidx[ci[cd[i]+1]++] = rd[i];
// Subsorts. We don't need a stable sort, all values are equal.
xcidx (0) = 0;
for (octave_idx_type j = 0; j < nc; j++)
{
std::sort (sidx + ci[j], sidx + ci[j+1]);
octave_idx_type l = -1;
octave_idx_type nzj = 0;
// Count.
for (octave_idx_type i = ci[j]; i < ci[j+1]; i++)
{
octave_idx_type k = sidx[i];
if (k != l)
{
l = k;
nzj++;
}
}
// Set column pointer.
xcidx (j+1) = xcidx (j) + nzj;
}
change_capacity (nzm > xcidx (nc) ? nzm : xcidx (nc));
octave_idx_type *rri = ridx ();
T *rrd = data ();
// Fill-in data.
for (octave_idx_type j = 0, jj = -1; j < nc; j++)
{
octave_quit ();
octave_idx_type l = -1;
if (sum_terms)
{
// Sum adjacent terms.
for (octave_idx_type i = ci[j]; i < ci[j+1]; i++)
{
octave_idx_type k = sidx[i];
if (k != l)
{
l = k;
rrd[++jj] = a0;
rri[jj] = k;
}
else
rrd[jj] += a0;
}
}
else
{
// Use the last one.
for (octave_idx_type i = ci[j]; i < ci[j+1]; i++)
{
octave_idx_type k = sidx[i];
if (k != l)
{
l = k;
rrd[++jj] = a0;
rri[jj] = k;
}
}
}
}
}
}
else if (cl == 1)
{
// Sparse column vector. Sort row indices.
Array rsi;
octave::idx_vector rs = r.sorted (rsi);
octave_quit ();
const octave_idx_type *rd = rs.raw ();
const octave_idx_type *rdi = rsi.data ();
// Count unique indices.
octave_idx_type new_nz = 1;
for (octave_idx_type i = 1; i < n; i++)
new_nz += rd[i-1] != rd[i];
// Allocate result.
change_capacity (nzm > new_nz ? nzm : new_nz);
std::fill_n (xcidx () + c(0) + 1, nc - c(0), new_nz);
octave_idx_type *rri = ridx ();
T *rrd = data ();
octave_quit ();
octave_idx_type k = 0;
rri[k] = rd[0];
rrd[k] = a(rdi[0]);
if (sum_terms)
{
// Sum repeated indices.
for (octave_idx_type i = 1; i < n; i++)
{
if (rd[i] != rd[i-1])
{
rri[++k] = rd[i];
rrd[k] = a(rdi[i]);
}
else
rrd[k] += a(rdi[i]);
}
}
else
{
// Pick the last one.
for (octave_idx_type i = 1; i < n; i++)
{
if (rd[i] != rd[i-1])
rri[++k] = rd[i];
rrd[k] = a(rdi[i]);
}
}
maybe_compress (true);
}
else
{
octave::idx_vector rr = r;
octave::idx_vector cc = c;
const octave_idx_type *rd = rr.raw ();
const octave_idx_type *cd = cc.raw ();
OCTAVE_LOCAL_BUFFER_INIT (octave_idx_type, ci, nc+1, 0);
ci[0] = 0;
// Bin counts of column indices.
for (octave_idx_type i = 0; i < n; i++)
ci[cd[i]+1]++;
// Make them cumulative, shifted one to right.
for (octave_idx_type i = 1, s = 0; i <= nc; i++)
{
octave_idx_type s1 = s + ci[i];
ci[i] = s;
s = s1;
}
octave_quit ();
typedef std::pair idx_pair;
// Bucket sort.
OCTAVE_LOCAL_BUFFER (idx_pair, spairs, n);
for (octave_idx_type i = 0; i < n; i++)
{
idx_pair& p = spairs[ci[cd[i]+1]++];
if (rl == 1)
p.first = rd[0];
else
p.first = rd[i];
p.second = i;
}
// Subsorts. We don't need a stable sort, the second index stabilizes it.
xcidx (0) = 0;
for (octave_idx_type j = 0; j < nc; j++)
{
std::sort (spairs + ci[j], spairs + ci[j+1]);
octave_idx_type l = -1;
octave_idx_type nzj = 0;
// Count.
for (octave_idx_type i = ci[j]; i < ci[j+1]; i++)
{
octave_idx_type k = spairs[i].first;
if (k != l)
{
l = k;
nzj++;
}
}
// Set column pointer.
xcidx (j+1) = xcidx (j) + nzj;
}
change_capacity (nzm > xcidx (nc) ? nzm : xcidx (nc));
octave_idx_type *rri = ridx ();
T *rrd = data ();
// Fill-in data.
for (octave_idx_type j = 0, jj = -1; j < nc; j++)
{
octave_quit ();
octave_idx_type l = -1;
if (sum_terms)
{
// Sum adjacent terms.
for (octave_idx_type i = ci[j]; i < ci[j+1]; i++)
{
octave_idx_type k = spairs[i].first;
if (k != l)
{
l = k;
rrd[++jj] = a(spairs[i].second);
rri[jj] = k;
}
else
rrd[jj] += a(spairs[i].second);
}
}
else
{
// Use the last one.
for (octave_idx_type i = ci[j]; i < ci[j+1]; i++)
{
octave_idx_type k = spairs[i].first;
if (k != l)
{
l = k;
rri[++jj] = k;
}
rrd[jj] = a(spairs[i].second);
}
}
}
maybe_compress (true);
}
}
/*
%!assert <*51880> (sparse (1:2, 2, 1:2, 2, 2), sparse ([0, 1; 0, 2]))
%!assert <*51880> (sparse (1:2, 1, 1:2, 2, 2), sparse ([1, 0; 2, 0]))
%!assert <*51880> (sparse (1:2, 2, 1:2, 2, 3), sparse ([0, 1, 0; 0, 2, 0]))
*/
template
OCTAVE_API
Sparse::Sparse (const Array& a)
: m_rep (nullptr), m_dimensions (a.dims ())
{
if (m_dimensions.ndims () > 2)
(*current_liboctave_error_handler)
("Sparse::Sparse (const Array&): dimension mismatch");
octave_idx_type nr = rows ();
octave_idx_type nc = cols ();
octave_idx_type len = a.numel ();
octave_idx_type new_nzmax = 0;
// First count the number of nonzero terms
for (octave_idx_type i = 0; i < len; i++)
if (a(i) != T ())
new_nzmax++;
m_rep = new typename Sparse::SparseRep (nr, nc, new_nzmax);
octave_idx_type ii = 0;
xcidx (0) = 0;
for (octave_idx_type j = 0; j < nc; j++)
{
for (octave_idx_type i = 0; i < nr; i++)
if (a.elem (i, j) != T ())
{
xdata (ii) = a.elem (i, j);
xridx (ii++) = i;
}
xcidx (j+1) = ii;
}
}
template
OCTAVE_API
Sparse::~Sparse (void)
{
if (--m_rep->m_count == 0)
delete m_rep;
}
template
Sparse&
Sparse::operator = (const Sparse& a)
{
if (this != &a)
{
if (--m_rep->m_count == 0)
delete m_rep;
m_rep = a.m_rep;
m_rep->m_count++;
m_dimensions = a.m_dimensions;
}
return *this;
}
template
OCTAVE_API
octave_idx_type
Sparse::compute_index (const Array& ra_idx) const
{
octave_idx_type n = m_dimensions.ndims ();
if (n <= 0 || n != ra_idx.numel ())
(*current_liboctave_error_handler)
("Sparse::compute_index: invalid ra_idxing operation");
octave_idx_type retval = -1;
retval = ra_idx(--n);
while (--n >= 0)
{
retval *= m_dimensions(n);
retval += ra_idx(n);
}
return retval;
}
template
OCTAVE_API
T
Sparse::range_error (const char *fcn, octave_idx_type n) const
{
(*current_liboctave_error_handler) ("%s (%" OCTAVE_IDX_TYPE_FORMAT "): "
"range error", fcn, n);
}
template
OCTAVE_API
T&
Sparse::range_error (const char *fcn, octave_idx_type n)
{
(*current_liboctave_error_handler) ("%s (%" OCTAVE_IDX_TYPE_FORMAT "): "
"range error", fcn, n);
}
template
OCTAVE_API
T
Sparse::range_error (const char *fcn, octave_idx_type i,
octave_idx_type j) const
{
(*current_liboctave_error_handler)
("%s (%" OCTAVE_IDX_TYPE_FORMAT ", %" OCTAVE_IDX_TYPE_FORMAT "): "
"range error", fcn, i, j);
}
template
OCTAVE_API
T&
Sparse::range_error (const char *fcn, octave_idx_type i, octave_idx_type j)
{
(*current_liboctave_error_handler)
("%s (%" OCTAVE_IDX_TYPE_FORMAT ", %" OCTAVE_IDX_TYPE_FORMAT "): "
"range error", fcn, i, j);
}
template
OCTAVE_API
T
Sparse::range_error (const char *fcn,
const Array& ra_idx) const
{
std::ostringstream buf;
buf << fcn << " (";
octave_idx_type n = ra_idx.numel ();
if (n > 0)
buf << ra_idx(0);
for (octave_idx_type i = 1; i < n; i++)
buf << ", " << ra_idx(i);
buf << "): range error";
std::string buf_str = buf.str ();
(*current_liboctave_error_handler) ("%s", buf_str.c_str ());
}
template
OCTAVE_API
T&
Sparse::range_error (const char *fcn, const Array& ra_idx)
{
std::ostringstream buf;
buf << fcn << " (";
octave_idx_type n = ra_idx.numel ();
if (n > 0)
buf << ra_idx(0);
for (octave_idx_type i = 1; i < n; i++)
buf << ", " << ra_idx(i);
buf << "): range error";
std::string buf_str = buf.str ();
(*current_liboctave_error_handler) ("%s", buf_str.c_str ());
}
template
OCTAVE_API
Sparse
Sparse::reshape (const dim_vector& new_dims) const
{
Sparse retval;
dim_vector dims2 = new_dims;
if (dims2.ndims () > 2)
{
(*current_liboctave_warning_with_id_handler)
("Octave:reshape-smashes-dims",
"reshape: sparse reshape to N-D array smashes dims");
for (octave_idx_type i = 2; i < dims2.ndims (); i++)
dims2(1) *= dims2(i);
dims2.resize (2);
}
if (m_dimensions != dims2)
{
if (m_dimensions.numel () == dims2.numel ())
{
octave_idx_type new_nnz = nnz ();
octave_idx_type new_nr = dims2 (0);
octave_idx_type new_nc = dims2 (1);
octave_idx_type old_nr = rows ();
octave_idx_type old_nc = cols ();
retval = Sparse (new_nr, new_nc, new_nnz);
octave_idx_type kk = 0;
retval.xcidx (0) = 0;
// Quotient and remainder of i * old_nr divided by new_nr.
// Track them individually to avoid overflow (bug #42850).
octave_idx_type i_old_qu = 0;
octave_idx_type i_old_rm = static_cast (-old_nr);
for (octave_idx_type i = 0; i < old_nc; i++)
{
i_old_rm += old_nr;
if (i_old_rm >= new_nr)
{
i_old_qu += i_old_rm / new_nr;
i_old_rm = i_old_rm % new_nr;
}
for (octave_idx_type j = cidx (i); j < cidx (i+1); j++)
{
octave_idx_type ii, jj;
ii = (i_old_rm + ridx (j)) % new_nr;
jj = i_old_qu + (i_old_rm + ridx (j)) / new_nr;
// Original calculation subject to overflow
// ii = (i*old_nr + ridx (j)) % new_nr
// jj = (i*old_nr + ridx (j)) / new_nr
for (octave_idx_type k = kk; k < jj; k++)
retval.xcidx (k+1) = j;
kk = jj;
retval.xdata (j) = data (j);
retval.xridx (j) = ii;
}
}
for (octave_idx_type k = kk; k < new_nc; k++)
retval.xcidx (k+1) = new_nnz;
}
else
{
std::string dimensions_str = m_dimensions.str ();
std::string new_dims_str = new_dims.str ();
(*current_liboctave_error_handler)
("reshape: can't reshape %s array to %s array",
dimensions_str.c_str (), new_dims_str.c_str ());
}
}
else
retval = *this;
return retval;
}
template
OCTAVE_API
Sparse
Sparse::permute (const Array& perm_vec, bool) const
{
// The only valid permutations of a sparse array are [1, 2] and [2, 1].
bool fail = false;
bool trans = false;
if (perm_vec.numel () == 2)
{
if (perm_vec(0) == 0 && perm_vec(1) == 1)
/* do nothing */;
else if (perm_vec(0) == 1 && perm_vec(1) == 0)
trans = true;
else
fail = true;
}
else
fail = true;
if (fail)
(*current_liboctave_error_handler)
("permutation vector contains an invalid element");
return trans ? this->transpose () : *this;
}
template
OCTAVE_API
void
Sparse::resize1 (octave_idx_type n)
{
octave_idx_type nr = rows ();
octave_idx_type nc = cols ();
if (nr == 0)
resize (1, std::max (nc, n));
else if (nc == 0)
resize (nr, (n + nr - 1) / nr); // Ain't it wicked?
else if (nr == 1)
resize (1, n);
else if (nc == 1)
resize (n, 1);
else
octave::err_invalid_resize ();
}
template
OCTAVE_API
void
Sparse::resize (const dim_vector& dv)
{
octave_idx_type n = dv.ndims ();
if (n != 2)
(*current_liboctave_error_handler) ("sparse array must be 2-D");
resize (dv(0), dv(1));
}
template
OCTAVE_API
void
Sparse::resize (octave_idx_type r, octave_idx_type c)
{
if (r < 0 || c < 0)
(*current_liboctave_error_handler) ("can't resize to negative dimension");
if (r == dim1 () && c == dim2 ())
return;
// This wouldn't be necessary for r >= rows () if m_nrows wasn't part of the
// Sparse rep. It is not good for anything in there.
make_unique ();
if (r < rows ())
{
octave_idx_type i = 0;
octave_idx_type k = 0;
for (octave_idx_type j = 1; j <= m_rep->m_ncols; j++)
{
octave_idx_type u = xcidx (j);
for (; i < u; i++)
if (xridx (i) < r)
{
xdata (k) = xdata (i);
xridx (k++) = xridx (i);
}
xcidx (j) = k;
}
}
m_rep->m_nrows = m_dimensions(0) = r;
if (c != m_rep->m_ncols)
{
octave_idx_type *new_cidx = new octave_idx_type [c+1];
std::copy_n (m_rep->m_cidx, std::min (c, m_rep->m_ncols) + 1, new_cidx);
m_rep->m_cidx_deleter (m_rep->m_cidx);
m_rep->m_cidx = new_cidx;
m_rep->m_cidx_deleter = delete_array;
if (c > m_rep->m_ncols)
std::fill_n (m_rep->m_cidx + m_rep->m_ncols + 1, c - m_rep->m_ncols,
m_rep->m_cidx[m_rep->m_ncols]);
}
m_rep->m_ncols = m_dimensions(1) = c;
m_rep->change_length (m_rep->nnz ());
}
template
OCTAVE_API
Sparse&
Sparse::insert (const Sparse& a, octave_idx_type r, octave_idx_type c)
{
octave_idx_type a_rows = a.rows ();
octave_idx_type a_cols = a.cols ();
octave_idx_type nr = rows ();
octave_idx_type nc = cols ();
if (r < 0 || r + a_rows > rows () || c < 0 || c + a_cols > cols ())
(*current_liboctave_error_handler) ("range error for insert");
// First count the number of elements in the final array
octave_idx_type nel = cidx (c) + a.nnz ();
if (c + a_cols < nc)
nel += cidx (nc) - cidx (c + a_cols);
for (octave_idx_type i = c; i < c + a_cols; i++)
for (octave_idx_type j = cidx (i); j < cidx (i+1); j++)
if (ridx (j) < r || ridx (j) >= r + a_rows)
nel++;
Sparse tmp (*this);
--m_rep->m_count;
m_rep = new typename Sparse::SparseRep (nr, nc, nel);
for (octave_idx_type i = 0; i < tmp.cidx (c); i++)
{
data (i) = tmp.data (i);
ridx (i) = tmp.ridx (i);
}
for (octave_idx_type i = 0; i < c + 1; i++)
cidx (i) = tmp.cidx (i);
octave_idx_type ii = cidx (c);
for (octave_idx_type i = c; i < c + a_cols; i++)
{
octave_quit ();
for (octave_idx_type j = tmp.cidx (i); j < tmp.cidx (i+1); j++)
if (tmp.ridx (j) < r)
{
data (ii) = tmp.data (j);
ridx (ii++) = tmp.ridx (j);
}
octave_quit ();
for (octave_idx_type j = a.cidx (i-c); j < a.cidx (i-c+1); j++)
{
data (ii) = a.data (j);
ridx (ii++) = r + a.ridx (j);
}
octave_quit ();
for (octave_idx_type j = tmp.cidx (i); j < tmp.cidx (i+1); j++)
if (tmp.ridx (j) >= r + a_rows)
{
data (ii) = tmp.data (j);
ridx (ii++) = tmp.ridx (j);
}
cidx (i+1) = ii;
}
for (octave_idx_type i = c + a_cols; i < nc; i++)
{
for (octave_idx_type j = tmp.cidx (i); j < tmp.cidx (i+1); j++)
{
data (ii) = tmp.data (j);
ridx (ii++) = tmp.ridx (j);
}
cidx (i+1) = ii;
}
return *this;
}
template
OCTAVE_API
Sparse&
Sparse::insert (const Sparse& a, const Array& ra_idx)
{
if (ra_idx.numel () != 2)
(*current_liboctave_error_handler) ("range error for insert");
return insert (a, ra_idx(0), ra_idx(1));
}
template
OCTAVE_API
Sparse
Sparse::transpose (void) const
{
assert (ndims () == 2);
octave_idx_type nr = rows ();
octave_idx_type nc = cols ();
octave_idx_type nz = nnz ();
Sparse retval (nc, nr, nz);
for (octave_idx_type i = 0; i < nz; i++)
retval.xcidx (ridx (i) + 1)++;
// retval.xcidx[1:nr] holds the row degrees for rows 0:(nr-1)
nz = 0;
for (octave_idx_type i = 1; i <= nr; i++)
{
const octave_idx_type tmp = retval.xcidx (i);
retval.xcidx (i) = nz;
nz += tmp;
}
// retval.xcidx[1:nr] holds row entry *start* offsets for rows 0:(nr-1)
for (octave_idx_type j = 0; j < nc; j++)
for (octave_idx_type k = cidx (j); k < cidx (j+1); k++)
{
octave_idx_type q = retval.xcidx (ridx (k) + 1)++;
retval.xridx (q) = j;
retval.xdata (q) = data (k);
}
assert (nnz () == retval.xcidx (nr));
// retval.xcidx[1:nr] holds row entry *end* offsets for rows 0:(nr-1)
// and retval.xcidx[0:(nr-1)] holds their row entry *start* offsets
return retval;
}
// Lower bound lookup. Could also use octave_sort, but that has upper bound
// semantics, so requires some manipulation to set right. Uses a plain loop
// for small columns.
static
octave_idx_type
lblookup (const octave_idx_type *ridx, octave_idx_type nr,
octave_idx_type ri)
{
if (nr <= 8)
{
octave_idx_type l;
for (l = 0; l < nr; l++)
if (ridx[l] >= ri)
break;
return l;
}
else
return std::lower_bound (ridx, ridx + nr, ri) - ridx;
}
template
OCTAVE_API
void
Sparse::delete_elements (const octave::idx_vector& idx)
{
Sparse retval;
assert (ndims () == 2);
octave_idx_type nr = dim1 ();
octave_idx_type nc = dim2 ();
octave_idx_type nz = nnz ();
octave_idx_type nel = numel (); // Can throw.
const dim_vector idx_dims = idx.orig_dimensions ();
if (idx.extent (nel) > nel)
octave::err_del_index_out_of_range (true, idx.extent (nel), nel);
if (nc == 1)
{
// Sparse column vector.
const Sparse tmp = *this; // constant copy to prevent COW.
octave_idx_type lb, ub;
if (idx.is_cont_range (nel, lb, ub))
{
// Special-case a contiguous range.
// Look-up indices first.
octave_idx_type li = lblookup (tmp.ridx (), nz, lb);
octave_idx_type ui = lblookup (tmp.ridx (), nz, ub);
// Copy data and adjust indices.
octave_idx_type nz_new = nz - (ui - li);
*this = Sparse (nr - (ub - lb), 1, nz_new);
std::copy_n (tmp.data (), li, data ());
std::copy_n (tmp.ridx (), li, xridx ());
std::copy (tmp.data () + ui, tmp.data () + nz, xdata () + li);
mx_inline_sub (nz - ui, xridx () + li, tmp.ridx () + ui, ub - lb);
xcidx (1) = nz_new;
}
else
{
OCTAVE_LOCAL_BUFFER (octave_idx_type, ridx_new, nz);
OCTAVE_LOCAL_BUFFER (T, data_new, nz);
octave::idx_vector sidx = idx.sorted (true);
const octave_idx_type *sj = sidx.raw ();
octave_idx_type sl = sidx.length (nel);
octave_idx_type nz_new = 0;
octave_idx_type j = 0;
for (octave_idx_type i = 0; i < nz; i++)
{
octave_idx_type r = tmp.ridx (i);
for (; j < sl && sj[j] < r; j++) ;
if (j == sl || sj[j] > r)
{
data_new[nz_new] = tmp.data (i);
ridx_new[nz_new++] = r - j;
}
}
*this = Sparse (nr - sl, 1, nz_new);
std::copy_n (ridx_new, nz_new, ridx ());
std::copy_n (data_new, nz_new, xdata ());
xcidx (1) = nz_new;
}
}
else if (nr == 1)
{
// Sparse row vector.
octave_idx_type lb, ub;
if (idx.is_cont_range (nc, lb, ub))
{
const Sparse tmp = *this;
octave_idx_type lbi = tmp.cidx (lb);
octave_idx_type ubi = tmp.cidx (ub);
octave_idx_type new_nz = nz - (ubi - lbi);
*this = Sparse (1, nc - (ub - lb), new_nz);
std::copy_n (tmp.data (), lbi, data ());
std::copy (tmp.data () + ubi, tmp.data () + nz, xdata () + lbi);
std::fill_n (ridx (), new_nz, static_cast (0));
std::copy_n (tmp.cidx () + 1, lb, cidx () + 1);
mx_inline_sub (nc - ub, xcidx () + 1, tmp.cidx () + ub + 1,
ubi - lbi);
}
else
*this = index (idx.complement (nc));
}
else if (idx.length (nel) != 0)
{
if (idx.is_colon_equiv (nel))
*this = Sparse ();
else
{
*this = index (octave::idx_vector::colon);
delete_elements (idx);
*this = transpose (); // We want a row vector.
}
}
}
template
OCTAVE_API
void
Sparse::delete_elements (const octave::idx_vector& idx_i, const octave::idx_vector& idx_j)
{
assert (ndims () == 2);
octave_idx_type nr = dim1 ();
octave_idx_type nc = dim2 ();
octave_idx_type nz = nnz ();
if (idx_i.is_colon ())
{
// Deleting columns.
octave_idx_type lb, ub;
if (idx_j.extent (nc) > nc)
octave::err_del_index_out_of_range (false, idx_j.extent (nc), nc);
else if (idx_j.is_cont_range (nc, lb, ub))
{
if (lb == 0 && ub == nc)
{
// Delete all rows and columns.
*this = Sparse (nr, 0);
}
else if (nz == 0)
{
// No elements to preserve; adjust dimensions.
*this = Sparse (nr, nc - (ub - lb));
}
else
{
const Sparse tmp = *this;
octave_idx_type lbi = tmp.cidx (lb);
octave_idx_type ubi = tmp.cidx (ub);
octave_idx_type new_nz = nz - (ubi - lbi);
*this = Sparse (nr, nc - (ub - lb), new_nz);
std::copy_n (tmp.data (), lbi, data ());
std::copy_n (tmp.ridx (), lbi, ridx ());
std::copy (tmp.data () + ubi, tmp.data () + nz, xdata () + lbi);
std::copy (tmp.ridx () + ubi, tmp.ridx () + nz, xridx () + lbi);
std::copy_n (tmp.cidx () + 1, lb, cidx () + 1);
mx_inline_sub (nc - ub, xcidx () + lb + 1,
tmp.cidx () + ub + 1, ubi - lbi);
}
}
else
*this = index (idx_i, idx_j.complement (nc));
}
else if (idx_j.is_colon ())
{
// Deleting rows.
octave_idx_type lb, ub;
if (idx_i.extent (nr) > nr)
octave::err_del_index_out_of_range (false, idx_i.extent (nr), nr);
else if (idx_i.is_cont_range (nr, lb, ub))
{
if (lb == 0 && ub == nr)
{
// Delete all rows and columns.
*this = Sparse (0, nc);
}
else if (nz == 0)
{
// No elements to preserve; adjust dimensions.
*this = Sparse (nr - (ub - lb), nc);
}
else
{
// This is more memory-efficient than the approach below.
const Sparse tmpl = index (octave::idx_vector (0, lb), idx_j);
const Sparse tmpu = index (octave::idx_vector (ub, nr), idx_j);
*this = Sparse (nr - (ub - lb), nc,
tmpl.nnz () + tmpu.nnz ());
for (octave_idx_type j = 0, k = 0; j < nc; j++)
{
for (octave_idx_type i = tmpl.cidx (j); i < tmpl.cidx (j+1);
i++)
{
xdata (k) = tmpl.data (i);
xridx (k++) = tmpl.ridx (i);
}
for (octave_idx_type i = tmpu.cidx (j); i < tmpu.cidx (j+1);
i++)
{
xdata (k) = tmpu.data (i);
xridx (k++) = tmpu.ridx (i) + lb;
}
xcidx (j+1) = k;
}
}
}
else
{
// This is done by transposing, deleting columns, then transposing
// again.
Sparse tmp = transpose ();
tmp.delete_elements (idx_j, idx_i);
*this = tmp.transpose ();
}
}
else
{
// Empty assignment (no elements to delete) is OK if there is at
// least one zero-length index and at most one other index that is
// non-colon (or equivalent) index. Since we only have two
// indices, we just need to check that we have at least one zero
// length index. Matlab considers "[]" to be an empty index but
// not "false". We accept both.
bool empty_assignment
= (idx_i.length (nr) == 0 || idx_j.length (nc) == 0);
if (! empty_assignment)
(*current_liboctave_error_handler)
("a null assignment can only have one non-colon index");
}
}
template
OCTAVE_API
void
Sparse::delete_elements (int dim, const octave::idx_vector& idx)
{
if (dim == 0)
delete_elements (idx, octave::idx_vector::colon);
else if (dim == 1)
delete_elements (octave::idx_vector::colon, idx);
else
(*current_liboctave_error_handler) ("invalid dimension in delete_elements");
}
template
OCTAVE_API
Sparse
Sparse::index (const octave::idx_vector& idx, bool resize_ok) const
{
Sparse retval;
assert (ndims () == 2);
octave_idx_type nr = dim1 ();
octave_idx_type nc = dim2 ();
octave_idx_type nz = nnz ();
octave_idx_type nel = numel (); // Can throw.
const dim_vector idx_dims = idx.orig_dimensions ().redim (2);
if (idx.is_colon ())
{
if (nc == 1)
retval = *this;
else
{
// Fast magic colon processing.
retval = Sparse (nel, 1, nz);
for (octave_idx_type i = 0; i < nc; i++)
{
for (octave_idx_type j = cidx (i); j < cidx (i+1); j++)
{
retval.xdata (j) = data (j);
retval.xridx (j) = ridx (j) + i * nr;
}
}
retval.xcidx (0) = 0;
retval.xcidx (1) = nz;
}
}
else if (idx.extent (nel) > nel)
{
if (! resize_ok)
octave::err_index_out_of_range (1, 1, idx.extent (nel), nel, dims ());
// resize_ok is completely handled here.
octave_idx_type ext = idx.extent (nel);
Sparse tmp = *this;
tmp.resize1 (ext);
retval = tmp.index (idx);
}
else if (nr == 1 && nc == 1)
{
// You have to be pretty sick to get to this bit of code,
// since you have a scalar stored as a sparse matrix, and
// then want to make a dense matrix with sparse representation.
// Ok, we'll do it, but you deserve what you get!!
retval = (Sparse (idx_dims(0), idx_dims(1), nz ? data (0) : T ()));
}
else if (nc == 1)
{
// Sparse column vector.
octave_idx_type lb, ub;
if (idx.is_scalar ())
{
// Scalar index - just a binary lookup.
octave_idx_type i = lblookup (ridx (), nz, idx(0));
if (i < nz && ridx (i) == idx(0))
retval = Sparse (1, 1, data (i));
else
retval = Sparse (1, 1);
}
else if (idx.is_cont_range (nel, lb, ub))
{
// Special-case a contiguous range.
// Look-up indices first.
octave_idx_type li = lblookup (ridx (), nz, lb);
octave_idx_type ui = lblookup (ridx (), nz, ub);
// Copy data and adjust indices.
octave_idx_type nz_new = ui - li;
retval = Sparse (ub - lb, 1, nz_new);
std::copy_n (data () + li, nz_new, retval.data ());
mx_inline_sub (nz_new, retval.xridx (), ridx () + li, lb);
retval.xcidx (1) = nz_new;
}
else if (idx.is_permutation (nel) && idx.isvector ())
{
if (idx.is_range () && idx.increment () == -1)
{
retval = Sparse (nr, 1, nz);
for (octave_idx_type j = 0; j < nz; j++)
retval.ridx (j) = nr - ridx (nz - j - 1) - 1;
std::copy_n (cidx (), 2, retval.cidx ());
std::reverse_copy (data (), data () + nz, retval.data ());
}
else
{
Array tmp = array_value ();
tmp = tmp.index (idx);
retval = Sparse (tmp);
}
}
else
{
// If indexing a sparse column vector by a vector, the result is a
// sparse column vector, otherwise it inherits the shape of index.
// Vector transpose is cheap, so do it right here.
Array tmp_idx = idx.as_array ().as_matrix ();
const Array idxa = (idx_dims(0) == 1
? tmp_idx.transpose ()
: tmp_idx);
octave_idx_type new_nr = idxa.rows ();
octave_idx_type new_nc = idxa.cols ();
// Lookup.
// FIXME: Could specialize for sorted idx?
Array lidx (dim_vector (new_nr, new_nc));
for (octave_idx_type i = 0; i < new_nr*new_nc; i++)
lidx.xelem (i) = lblookup (ridx (), nz, idxa(i));
// Count matches.
retval = Sparse (idxa.rows (), idxa.cols ());
for (octave_idx_type j = 0; j < new_nc; j++)
{
octave_idx_type nzj = 0;
for (octave_idx_type i = 0; i < new_nr; i++)
{
octave_idx_type l = lidx.xelem (i, j);
if (l < nz && ridx (l) == idxa(i, j))
nzj++;
else
lidx.xelem (i, j) = nz;
}
retval.xcidx (j+1) = retval.xcidx (j) + nzj;
}
retval.change_capacity (retval.xcidx (new_nc));
// Copy data and set row indices.
octave_idx_type k = 0;
for (octave_idx_type j = 0; j < new_nc; j++)
for (octave_idx_type i = 0; i < new_nr; i++)
{
octave_idx_type l = lidx.xelem (i, j);
if (l < nz)
{
retval.data (k) = data (l);
retval.xridx (k++) = i;
}
}
}
}
else if (nr == 1)
{
octave_idx_type lb, ub;
if (idx.is_scalar ())
retval = Sparse (1, 1, elem (0, idx(0)));
else if (idx.is_cont_range (nel, lb, ub))
{
// Special-case a contiguous range.
octave_idx_type lbi = cidx (lb);
octave_idx_type ubi = cidx (ub);
octave_idx_type new_nz = ubi - lbi;
retval = Sparse (1, ub - lb, new_nz);
std::copy_n (data () + lbi, new_nz, retval.data ());
std::fill_n (retval.ridx (), new_nz, static_cast (0));
mx_inline_sub (ub - lb + 1, retval.cidx (), cidx () + lb, lbi);
}
else
{
// Sparse row vectors occupy O(nr) storage anyway, so let's just
// convert the matrix to full, index, and sparsify the result.
retval = Sparse (array_value ().index (idx));
}
}
else
{
if (nr != 0 && idx.is_scalar ())
retval = Sparse (1, 1, elem (idx(0) % nr, idx(0) / nr));
else
{
// Indexing a non-vector sparse matrix by linear indexing.
// I suppose this is rare (and it may easily overflow), so let's take
// the easy way, and reshape first to column vector, which is already
// handled above.
retval = index (octave::idx_vector::colon).index (idx);
// In this case we're supposed to always inherit the shape, but
// column(row) doesn't do it, so we'll do it instead.
if (idx_dims(0) == 1 && idx_dims(1) != 1)
retval = retval.transpose ();
}
}
return retval;
}
template
OCTAVE_API
Sparse
Sparse::index (const octave::idx_vector& idx_i, const octave::idx_vector& idx_j,
bool resize_ok) const
{
Sparse retval;
assert (ndims () == 2);
octave_idx_type nr = dim1 ();
octave_idx_type nc = dim2 ();
octave_idx_type n = idx_i.length (nr);
octave_idx_type m = idx_j.length (nc);
octave_idx_type lb, ub;
if (idx_i.extent (nr) > nr || idx_j.extent (nc) > nc)
{
// resize_ok is completely handled here.
if (resize_ok)
{
octave_idx_type ext_i = idx_i.extent (nr);
octave_idx_type ext_j = idx_j.extent (nc);
Sparse tmp = *this;
tmp.resize (ext_i, ext_j);
retval = tmp.index (idx_i, idx_j);
}
else if (idx_i.extent (nr) > nr)
octave::err_index_out_of_range (2, 1, idx_i.extent (nr), nr, dims ());
else
octave::err_index_out_of_range (2, 2, idx_j.extent (nc), nc, dims ());
}
else if (nr == 1 && nc == 1)
{
// Scalars stored as sparse matrices occupy more memory than
// a scalar, so let's just convert the matrix to full, index,
// and sparsify the result.
retval = Sparse (array_value ().index (idx_i, idx_j));
}
else if (idx_i.is_colon ())
{
// Great, we're just manipulating columns. This is going to be quite
// efficient, because the columns can stay compressed as they are.
if (idx_j.is_colon ())
retval = *this; // Shallow copy.
else if (idx_j.is_cont_range (nc, lb, ub))
{
// Special-case a contiguous range.
octave_idx_type lbi = cidx (lb);
octave_idx_type ubi = cidx (ub);
octave_idx_type new_nz = ubi - lbi;
retval = Sparse (nr, ub - lb, new_nz);
std::copy_n (data () + lbi, new_nz, retval.data ());
std::copy_n (ridx () + lbi, new_nz, retval.ridx ());
mx_inline_sub (ub - lb + 1, retval.cidx (), cidx () + lb, lbi);
}
else
{
// Count new nonzero elements.
retval = Sparse (nr, m);
for (octave_idx_type j = 0; j < m; j++)
{
octave_idx_type jj = idx_j(j);
retval.xcidx (j+1) = retval.xcidx (j) + (cidx (jj+1) - cidx (jj));
}
retval.change_capacity (retval.xcidx (m));
// Copy data & indices.
for (octave_idx_type j = 0; j < m; j++)
{
octave_idx_type ljj = cidx (idx_j(j));
octave_idx_type lj = retval.xcidx (j);
octave_idx_type nzj = retval.xcidx (j+1) - lj;
std::copy_n (data () + ljj, nzj, retval.data () + lj);
std::copy_n (ridx () + ljj, nzj, retval.ridx () + lj);
}
}
}
else if (nc == 1 && idx_j.is_colon_equiv (nc) && idx_i.isvector ())
{
// It's actually vector indexing. The 1D index is specialized for that.
retval = index (idx_i);
// If nr == 1 then the vector indexing will return a column vector!!
if (nr == 1)
retval.transpose ();
}
else if (idx_i.is_scalar ())
{
octave_idx_type ii = idx_i(0);
retval = Sparse (1, m);
OCTAVE_LOCAL_BUFFER (octave_idx_type, ij, m);
for (octave_idx_type j = 0; j < m; j++)
{
octave_quit ();
octave_idx_type jj = idx_j(j);
octave_idx_type lj = cidx (jj);
octave_idx_type nzj = cidx (jj+1) - cidx (jj);
// Scalar index - just a binary lookup.
octave_idx_type i = lblookup (ridx () + lj, nzj, ii);
if (i < nzj && ridx (i+lj) == ii)
{
ij[j] = i + lj;
retval.xcidx (j+1) = retval.xcidx (j) + 1;
}
else
retval.xcidx (j+1) = retval.xcidx (j);
}
retval.change_capacity (retval.xcidx (m));
// Copy data, adjust row indices.
for (octave_idx_type j = 0; j < m; j++)
{
octave_idx_type i = retval.xcidx (j);
if (retval.xcidx (j+1) > i)
{
retval.xridx (i) = 0;
retval.xdata (i) = data (ij[j]);
}
}
}
else if (idx_i.is_cont_range (nr, lb, ub))
{
retval = Sparse (n, m);
OCTAVE_LOCAL_BUFFER (octave_idx_type, li, m);
OCTAVE_LOCAL_BUFFER (octave_idx_type, ui, m);
for (octave_idx_type j = 0; j < m; j++)
{
octave_quit ();
octave_idx_type jj = idx_j(j);
octave_idx_type lj = cidx (jj);
octave_idx_type nzj = cidx (jj+1) - cidx (jj);
octave_idx_type lij, uij;
// Lookup indices.
li[j] = lij = lblookup (ridx () + lj, nzj, lb) + lj;
ui[j] = uij = lblookup (ridx () + lj, nzj, ub) + lj;
retval.xcidx (j+1) = retval.xcidx (j) + ui[j] - li[j];
}
retval.change_capacity (retval.xcidx (m));
// Copy data, adjust row indices.
for (octave_idx_type j = 0, k = 0; j < m; j++)
{
octave_quit ();
for (octave_idx_type i = li[j]; i < ui[j]; i++)
{
retval.xdata (k) = data (i);
retval.xridx (k++) = ridx (i) - lb;
}
}
}
else if (idx_i.is_permutation (nr))
{
// The columns preserve their length, just need to renumber and sort them.
// Count new nonzero elements.
retval = Sparse (nr, m);
for (octave_idx_type j = 0; j < m; j++)
{
octave_idx_type jj = idx_j(j);
retval.xcidx (j+1) = retval.xcidx (j) + (cidx (jj+1) - cidx (jj));
}
retval.change_capacity (retval.xcidx (m));
octave_quit ();
if (idx_i.is_range () && idx_i.increment () == -1)
{
// It's nr:-1:1. Just flip all columns.
for (octave_idx_type j = 0; j < m; j++)
{
octave_quit ();
octave_idx_type jj = idx_j(j);
octave_idx_type lj = cidx (jj);
octave_idx_type nzj = cidx (jj+1) - cidx (jj);
octave_idx_type li = retval.xcidx (j);
octave_idx_type uj = lj + nzj - 1;
for (octave_idx_type i = 0; i < nzj; i++)
{
retval.xdata (li + i) = data (uj - i); // Copy in reverse order.
retval.xridx (li + i) = nr - 1 - ridx (uj - i); // Ditto with transform.
}
}
}
else
{
// Get inverse permutation.
octave::idx_vector idx_iinv = idx_i.inverse_permutation (nr);
const octave_idx_type *iinv = idx_iinv.raw ();
// Scatter buffer.
OCTAVE_LOCAL_BUFFER (T, scb, nr);
octave_idx_type *rri = retval.ridx ();
for (octave_idx_type j = 0; j < m; j++)
{
octave_quit ();
octave_idx_type jj = idx_j(j);
octave_idx_type lj = cidx (jj);
octave_idx_type nzj = cidx (jj+1) - cidx (jj);
octave_idx_type li = retval.xcidx (j);
// Scatter the column, transform indices.
for (octave_idx_type i = 0; i < nzj; i++)
scb[rri[li + i] = iinv[ridx (lj + i)]] = data (lj + i);
octave_quit ();
// Sort them.
std::sort (rri + li, rri + li + nzj);
// Gather.
for (octave_idx_type i = 0; i < nzj; i++)
retval.xdata (li + i) = scb[rri[li + i]];
}
}
}
else if (idx_j.is_colon ())
{
// This requires uncompressing columns, which is generally costly,
// so we rely on the efficient transpose to handle this.
// It may still make sense to optimize some cases here.
retval = transpose ();
retval = retval.index (octave::idx_vector::colon, idx_i);
retval = retval.transpose ();
}
else
{
// A(I, J) is decomposed into A(:, J)(I, :).
retval = index (octave::idx_vector::colon, idx_j);
retval = retval.index (idx_i, octave::idx_vector::colon);
}
return retval;
}
template
OCTAVE_API
void
Sparse::assign (const octave::idx_vector& idx, const Sparse& rhs)
{
Sparse retval;
assert (ndims () == 2);
octave_idx_type nr = dim1 ();
octave_idx_type nc = dim2 ();
octave_idx_type nz = nnz ();
octave_idx_type n = numel (); // Can throw.
octave_idx_type rhl = rhs.numel ();
if (idx.length (n) == rhl)
{
if (rhl == 0)
return;
octave_idx_type nx = idx.extent (n);
// Try to resize first if necessary.
if (nx != n)
{
resize1 (nx);
n = numel ();
nr = rows ();
nc = cols ();
// nz is preserved.
}
if (idx.is_colon ())
{
*this = rhs.reshape (m_dimensions);
}
else if (nc == 1 && rhs.cols () == 1)
{
// Sparse column vector to sparse column vector assignment.
octave_idx_type lb, ub;
if (idx.is_cont_range (nr, lb, ub))
{
// Special-case a contiguous range.
// Look-up indices first.
octave_idx_type li = lblookup (ridx (), nz, lb);
octave_idx_type ui = lblookup (ridx (), nz, ub);
octave_idx_type rnz = rhs.nnz ();
octave_idx_type new_nz = nz - (ui - li) + rnz;
if (new_nz >= nz && new_nz <= nzmax ())
{
// Adding/overwriting elements, enough capacity allocated.
if (new_nz > nz)
{
// Make room first.
std::copy_backward (data () + ui, data () + nz,
data () + nz + rnz);
std::copy_backward (ridx () + ui, ridx () + nz,
ridx () + nz + rnz);
}
// Copy data and adjust indices from rhs.
std::copy_n (rhs.data (), rnz, data () + li);
mx_inline_add (rnz, ridx () + li, rhs.ridx (), lb);
}
else
{
// Clearing elements or exceeding capacity, allocate afresh
// and paste pieces.
const Sparse tmp = *this;
*this = Sparse (nr, 1, new_nz);
// Head ...
std::copy_n (tmp.data (), li, data ());
std::copy_n (tmp.ridx (), li, ridx ());
// new stuff ...
std::copy_n (rhs.data (), rnz, data () + li);
mx_inline_add (rnz, ridx () + li, rhs.ridx (), lb);
// ...tail
std::copy (tmp.data () + ui, tmp.data () + nz,
data () + li + rnz);
std::copy (tmp.ridx () + ui, tmp.ridx () + nz,
ridx () + li + rnz);
}
cidx (1) = new_nz;
}
else if (idx.is_range () && idx.increment () == -1)
{
// It's s(u:-1:l) = r. Reverse the assignment.
assign (idx.sorted (), rhs.index (octave::idx_vector (rhl - 1, 0, -1)));
}
else if (idx.is_permutation (n))
{
*this = rhs.index (idx.inverse_permutation (n));
}
else if (rhs.nnz () == 0)
{
// Elements are being zeroed.
octave_idx_type *ri = ridx ();
for (octave_idx_type i = 0; i < rhl; i++)
{
octave_idx_type iidx = idx(i);
octave_idx_type li = lblookup (ri, nz, iidx);
if (li != nz && ri[li] == iidx)
xdata (li) = T ();
}
maybe_compress (true);
}
else
{
const Sparse tmp = *this;
octave_idx_type new_nz = nz + rhl;
// Disassembly our matrix...
Array new_ri (dim_vector (new_nz, 1));
Array new_data (dim_vector (new_nz, 1));
std::copy_n (tmp.ridx (), nz, new_ri.fortran_vec ());
std::copy_n (tmp.data (), nz, new_data.fortran_vec ());
// ... insert new data (densified) ...
idx.copy_data (new_ri.fortran_vec () + nz);
new_data.assign (octave::idx_vector (nz, new_nz), rhs.array_value ());
// ... reassembly.
*this = Sparse (new_data, new_ri, 0, nr, nc, false);
}
}
else
{
dim_vector save_dims = m_dimensions;
*this = index (octave::idx_vector::colon);
assign (idx, rhs.index (octave::idx_vector::colon));
*this = reshape (save_dims);
}
}
else if (rhl == 1)
{
rhl = idx.length (n);
if (rhs.nnz () != 0)
assign (idx, Sparse (rhl, 1, rhs.data (0)));
else
assign (idx, Sparse (rhl, 1));
}
else
octave::err_nonconformant ("=", dim_vector(idx.length (n), 1), rhs.dims());
}
template
OCTAVE_API
void
Sparse::assign (const octave::idx_vector& idx, const T& rhs)
{
// FIXME: Converting the RHS and forwarding to the sparse matrix
// assignment function is simpler, but it might be good to have a
// specialization...
assign (idx, Sparse (1, 1, rhs));
}
template
OCTAVE_API
void
Sparse::assign (const octave::idx_vector& idx_i,
const octave::idx_vector& idx_j, const Sparse& rhs)
{
Sparse retval;
assert (ndims () == 2);
octave_idx_type nr = dim1 ();
octave_idx_type nc = dim2 ();
octave_idx_type nz = nnz ();
octave_idx_type n = rhs.rows ();
octave_idx_type m = rhs.columns ();
// FIXME: this should probably be written more like the
// Array::assign function...
bool orig_zero_by_zero = (nr == 0 && nc == 0);
if (orig_zero_by_zero || (idx_i.length (nr) == n && idx_j.length (nc) == m))
{
octave_idx_type nrx;
octave_idx_type ncx;
if (orig_zero_by_zero)
{
if (idx_i.is_colon ())
{
nrx = n;
if (idx_j.is_colon ())
ncx = m;
else
ncx = idx_j.extent (nc);
}
else if (idx_j.is_colon ())
{
nrx = idx_i.extent (nr);
ncx = m;
}
else
{
nrx = idx_i.extent (nr);
ncx = idx_j.extent (nc);
}
}
else
{
nrx = idx_i.extent (nr);
ncx = idx_j.extent (nc);
}
// Try to resize first if necessary.
if (nrx != nr || ncx != nc)
{
resize (nrx, ncx);
nr = rows ();
nc = cols ();
// nz is preserved.
}
if (n == 0 || m == 0)
return;
if (idx_i.is_colon ())
{
octave_idx_type lb, ub;
// Great, we're just manipulating columns. This is going to be quite
// efficient, because the columns can stay compressed as they are.
if (idx_j.is_colon ())
*this = rhs; // Shallow copy.
else if (idx_j.is_cont_range (nc, lb, ub))
{
// Special-case a contiguous range.
octave_idx_type li = cidx (lb);
octave_idx_type ui = cidx (ub);
octave_idx_type rnz = rhs.nnz ();
octave_idx_type new_nz = nz - (ui - li) + rnz;
if (new_nz >= nz && new_nz <= nzmax ())
{
// Adding/overwriting elements, enough capacity allocated.
if (new_nz > nz)
{
// Make room first.
std::copy_backward (data () + ui, data () + nz,
data () + new_nz);
std::copy_backward (ridx () + ui, ridx () + nz,
ridx () + new_nz);
mx_inline_add2 (nc - ub, cidx () + ub + 1, new_nz - nz);
}
// Copy data and indices from rhs.
std::copy_n (rhs.data (), rnz, data () + li);
std::copy_n (rhs.ridx (), rnz, ridx () + li);
mx_inline_add (ub - lb, cidx () + lb + 1, rhs.cidx () + 1,
li);
assert (nnz () == new_nz);
}
else
{
// Clearing elements or exceeding capacity, allocate afresh
// and paste pieces.
const Sparse tmp = *this;
*this = Sparse (nr, nc, new_nz);
// Head...
std::copy_n (tmp.data (), li, data ());
std::copy_n (tmp.ridx (), li, ridx ());
std::copy_n (tmp.cidx () + 1, lb, cidx () + 1);
// new stuff...
std::copy_n (rhs.data (), rnz, data () + li);
std::copy_n (rhs.ridx (), rnz, ridx () + li);
mx_inline_add (ub - lb, cidx () + lb + 1, rhs.cidx () + 1,
li);
// ...tail.
std::copy (tmp.data () + ui, tmp.data () + nz,
data () + li + rnz);
std::copy (tmp.ridx () + ui, tmp.ridx () + nz,
ridx () + li + rnz);
mx_inline_add (nc - ub, cidx () + ub + 1,
tmp.cidx () + ub + 1, new_nz - nz);
assert (nnz () == new_nz);
}
}
else if (idx_j.is_range () && idx_j.increment () == -1)
{
// It's s(:,u:-1:l) = r. Reverse the assignment.
assign (idx_i, idx_j.sorted (),
rhs.index (idx_i, octave::idx_vector (m - 1, 0, -1)));
}
else if (idx_j.is_permutation (nc))
{
*this = rhs.index (idx_i, idx_j.inverse_permutation (nc));
}
else
{
const Sparse tmp = *this;
*this = Sparse (nr, nc);
OCTAVE_LOCAL_BUFFER_INIT (octave_idx_type, jsav, nc, -1);
// Assemble column lengths.
for (octave_idx_type i = 0; i < nc; i++)
xcidx (i+1) = tmp.cidx (i+1) - tmp.cidx (i);
for (octave_idx_type i = 0; i < m; i++)
{
octave_idx_type j =idx_j(i);
jsav[j] = i;
xcidx (j+1) = rhs.cidx (i+1) - rhs.cidx (i);
}
// Make cumulative.
for (octave_idx_type i = 0; i < nc; i++)
xcidx (i+1) += xcidx (i);
change_capacity (nnz ());
// Merge columns.
for (octave_idx_type i = 0; i < nc; i++)
{
octave_idx_type l = xcidx (i);
octave_idx_type u = xcidx (i+1);
octave_idx_type j = jsav[i];
if (j >= 0)
{
// from rhs
octave_idx_type k = rhs.cidx (j);
std::copy_n (rhs.data () + k, u - l, xdata () + l);
std::copy_n (rhs.ridx () + k, u - l, xridx () + l);
}
else
{
// original
octave_idx_type k = tmp.cidx (i);
std::copy_n (tmp.data () + k, u - l, xdata () + l);
std::copy_n (tmp.ridx () + k, u - l, xridx () + l);
}
}
}
}
else if (nc == 1 && idx_j.is_colon_equiv (nc) && idx_i.isvector ())
{
// It's just vector indexing. The 1D assign is specialized for that.
assign (idx_i, rhs);
}
else if (idx_j.is_colon ())
{
if (idx_i.is_permutation (nr))
{
*this = rhs.index (idx_i.inverse_permutation (nr), idx_j);
}
else
{
// FIXME: optimize more special cases?
// In general this requires unpacking the columns, which is slow,
// especially for many small columns. OTOH, transpose is an
// efficient O(nr+nc+nnz) operation.
*this = transpose ();
assign (octave::idx_vector::colon, idx_i, rhs.transpose ());
*this = transpose ();
}
}
else
{
// Split it into 2 assignments and one indexing.
Sparse tmp = index (octave::idx_vector::colon, idx_j);
tmp.assign (idx_i, octave::idx_vector::colon, rhs);
assign (octave::idx_vector::colon, idx_j, tmp);
}
}
else if (m == 1 && n == 1)
{
n = idx_i.length (nr);
m = idx_j.length (nc);
if (rhs.nnz () != 0)
assign (idx_i, idx_j, Sparse (n, m, rhs.data (0)));
else
assign (idx_i, idx_j, Sparse (n, m));
}
else if (idx_i.length (nr) == m && idx_j.length (nc) == n
&& (n == 1 || m == 1))
{
assign (idx_i, idx_j, rhs.transpose ());
}
else
octave::err_nonconformant ("=", idx_i.length (nr), idx_j.length (nc), n, m);
}
template
OCTAVE_API
void
Sparse::assign (const octave::idx_vector& idx_i,
const octave::idx_vector& idx_j, const T& rhs)
{
// FIXME: Converting the RHS and forwarding to the sparse matrix
// assignment function is simpler, but it might be good to have a
// specialization...
assign (idx_i, idx_j, Sparse (1, 1, rhs));
}
// Can't use versions of these in Array.cc due to duplication of the
// instantiations for Array, etc.
template
OCTAVE_API
bool
sparse_ascending_compare (typename ref_param::type a,
typename ref_param::type b)
{
return (a < b);
}
template
OCTAVE_API
bool
sparse_descending_compare (typename ref_param::type a,
typename ref_param::type b)
{
return (a > b);
}
template
OCTAVE_API
Sparse
Sparse::sort (octave_idx_type dim, sortmode mode) const
{
Sparse m = *this;
octave_idx_type nr = m.rows ();
octave_idx_type nc = m.columns ();
if (m.numel () < 1 || dim > 1)
return m;
bool sort_by_column = (dim > 0);
if (sort_by_column)
{
m = m.transpose ();
std::swap (nr, nc);
}
octave_sort lsort;
if (mode == ASCENDING)
lsort.set_compare (sparse_ascending_compare);
else if (mode == DESCENDING)
lsort.set_compare (sparse_descending_compare);
else
(*current_liboctave_error_handler)
("Sparse::sort: invalid MODE");
T *v = m.data ();
octave_idx_type *mcidx = m.cidx ();
octave_idx_type *mridx = m.ridx ();
for (octave_idx_type j = 0; j < nc; j++)
{
octave_idx_type ns = mcidx[j + 1] - mcidx[j];
lsort.sort (v, ns);
octave_idx_type i;
if (mode == ASCENDING)
{
for (i = 0; i < ns; i++)
if (sparse_ascending_compare (static_cast (0), v[i]))
break;
}
else
{
for (i = 0; i < ns; i++)
if (sparse_descending_compare (static_cast (0), v[i]))
break;
}
for (octave_idx_type k = 0; k < i; k++)
mridx[k] = k;
for (octave_idx_type k = i; k < ns; k++)
mridx[k] = k - ns + nr;
v += ns;
mridx += ns;
}
if (sort_by_column)
m = m.transpose ();
return m;
}
template
OCTAVE_API
Sparse
Sparse::sort (Array& sidx, octave_idx_type dim,
sortmode mode) const
{
Sparse m = *this;
octave_idx_type nr = m.rows ();
octave_idx_type nc = m.columns ();
if (m.numel () < 1 || dim > 1)
{
sidx = Array (dim_vector (nr, nc), 1);
return m;
}
bool sort_by_column = (dim > 0);
if (sort_by_column)
{
m = m.transpose ();
std::swap (nr, nc);
}
octave_sort indexed_sort;
if (mode == ASCENDING)
indexed_sort.set_compare (sparse_ascending_compare);
else if (mode == DESCENDING)
indexed_sort.set_compare (sparse_descending_compare);
else
(*current_liboctave_error_handler)
("Sparse::sort: invalid MODE");
T *v = m.data ();
octave_idx_type *mcidx = m.cidx ();
octave_idx_type *mridx = m.ridx ();
sidx = Array (dim_vector (nr, nc));
OCTAVE_LOCAL_BUFFER (octave_idx_type, vi, nr);
for (octave_idx_type j = 0; j < nc; j++)
{
octave_idx_type ns = mcidx[j + 1] - mcidx[j];
octave_idx_type offset = j * nr;
if (ns == 0)
{
for (octave_idx_type k = 0; k < nr; k++)
sidx(offset + k) = k;
}
else
{
for (octave_idx_type i = 0; i < ns; i++)
vi[i] = mridx[i];
indexed_sort.sort (v, vi, ns);
octave_idx_type i;
if (mode == ASCENDING)
{
for (i = 0; i < ns; i++)
if (sparse_ascending_compare (static_cast (0), v[i]))
break;
}
else
{
for (i = 0; i < ns; i++)
if (sparse_descending_compare (static_cast (0), v[i]))
break;
}
octave_idx_type ii = 0;
octave_idx_type jj = i;
for (octave_idx_type k = 0; k < nr; k++)
{
if (ii < ns && mridx[ii] == k)
ii++;
else
sidx(offset + jj++) = k;
}
for (octave_idx_type k = 0; k < i; k++)
{
sidx(k + offset) = vi[k];
mridx[k] = k;
}
for (octave_idx_type k = i; k < ns; k++)
{
sidx(k - ns + nr + offset) = vi[k];
mridx[k] = k - ns + nr;
}
v += ns;
mridx += ns;
}
}
if (sort_by_column)
{
m = m.transpose ();
sidx = sidx.transpose ();
}
return m;
}
template
OCTAVE_API
Sparse
Sparse::diag (octave_idx_type k) const
{
octave_idx_type nnr = rows ();
octave_idx_type nnc = cols ();
Sparse d;
if (nnr == 0 || nnc == 0)
; // do nothing
else if (nnr != 1 && nnc != 1)
{
if (k > 0)
nnc -= k;
else if (k < 0)
nnr += k;
if (nnr > 0 && nnc > 0)
{
octave_idx_type ndiag = (nnr < nnc) ? nnr : nnc;
// Count the number of nonzero elements
octave_idx_type nel = 0;
if (k > 0)
{
for (octave_idx_type i = 0; i < ndiag; i++)
if (elem (i, i+k) != 0.)
nel++;
}
else if (k < 0)
{
for (octave_idx_type i = 0; i < ndiag; i++)
if (elem (i-k, i) != 0.)
nel++;
}
else
{
for (octave_idx_type i = 0; i < ndiag; i++)
if (elem (i, i) != 0.)
nel++;
}
d = Sparse (ndiag, 1, nel);
d.xcidx (0) = 0;
d.xcidx (1) = nel;
octave_idx_type ii = 0;
if (k > 0)
{
for (octave_idx_type i = 0; i < ndiag; i++)
{
T tmp = elem (i, i+k);
if (tmp != 0.)
{
d.xdata (ii) = tmp;
d.xridx (ii++) = i;
}
}
}
else if (k < 0)
{
for (octave_idx_type i = 0; i < ndiag; i++)
{
T tmp = elem (i-k, i);
if (tmp != 0.)
{
d.xdata (ii) = tmp;
d.xridx (ii++) = i;
}
}
}
else
{
for (octave_idx_type i = 0; i < ndiag; i++)
{
T tmp = elem (i, i);
if (tmp != 0.)
{
d.xdata (ii) = tmp;
d.xridx (ii++) = i;
}
}
}
}
else
{
// Matlab returns [] 0x1 for out-of-range diagonal
octave_idx_type nr = 0;
octave_idx_type nc = 1;
octave_idx_type nz = 0;
d = Sparse (nr, nc, nz);
}
}
else // one of dimensions == 1 (vector)
{
octave_idx_type roff = 0;
octave_idx_type coff = 0;
if (k > 0)
{
roff = 0;
coff = k;
}
else if (k < 0)
{
roff = -k;
coff = 0;
}
if (nnr == 1)
{
octave_idx_type n = nnc + std::abs (k);
octave_idx_type nz = nnz ();
d = Sparse (n, n, nz);
if (nnz () > 0)
{
for (octave_idx_type i = 0; i < coff+1; i++)
d.xcidx (i) = 0;
for (octave_idx_type j = 0; j < nnc; j++)
{
for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
{
d.xdata (i) = data (i);
d.xridx (i) = j + roff;
}
d.xcidx (j + coff + 1) = cidx (j+1);
}
for (octave_idx_type i = nnc + coff + 1; i < n + 1; i++)
d.xcidx (i) = nz;
}
}
else
{
octave_idx_type n = nnr + std::abs (k);
octave_idx_type nz = nnz ();
d = Sparse (n, n, nz);
if (nnz () > 0)
{
octave_idx_type ii = 0;
octave_idx_type ir = ridx (0);
for (octave_idx_type i = 0; i < coff+1; i++)
d.xcidx (i) = 0;
for (octave_idx_type i = 0; i < nnr; i++)
{
if (ir == i)
{
d.xdata (ii) = data (ii);
d.xridx (ii++) = ir + roff;
if (ii != nz)
ir = ridx (ii);
}
d.xcidx (i + coff + 1) = ii;
}
for (octave_idx_type i = nnr + coff + 1; i < n+1; i++)
d.xcidx (i) = nz;
}
}
}
return d;
}
template
OCTAVE_API
Sparse
Sparse::cat (int dim, octave_idx_type n, const Sparse *sparse_list)
{
// Default concatenation.
bool (dim_vector::*concat_rule) (const dim_vector&, int) = &dim_vector::concat;
if (dim == -1 || dim == -2)
{
concat_rule = &dim_vector::hvcat;
dim = -dim - 1;
}
else if (dim < 0)
(*current_liboctave_error_handler) ("cat: invalid dimension");
dim_vector dv;
octave_idx_type total_nz = 0;
if (dim != 0 && dim != 1)
(*current_liboctave_error_handler)
("cat: invalid dimension for sparse concatenation");
if (n == 1)
return sparse_list[0];
for (octave_idx_type i = 0; i < n; i++)
{
if (! (dv.*concat_rule) (sparse_list[i].dims (), dim))
(*current_liboctave_error_handler) ("cat: dimension mismatch");
total_nz += sparse_list[i].nnz ();
}
Sparse retval (dv, total_nz);
if (retval.isempty ())
return retval;
switch (dim)
{
case 0:
{
// sparse vertcat. This is not efficiently handled by assignment,
// so we'll do it directly.
octave_idx_type l = 0;
for (octave_idx_type j = 0; j < dv(1); j++)
{
octave_quit ();
octave_idx_type rcum = 0;
for (octave_idx_type i = 0; i < n; i++)
{
const Sparse& spi = sparse_list[i];
// Skipping empty matrices. See the comment in Array.cc.
if (spi.isempty ())
continue;
octave_idx_type kl = spi.cidx (j);
octave_idx_type ku = spi.cidx (j+1);
for (octave_idx_type k = kl; k < ku; k++, l++)
{
retval.xridx (l) = spi.ridx (k) + rcum;
retval.xdata (l) = spi.data (k);
}
rcum += spi.rows ();
}
retval.xcidx (j+1) = l;
}
break;
}
case 1:
{
octave_idx_type l = 0;
for (octave_idx_type i = 0; i < n; i++)
{
octave_quit ();
// Skipping empty matrices. See the comment in Array.cc.
if (sparse_list[i].isempty ())
continue;
octave_idx_type u = l + sparse_list[i].columns ();
retval.assign (octave::idx_vector::colon, octave::idx_vector (l, u),
sparse_list[i]);
l = u;
}
break;
}
default:
assert (false);
}
return retval;
}
template
OCTAVE_API
Array
Sparse::array_value () const
{
Array retval (dims (), T ());
if (rows () == 1)
{
octave_idx_type i = 0;
for (octave_idx_type j = 0, nc = cols (); j < nc; j++)
{
if (cidx (j+1) > i)
retval.xelem (j) = data (i++);
}
}
else
{
for (octave_idx_type j = 0, nc = cols (); j < nc; j++)
for (octave_idx_type i = cidx (j), iu = cidx (j+1); i < iu; i++)
retval.xelem (ridx (i), j) = data (i);
}
return retval;
}
template
OCTAVE_API
std::istream&
read_sparse_matrix (std::istream& is, Sparse& a,
T (*read_fcn) (std::istream&))
{
octave_idx_type nr = a.rows ();
octave_idx_type nc = a.cols ();
octave_idx_type nz = a.nzmax ();
if (nr > 0 && nc > 0)
{
octave_idx_type itmp;
octave_idx_type jtmp;
octave_idx_type iold = 0;
octave_idx_type jold = 0;
octave_idx_type ii = 0;
T tmp;
a.cidx (0) = 0;
for (octave_idx_type i = 0; i < nz; i++)
{
itmp = 0; jtmp = 0;
is >> itmp;
itmp--;
is >> jtmp;
jtmp--;
if (is.fail ())
{
is.clear();
std::string err_field;
is >> err_field;
(*current_liboctave_error_handler)
("invalid sparse matrix: element %" OCTAVE_IDX_TYPE_FORMAT ": "
"Symbols '%s' is not an integer format",
i+1, err_field.c_str ());
}
if (itmp < 0 || itmp >= nr)
{
is.setstate (std::ios::failbit);
(*current_liboctave_error_handler)
("invalid sparse matrix: element %" OCTAVE_IDX_TYPE_FORMAT ": "
"row index = %" OCTAVE_IDX_TYPE_FORMAT " out of range",
i+1, itmp + 1);
}
if (jtmp < 0 || jtmp >= nc)
{
is.setstate (std::ios::failbit);
(*current_liboctave_error_handler)
("invalid sparse matrix: element %" OCTAVE_IDX_TYPE_FORMAT ": "
"column index = %" OCTAVE_IDX_TYPE_FORMAT " out of range",
i+1, jtmp + 1);
}
if (jtmp < jold)
{
is.setstate (std::ios::failbit);
(*current_liboctave_error_handler)
("invalid sparse matrix: element %" OCTAVE_IDX_TYPE_FORMAT ":"
"column indices must appear in ascending order "
"(%" OCTAVE_IDX_TYPE_FORMAT " < %" OCTAVE_IDX_TYPE_FORMAT ")",
i+1, jtmp, jold);
}
else if (jtmp > jold)
{
for (octave_idx_type j = jold; j < jtmp; j++)
a.cidx (j+1) = ii;
}
else if (itmp < iold)
{
is.setstate (std::ios::failbit);
(*current_liboctave_error_handler)
("invalid sparse matrix: element %" OCTAVE_IDX_TYPE_FORMAT ": "
"row indices must appear in ascending order in each column "
"(%" OCTAVE_IDX_TYPE_FORMAT " < %" OCTAVE_IDX_TYPE_FORMAT ")",
i+1, iold, itmp);
}
iold = itmp;
jold = jtmp;
tmp = read_fcn (is);
if (! is)
return is; // Problem, return is in error state
a.data (ii) = tmp;
a.ridx (ii++) = itmp;
}
for (octave_idx_type j = jold; j < nc; j++)
a.cidx (j+1) = ii;
}
return is;
}
/*
* Tests
*
%!function x = set_slice (x, dim, slice, arg)
%! switch (dim)
%! case 11
%! x(slice) = 2;
%! case 21
%! x(slice, :) = 2;
%! case 22
%! x(:, slice) = 2;
%! otherwise
%! error ("invalid dim, '%d'", dim);
%! endswitch
%!endfunction
%!function x = set_slice2 (x, dim, slice)
%! switch (dim)
%! case 11
%! x(slice) = 2 * ones (size (slice));
%! case 21
%! x(slice, :) = 2 * ones (length (slice), columns (x));
%! case 22
%! x(:, slice) = 2 * ones (rows (x), length (slice));
%! otherwise
%! error ("invalid dim, '%d'", dim);
%! endswitch
%!endfunction
%!function test_sparse_slice (size, dim, slice)
%! x = ones (size);
%! s = set_slice (sparse (x), dim, slice);
%! f = set_slice (x, dim, slice);
%! assert (nnz (s), nnz (f));
%! assert (full (s), f);
%! s = set_slice2 (sparse (x), dim, slice);
%! f = set_slice2 (x, dim, slice);
%! assert (nnz (s), nnz (f));
%! assert (full (s), f);
%!endfunction
#### 1d indexing
## size = [2 0]
%!test test_sparse_slice ([2 0], 11, []);
%!assert (set_slice (sparse (ones ([2 0])), 11, 1), sparse ([2 0]')) # sparse different from full
%!assert (set_slice (sparse (ones ([2 0])), 11, 2), sparse ([0 2]')) # sparse different from full
%!assert (set_slice (sparse (ones ([2 0])), 11, 3), sparse ([0 0; 2 0]')) # sparse different from full
%!assert (set_slice (sparse (ones ([2 0])), 11, 4), sparse ([0 0; 0 2]')) # sparse different from full
## size = [0 2]
%!test test_sparse_slice ([0 2], 11, []);
%!assert (set_slice (sparse (ones ([0 2])), 11, 1), sparse ([2 0])) # sparse different from full
%!test test_sparse_slice ([0 2], 11, 2);
%!test test_sparse_slice ([0 2], 11, 3);
%!test test_sparse_slice ([0 2], 11, 4);
%!test test_sparse_slice ([0 2], 11, [4, 4]);
## size = [2 1]
%!test test_sparse_slice ([2 1], 11, []);
%!test test_sparse_slice ([2 1], 11, 1);
%!test test_sparse_slice ([2 1], 11, 2);
%!test test_sparse_slice ([2 1], 11, 3);
%!test test_sparse_slice ([2 1], 11, 4);
%!test test_sparse_slice ([2 1], 11, [4, 4]);
## size = [1 2]
%!test test_sparse_slice ([1 2], 11, []);
%!test test_sparse_slice ([1 2], 11, 1);
%!test test_sparse_slice ([1 2], 11, 2);
%!test test_sparse_slice ([1 2], 11, 3);
%!test test_sparse_slice ([1 2], 11, 4);
%!test test_sparse_slice ([1 2], 11, [4, 4]);
## size = [2 2]
%!test test_sparse_slice ([2 2], 11, []);
%!test test_sparse_slice ([2 2], 11, 1);
%!test test_sparse_slice ([2 2], 11, 2);
%!test test_sparse_slice ([2 2], 11, 3);
%!test test_sparse_slice ([2 2], 11, 4);
%!test test_sparse_slice ([2 2], 11, [4, 4]);
# These 2 errors are the same as in the full case
%!error id=Octave:invalid-resize set_slice (sparse (ones ([2 2])), 11, 5)
%!error id=Octave:invalid-resize set_slice (sparse (ones ([2 2])), 11, 6)
#### 2d indexing
## size = [2 0]
%!test test_sparse_slice ([2 0], 21, []);
%!test test_sparse_slice ([2 0], 21, 1);
%!test test_sparse_slice ([2 0], 21, 2);
%!test test_sparse_slice ([2 0], 21, [2,2]);
%!assert (set_slice (sparse (ones ([2 0])), 21, 3), sparse (3,0))
%!assert (set_slice (sparse (ones ([2 0])), 21, 4), sparse (4,0))
%!test test_sparse_slice ([2 0], 22, []);
%!test test_sparse_slice ([2 0], 22, 1);
%!test test_sparse_slice ([2 0], 22, 2);
%!test test_sparse_slice ([2 0], 22, [2,2]);
%!assert (set_slice (sparse (ones ([2 0])), 22, 3), sparse ([0 0 2;0 0 2])) # sparse different from full
%!assert (set_slice (sparse (ones ([2 0])), 22, 4), sparse ([0 0 0 2;0 0 0 2])) # sparse different from full
## size = [0 2]
%!test test_sparse_slice ([0 2], 21, []);
%!test test_sparse_slice ([0 2], 21, 1);
%!test test_sparse_slice ([0 2], 21, 2);
%!test test_sparse_slice ([0 2], 21, [2,2]);
%!assert (set_slice (sparse (ones ([0 2])), 21, 3), sparse ([0 0;0 0;2 2])) # sparse different from full
%!assert (set_slice (sparse (ones ([0 2])), 21, 4), sparse ([0 0;0 0;0 0;2 2])) # sparse different from full
%!test test_sparse_slice ([0 2], 22, []);
%!test test_sparse_slice ([0 2], 22, 1);
%!test test_sparse_slice ([0 2], 22, 2);
%!test test_sparse_slice ([0 2], 22, [2,2]);
%!assert (set_slice (sparse (ones ([0 2])), 22, 3), sparse (0,3))
%!assert (set_slice (sparse (ones ([0 2])), 22, 4), sparse (0,4))
## size = [2 1]
%!test test_sparse_slice ([2 1], 21, []);
%!test test_sparse_slice ([2 1], 21, 1);
%!test test_sparse_slice ([2 1], 21, 2);
%!test test_sparse_slice ([2 1], 21, [2,2]);
%!test test_sparse_slice ([2 1], 21, 3);
%!test test_sparse_slice ([2 1], 21, 4);
%!test test_sparse_slice ([2 1], 22, []);
%!test test_sparse_slice ([2 1], 22, 1);
%!test test_sparse_slice ([2 1], 22, 2);
%!test test_sparse_slice ([2 1], 22, [2,2]);
%!test test_sparse_slice ([2 1], 22, 3);
%!test test_sparse_slice ([2 1], 22, 4);
## size = [1 2]
%!test test_sparse_slice ([1 2], 21, []);
%!test test_sparse_slice ([1 2], 21, 1);
%!test test_sparse_slice ([1 2], 21, 2);
%!test test_sparse_slice ([1 2], 21, [2,2]);
%!test test_sparse_slice ([1 2], 21, 3);
%!test test_sparse_slice ([1 2], 21, 4);
%!test test_sparse_slice ([1 2], 22, []);
%!test test_sparse_slice ([1 2], 22, 1);
%!test test_sparse_slice ([1 2], 22, 2);
%!test test_sparse_slice ([1 2], 22, [2,2]);
%!test test_sparse_slice ([1 2], 22, 3);
%!test test_sparse_slice ([1 2], 22, 4);
## size = [2 2]
%!test test_sparse_slice ([2 2], 21, []);
%!test test_sparse_slice ([2 2], 21, 1);
%!test test_sparse_slice ([2 2], 21, 2);
%!test test_sparse_slice ([2 2], 21, [2,2]);
%!test test_sparse_slice ([2 2], 21, 3);
%!test test_sparse_slice ([2 2], 21, 4);
%!test test_sparse_slice ([2 2], 22, []);
%!test test_sparse_slice ([2 2], 22, 1);
%!test test_sparse_slice ([2 2], 22, 2);
%!test test_sparse_slice ([2 2], 22, [2,2]);
%!test test_sparse_slice ([2 2], 22, 3);
%!test test_sparse_slice ([2 2], 22, 4);
%!assert <*35570> (speye (3,1)(3:-1:1), sparse ([0; 0; 1]))
## Test removing columns
%!test <*36656>
%! s = sparse (magic (5));
%! s(:,2:4) = [];
%! assert (s, sparse (magic (5)(:, [1,5])));
%!test
%! s = sparse ([], [], [], 1, 1);
%! s(1,:) = [];
%! assert (s, sparse ([], [], [], 0, 1));
## Test (bug #37321)
%!test <*37321> a=sparse (0,0); assert (all (a) == sparse ([1]));
%!test <*37321> a=sparse (0,1); assert (all (a) == sparse ([1]));
%!test <*37321> a=sparse (1,0); assert (all (a) == sparse ([1]));
%!test <*37321> a=sparse (1,0); assert (all (a,2) == sparse ([1]));
%!test <*37321> a=sparse (1,0); assert (size (all (a,1)), [1 0]);
%!test <*37321> a=sparse (1,1);
%! assert (all (a) == sparse ([0]));
%! assert (size (all (a)), [1 1]);
%!test <*37321> a=sparse (2,1);
%! assert (all (a) == sparse ([0]));
%! assert (size (all (a)), [1 1]);
%!test <*37321> a=sparse (1,2);
%! assert (all (a) == sparse ([0]));
%! assert (size (all (a)), [1 1]);
%!test <*37321> a=sparse (2,2); assert (isequal (all (a), sparse ([0 0])));
## Test assigning row to a column slice
%!test <45589>
%! a = sparse (magic (3));
%! b = a;
%! a(1,:) = 1:3;
%! b(1,:) = (1:3)';
%! assert (a, b);
*/
template
OCTAVE_API
void
Sparse::print_info (std::ostream& os, const std::string& prefix) const
{
os << prefix << "m_rep address: " << m_rep << "\n"
<< prefix << "m_rep->m_nzmax: " << m_rep->m_nzmax << "\n"
<< prefix << "m_rep->m_nrows: " << m_rep->m_nrows << "\n"
<< prefix << "m_rep->m_ncols: " << m_rep->m_ncols << "\n"
<< prefix << "m_rep->m_data: " << m_rep->m_data << "\n"
<< prefix << "m_rep->m_ridx: " << m_rep->m_ridx << "\n"
<< prefix << "m_rep->m_cidx: " << m_rep->m_cidx << "\n"
<< prefix << "m_rep->m_count: " << m_rep->m_count << "\n";
}
#if defined (__clang__)
# define INSTANTIATE_SPARSE(T) \
template class OCTAVE_API Sparse; \
template OCTAVE_API std::istream& \
read_sparse_matrix (std::istream& is, Sparse