List of all members | Classes | Public Types | Public Member Functions | Static Public Member Functions | Protected Member Functions
casadi::SparsityInternal Class Reference

#include <sparsity_internal.hpp>

Detailed Description

Definition at line 44 of file sparsity_internal.hpp.

Inheritance diagram for casadi::SparsityInternal:
Inheritance graph
[legend]
Collaboration diagram for casadi::SparsityInternal:
Collaboration graph
[legend]

Public Types

using weak_ref_type = WeakRefInternal
 

Public Member Functions

 SparsityInternal (casadi_int nrow, casadi_int ncol, const casadi_int *colind, const casadi_int *row)
 Construct a sparsity pattern from arrays. More...
 
 ~SparsityInternal () override
 Destructor. More...
 
const std::vector< casadi_int > & sp () const
 Get number of rows (see public class) More...
 
casadi_int size1 () const
 Get number of rows (see public class) More...
 
casadi_int size2 () const
 Get number of columns (see public class) More...
 
const casadi_int * colind () const
 Get column offsets (see public class) More...
 
const casadi_int * row () const
 Get row indices (see public class) More...
 
casadi_int nnz () const
 Number of structural non-zeros. More...
 
Sparsity get_diag (std::vector< casadi_int > &mapping) const
 Get the diagonal of the matrix/create a diagonal matrix. More...
 
bool has_diag () const
 has diagonal entries? More...
 
Sparsity drop_diag () const
 Drop diagonal entries. More...
 
casadi_int dfs (casadi_int j, casadi_int top, std::vector< casadi_int > &xi, std::vector< casadi_int > &pstack, const std::vector< casadi_int > &pinv, std::vector< bool > &marked) const
 Depth-first search. More...
 
casadi_int scc (std::vector< casadi_int > &p, std::vector< casadi_int > &r) const
 Find the strongly connected components of a square matrix. More...
 
std::vector< casadi_int > amd () const
 Approximate minimal degree preordering. More...
 
Sparsity T () const
 Transpose the matrix. More...
 
Sparsity transpose (std::vector< casadi_int > &mapping, bool invert_mapping=false) const
 Transpose the matrix and get the reordering of the non-zero entries,. More...
 
bool is_transpose (const SparsityInternal &y) const
 Check if the sparsity is the transpose of another. More...
 
bool is_reshape (const SparsityInternal &y) const
 Check if the sparsity is a reshape of another. More...
 
void bfs (casadi_int n, std::vector< casadi_int > &wi, std::vector< casadi_int > &wj, std::vector< casadi_int > &queue, const std::vector< casadi_int > &imatch, const std::vector< casadi_int > &jmatch, casadi_int mark) const
 Breadth-first search for coarse decomposition. More...
 
casadi_int btf (std::vector< casadi_int > &rowperm, std::vector< casadi_int > &colperm, std::vector< casadi_int > &rowblock, std::vector< casadi_int > &colblock, std::vector< casadi_int > &coarse_rowblock, std::vector< casadi_int > &coarse_colblock) const
 Compute the Dulmage-Mendelsohn decomposition. More...
 
const Btf & btf () const
 Get cached block triangular form. More...
 
void dmperm (std::vector< casadi_int > &rowperm, std::vector< casadi_int > &colperm, std::vector< casadi_int > &rowblock, std::vector< casadi_int > &colblock, std::vector< casadi_int > &coarse_rowblock, std::vector< casadi_int > &coarse_colblock) const
 Compute the Dulmage-Mendelsohn decomposition. More...
 
void maxtrans (std::vector< casadi_int > &imatch, std::vector< casadi_int > &jmatch, Sparsity &trans, casadi_int seed) const
 Compute the maximum transversal (maximum matching) More...
 
void augment (casadi_int k, std::vector< casadi_int > &jmatch, casadi_int *cheap, std::vector< casadi_int > &w, casadi_int *js, casadi_int *is, casadi_int *ps) const
 Find an augmenting path. More...
 
Sparsity permute (const std::vector< casadi_int > &pinv, const std::vector< casadi_int > &q, casadi_int values) const
 C = A(p, q) where p and q are permutations of 0..m-1 and 0..n-1. More...
 
void permute (const std::vector< casadi_int > &pinv, const std::vector< casadi_int > &q, casadi_int values, std::vector< casadi_int > &colind_C, std::vector< casadi_int > &row_C) const
 C = A(p, q) where p and q are permutations of 0..m-1 and 0..n-1. More...
 
Sparsity multiply (const Sparsity &B) const
 C = A*B. More...
 
casadi_int scatter (casadi_int j, std::vector< casadi_int > &w, casadi_int mark, casadi_int *Ci, casadi_int nz) const
 x = x + beta * A(:, j), where x is a dense vector and A(:, j) is sparse More...
 
std::vector< casadi_int > get_row () const
 Get row() as a vector. More...
 
std::vector< casadi_int > get_colind () const
 Get colind() as a vector. More...
 
std::vector< casadi_int > get_col () const
 Get the column for each nonzero. More...
 
Sparsity _resize (casadi_int nrow, casadi_int ncol) const
 Resize. More...
 
Sparsity _reshape (casadi_int nrow, casadi_int ncol) const
 Reshape a sparsity, order of nonzeros remains the same. More...
 
casadi_int numel () const
 Number of elements. More...
 
casadi_int nnz_lower (bool strictly=false) const
 Number of non-zeros in the lower triangular half. More...
 
casadi_int nnz_upper (bool strictly=false) const
 Number of non-zeros in the upper triangular half. More...
 
casadi_int nnz_diag () const
 Number of non-zeros on the diagonal. More...
 
casadi_int bw_upper () const
 Upper half-bandwidth. More...
 
casadi_int bw_lower () const
 Lower half-bandwidth. More...
 
std::pair< casadi_int, casadi_int > size () const
 Shape. More...
 
bool is_scalar (bool scalar_and_dense) const
 Is scalar? More...
 
bool is_empty (bool both=false) const
 Check if the sparsity is empty. More...
 
bool is_dense () const
 Is dense? More...
 
bool is_row () const
 Check if the pattern is a row vector (i.e. size1()==1) More...
 
bool is_column () const
 Check if the pattern is a column vector (i.e. size2()==1) More...
 
bool is_vector () const
 Check if the pattern is a row or column vector. More...
 
bool is_diag () const
 Is diagonal? More...
 
bool is_square () const
 Is square? More...
 
bool is_permutation () const
 Is this a permutation matrix? More...
 
bool is_selection (bool allow_empty=false) const
 Is this a selection matrix. More...
 
bool is_orthonormal (bool allow_empty=false) const
 Are the rows and columns of the pattern orthonormal ? More...
 
bool is_orthonormal_rows (bool allow_empty=false) const
 Are the rows of the pattern orthonormal ? More...
 
bool is_orthonormal_columns (bool allow_empty=false) const
 Are the columns of the pattern orthonormal ? More...
 
bool is_symmetric () const
 Is symmetric? More...
 
bool is_tril (bool strictly) const
 Is lower triangular? More...
 
bool is_triu (bool strictly) const
 is upper triangular? More...
 
Sparsity _triu (bool includeDiagonal) const
 Get upper triangular part. More...
 
Sparsity _tril (bool includeDiagonal) const
 Get lower triangular part. More...
 
std::vector< casadi_int > get_lower () const
 Get nonzeros in lower triangular part. More...
 
std::vector< casadi_int > get_upper () const
 Get nonzeros in upper triangular part. More...
 
std::string dim (bool with_nz=false) const
 Get the dimension as a string. More...
 
std::string repr_el (casadi_int k) const
 Describe the nonzero location k as a string. More...
 
Sparsity _mtimes (const Sparsity &y) const
 Sparsity pattern for a matrix-matrix product (details in public class) More...
 
bool is_subset (const Sparsity &rhs) const
 Is subset? More...
 
Sparsity pattern_inverse () const
 Take the inverse of a sparsity pattern; flip zeros and non-zeros. More...
 
bool is_equal (const Sparsity &y) const
 Check if two sparsity patterns are the same. More...
 
bool is_equal (casadi_int y_nrow, casadi_int y_ncol, const std::vector< casadi_int > &y_colind, const std::vector< casadi_int > &y_row) const
 Check if two sparsity patterns are the same. More...
 
bool is_equal (casadi_int y_nrow, casadi_int y_ncol, const casadi_int *y_colind, const casadi_int *y_row) const
 Check if two sparsity patterns are the same. More...
 
bool is_stacked (const Sparsity &y, casadi_int n) const
 Check if pattern is repeated. More...
 
Sparsity _enlargeRows (casadi_int nrow, const std::vector< casadi_int > &rr, bool ind1) const
 Enlarge the matrix along the first dimension (i.e. insert rows) More...
 
Sparsity _enlargeColumns (casadi_int ncol, const std::vector< casadi_int > &cc, bool ind1) const
 Enlarge the matrix along the second dimension (i.e. insert columns) More...
 
Sparsity makeDense (std::vector< casadi_int > &mapping) const
 Make a patten dense. More...
 
Sparsity _erase (const std::vector< casadi_int > &rr, const std::vector< casadi_int > &cc, bool ind1, std::vector< casadi_int > &mapping) const
 Erase rows and/or columns - does bounds checking. More...
 
Sparsity _erase (const std::vector< casadi_int > &rr, bool ind1, std::vector< casadi_int > &mapping) const
 Erase elements. More...
 
Sparsity _appendVector (const SparsityInternal &sp) const
 Append another sparsity patten vertically (vectors only) More...
 
Sparsity _appendColumns (const SparsityInternal &sp) const
 Append another sparsity patten horizontally. More...
 
Sparsity sub (const std::vector< casadi_int > &rr, const std::vector< casadi_int > &cc, std::vector< casadi_int > &mapping, bool ind1) const
 Get a submatrix. More...
 
Sparsity sub (const std::vector< casadi_int > &rr, const SparsityInternal &sp, std::vector< casadi_int > &mapping, bool ind1) const
 Get a set of elements. More...
 
casadi_int get_nz (casadi_int rr, casadi_int cc) const
 Get the index of an existing non-zero element. More...
 
std::vector< casadi_int > get_nz (const std::vector< casadi_int > &rr, const std::vector< casadi_int > &cc) const
 Get a set of non-zero element - does bounds checking. More...
 
void get_nz (std::vector< casadi_int > &indices) const
 Get the nonzero index for a set of elements (see description in public class) More...
 
bool rowsSequential (bool strictly) const
 Does the rows appear sequentially on each col. More...
 
Sparsity _removeDuplicates (std::vector< casadi_int > &mapping) const
 Remove duplicate entries. More...
 
void find (std::vector< casadi_int > &loc, bool ind1) const
 Get element index for each nonzero. More...
 
std::size_t hash () const
 Hash the sparsity pattern. More...
 
std::string class_name () const override
 Readable name of the internal class. More...
 
void disp (std::ostream &stream, bool more) const override
 Print description. More...
 
Sparsity uni_coloring (const Sparsity &AT, casadi_int cutoff) const
 Perform a unidirectional coloring. More...
 
Sparsity star_coloring (casadi_int ordering, casadi_int cutoff) const
 A greedy distance-2 coloring algorithm. More...
 
Sparsity star_coloring2 (casadi_int ordering, casadi_int cutoff) const
 An improved distance-2 coloring algorithm. More...
 
std::vector< casadi_int > largest_first () const
 Order the columns by decreasing degree. More...
 
Sparsity pmult (const std::vector< casadi_int > &p, bool permute_rows=true, bool permute_cols=true, bool invert_permutation=false) const
 Permute rows and/or columns. More...
 
void spy (std::ostream &stream) const
 Print a textual representation of sparsity. More...
 
void spy_matlab (const std::string &mfile) const
 Generate a script for Matlab or Octave which visualizes the sparsity using the spy command. More...
 
void export_code (const std::string &lang, std::ostream &stream, const Dict &options) const
 Export sparsity in Matlab format. More...
 
void spsolve (bvec_t *X, bvec_t *B, bool tr) const
 Propagate sparsity through a linear solve. More...
 
casadi_int getCount () const
 Get the reference count. More...
 
std::string debug_repr (const SharedObjectInternal *) const
 
GenericWeakRef< SharedObject, SharedObjectInternal > * weak ()
 Get a weak reference to the object. More...
 
Sparsity combine (const Sparsity &y, bool f0x_is_zero, bool function0_is_zero, std::vector< unsigned char > &mapping) const
 
Sparsity combine (const Sparsity &y, bool f0x_is_zero, bool function0_is_zero) const
 
template<bool with_mapping>
Sparsity combineGen1 (const Sparsity &y, bool f0x_is_zero, bool function0_is_zero, std::vector< unsigned char > &mapping) const
 
template<bool with_mapping, bool f0x_is_zero, bool function0_is_zero>
Sparsity combineGen (const Sparsity &y, std::vector< unsigned char > &mapping) const
 

Static Public Member Functions

static void etree (const casadi_int *sp, casadi_int *parent, casadi_int *w, casadi_int ata)
 Calculate the elimination tree for a matrix. More...
 
static casadi_int postorder_dfs (casadi_int j, casadi_int k, casadi_int *head, const casadi_int *next, casadi_int *post, casadi_int *stack)
 Traverse an elimination tree using depth first search. More...
 
static void postorder (const casadi_int *parent, casadi_int n, casadi_int *post, casadi_int *w)
 Calculate the postorder permuation. More...
 
static casadi_int leaf (casadi_int i, casadi_int j, const casadi_int *first, casadi_int *maxfirst, casadi_int *prevleaf, casadi_int *ancestor, casadi_int *jleaf)
 Needed by casadi_qr_colind. More...
 
static casadi_int qr_counts (const casadi_int *tr_sp, const casadi_int *parent, const casadi_int *post, casadi_int *counts, casadi_int *w)
 Calculate the column offsets for the QR R matrix. More...
 
static casadi_int qr_nnz (const casadi_int *sp, casadi_int *pinv, casadi_int *leftmost, const casadi_int *parent, casadi_int *nrow_ext, casadi_int *w)
 Calculate the number of nonzeros in the QR V matrix. More...
 
static void qr_init (const casadi_int *sp, const casadi_int *sp_tr, casadi_int *leftmost, casadi_int *parent, casadi_int *pinv, casadi_int *nrow_ext, casadi_int *v_nnz, casadi_int *r_nnz, casadi_int *w)
 Setup QP solver. More...
 
static void qr_sparsities (const casadi_int *sp_a, casadi_int nrow_ext, casadi_int *sp_v, casadi_int *sp_r, const casadi_int *leftmost, const casadi_int *parent, const casadi_int *pinv, casadi_int *iw)
 Get the row indices for V and R in QR factorization. More...
 
static void ldl_colind (const casadi_int *sp, casadi_int *parent, casadi_int *l_colind, casadi_int *w)
 Calculate the column offsets for the L factor of an LDL^T factorization. More...
 
static void ldl_row (const casadi_int *sp, const casadi_int *parent, casadi_int *l_colind, casadi_int *l_row, casadi_int *w)
 Calculate the row indices for the L factor of an LDL^T factorization. More...
 
static void matched (casadi_int n, const std::vector< casadi_int > &wj, const std::vector< casadi_int > &imatch, std::vector< casadi_int > &p, std::vector< casadi_int > &q, std::vector< casadi_int > &cc, std::vector< casadi_int > &rr, casadi_int set, casadi_int mark)
 Collect matched columns and rows into p and q. More...
 
static void unmatched (casadi_int m, const std::vector< casadi_int > &wi, std::vector< casadi_int > &p, std::vector< casadi_int > &rr, casadi_int set)
 Collect unmatched columns into the permutation vector p. More...
 
static casadi_int rprune (casadi_int i, casadi_int j, double aij, void *other)
 return 1 if column i is in R2 More...
 
static casadi_int drop (casadi_int(*fkeep)(casadi_int, casadi_int, double, void *), void *other, casadi_int nrow, casadi_int ncol, std::vector< casadi_int > &colind, std::vector< casadi_int > &row)
 drop entries for which fkeep(A(i, j)) is false; return nz if OK, else -1 More...
 
static std::vector< casadi_int > randperm (casadi_int n, casadi_int seed)
 return a random permutation vector More...
 
static std::vector< casadi_int > invertPermutation (const std::vector< casadi_int > &p)
 Invert a permutation vector. More...
 
static casadi_int wclear (casadi_int mark, casadi_int lemax, casadi_int *w, casadi_int n)
 clear w More...
 
static casadi_int diag (casadi_int i, casadi_int j, double aij, void *other)
 keep off-diagonal entries; drop diagonal entries More...
 

Protected Member Functions

void initSingleton ()
 
void destroySingleton ()
 
B shared_from_this ()
 Get a shared object from the current internal object. More...
 
const B shared_from_this () const
 Get a shared object from the current internal object. More...
 

Member Typedef Documentation

◆ weak_ref_type

using casadi::SharedObjectInternal::weak_ref_type = WeakRefInternal
inherited

Definition at line 152 of file shared_object.hpp.

Constructor & Destructor Documentation

◆ SparsityInternal()

casadi::SparsityInternal::SparsityInternal ( casadi_int  nrow,
casadi_int  ncol,
const casadi_int *  colind,
const casadi_int *  row 
)

Definition at line 449 of file sparsity_internal.cpp.

451  :
452  sp_(2 + ncol+1 + colind[ncol]), btf_(nullptr) {
453  sp_[0] = nrow;
454  sp_[1] = ncol;
455  std::copy(colind, colind+ncol+1, sp_.begin()+2);
456  std::copy(row, row+colind[ncol], sp_.begin()+2+ncol+1);
457  }
casadi::SparsityInternal::row
const casadi_int * row() const
Get row indices (see public class)
Definition: sparsity_internal.hpp:110
casadi::SparsityInternal::colind
const casadi_int * colind() const
Get column offsets (see public class)
Definition: sparsity_internal.hpp:105

References colind(), and row().

◆ ~SparsityInternal()

casadi::SparsityInternal::~SparsityInternal ( )
override

Definition at line 459 of file sparsity_internal.cpp.

459  {
460  delete btf_;
461  }

Member Function Documentation

◆ _appendColumns()

Sparsity casadi::SparsityInternal::_appendColumns ( const SparsityInternal &  sp) const

Definition at line 2577 of file sparsity_internal.cpp.

2577  {
2578  casadi_assert(size1()== sp.size1(),
2579  "_appendColumns(sp): row sizes must match but got " + str(size1())
2580  + " for lhs, and " + str(sp.size1()) + " for rhs.");
2581 
2582  // Append rows
2583  std::vector<casadi_int> new_row = get_row();
2584  const casadi_int* sp_row = sp.row();
2585  new_row.insert(new_row.end(), sp_row, sp_row+sp.nnz());
2586 
2587  // Get column indices
2588  std::vector<casadi_int> new_colind = get_colind();
2589  const casadi_int* sp_colind = sp.colind();
2590  new_colind.resize(size2() + sp.size2() + 1);
2591  for (casadi_int i = size2()+1; i<new_colind.size(); ++i)
2592  new_colind[i] = sp_colind[i-size2()] + nnz();
2593 
2594  return Sparsity(size1(), size2()+sp.size2(), new_colind, new_row);
2595  }
casadi::SparsityInternal::sp
const std::vector< casadi_int > & sp() const
Get number of rows (see public class)
Definition: sparsity_internal.hpp:90
casadi::SparsityInternal::size1
casadi_int size1() const
Get number of rows (see public class)
Definition: sparsity_internal.hpp:95
casadi::SparsityInternal::size2
casadi_int size2() const
Get number of columns (see public class)
Definition: sparsity_internal.hpp:100
casadi::SparsityInternal::get_colind
std::vector< casadi_int > get_colind() const
Get colind() as a vector.
Definition: sparsity_internal.cpp:3835
casadi::SparsityInternal::get_row
std::vector< casadi_int > get_row() const
Get row() as a vector.
Definition: sparsity_internal.cpp:3840
casadi::SparsityInternal::nnz
casadi_int nnz() const
Number of structural non-zeros.
Definition: sparsity_internal.hpp:113
casadi::str
std::string str(const T &v)
String representation, any type.
Definition: casadi_common.hpp:247

References get_colind(), get_row(), nnz(), size1(), size2(), sp(), and casadi::str().

Referenced by casadi::Sparsity::appendColumns().

◆ _appendVector()

Sparsity casadi::SparsityInternal::_appendVector ( const SparsityInternal &  sp) const

Definition at line 2556 of file sparsity_internal.cpp.

2556  {
2557  casadi_assert(size2() == 1 && sp.size2() == 1,
2558  "_appendVector(sp): Both arguments must be vectors but got "
2559  + str(size2()) + " columns for lhs, and " + str(sp.size2()) + " columns for rhs.");
2560 
2561  // Get current number of non-zeros
2562  casadi_int sz = nnz();
2563 
2564  // Add row indices
2565  std::vector<casadi_int> new_row = get_row();
2566  const casadi_int* sp_row = sp.row();
2567  new_row.resize(sz + sp.nnz());
2568  for (casadi_int i=sz; i<new_row.size(); ++i)
2569  new_row[i] = sp_row[i-sz] + size1();
2570 
2571  // New column indices
2572  std::vector<casadi_int> new_colind(2, 0);
2573  new_colind[1] = new_row.size();
2574  return Sparsity(size1()+sp.size1(), 1, new_colind, new_row);
2575  }

References get_row(), nnz(), size1(), size2(), sp(), and casadi::str().

Referenced by casadi::Sparsity::append().

◆ _enlargeColumns()

Sparsity casadi::SparsityInternal::_enlargeColumns ( casadi_int  ncol,
const std::vector< casadi_int > &  cc,
bool  ind1 
) const

Definition at line 2597 of file sparsity_internal.cpp.

2598  {
2599  casadi_assert_in_range(cc, -ncol+ind1, ncol+ind1);
2600 
2601  // Handle index-1, negative indices
2602  if (ind1 || has_negative(cc)) {
2603  std::vector<casadi_int> cc_mod = cc;
2604  for (std::vector<casadi_int>::iterator i=cc_mod.begin(); i!=cc_mod.end(); ++i) {
2605  if (ind1) (*i)--;
2606  if (*i<0) *i += ncol;
2607  }
2608  return _enlargeColumns(ncol, cc_mod, false); // Call recursively
2609  }
2610 
2611  // Sparsify the columns
2612  std::vector<casadi_int> new_colind = get_colind();
2613  new_colind.resize(ncol+1, nnz());
2614 
2615  casadi_int ik=cc.back(); // need only to update from the last new index
2616  casadi_int nz=nnz(); // number of nonzeros up till this column
2617  for (casadi_int i=cc.size()-1; i>=0; --i) {
2618  // Update colindex for new columns
2619  for (; ik>cc[i]; --ik) {
2620  new_colind[ik] = nz;
2621  }
2622 
2623  // Update non-zero counter
2624  nz = new_colind[i];
2625 
2626  // Update colindex for old colums
2627  new_colind[cc[i]] = nz;
2628  }
2629 
2630  // Append zeros to the beginning
2631  for (; ik>=0; --ik) {
2632  new_colind[ik] = 0;
2633  }
2634  return Sparsity(size1(), ncol, new_colind, get_row());
2635  }
casadi::SparsityInternal::_enlargeColumns
Sparsity _enlargeColumns(casadi_int ncol, const std::vector< casadi_int > &cc, bool ind1) const
Enlarge the matrix along the second dimension (i.e. insert columns)
Definition: sparsity_internal.cpp:2597
casadi::has_negative
bool has_negative(const std::vector< T > &v)
Check if the vector has negative entries.
Definition: casadi_misc.hpp:776

References get_colind(), get_row(), casadi::has_negative(), nnz(), and size1().

◆ _enlargeRows()

Sparsity casadi::SparsityInternal::_enlargeRows ( casadi_int  nrow,
const std::vector< casadi_int > &  rr,
bool  ind1 
) const

Definition at line 2637 of file sparsity_internal.cpp.

2638  {
2639  casadi_assert_in_range(rr, -nrow+ind1, nrow+ind1);
2640 
2641  // Handle index-1, negative indices
2642  if (ind1 || has_negative(rr)) {
2643  std::vector<casadi_int> rr_mod = rr;
2644  for (std::vector<casadi_int>::iterator i=rr_mod.begin(); i!=rr_mod.end(); ++i) {
2645  if (ind1) (*i)--;
2646  if (*i<0) *i += nrow;
2647  }
2648  return _enlargeRows(nrow, rr_mod, false); // Call recursively
2649  }
2650 
2651  // Assert dimensions
2652  casadi_assert_dev(rr.size() == size1());
2653 
2654  // Begin by sparsify the rows
2655  std::vector<casadi_int> new_row = get_row();
2656  for (casadi_int k=0; k<nnz(); ++k) {
2657  new_row[k] = rr[new_row[k]];
2658  }
2659  return Sparsity(nrow, size2(), get_colind(), new_row);
2660  }
casadi::SparsityInternal::_enlargeRows
Sparsity _enlargeRows(casadi_int nrow, const std::vector< casadi_int > &rr, bool ind1) const
Enlarge the matrix along the first dimension (i.e. insert rows)
Definition: sparsity_internal.cpp:2637

References get_colind(), get_row(), casadi::has_negative(), nnz(), size1(), and size2().

◆ _erase() [1/2]

Sparsity casadi::SparsityInternal::_erase ( const std::vector< casadi_int > &  rr,
bool  ind1,
std::vector< casadi_int > &  mapping 
) const

Definition at line 1923 of file sparsity_internal.cpp.

1924  {
1925  // Quick return if nothing to erase
1926  if (rr.empty()) {
1927  mapping = range(nnz());
1928  return shared_from_this<Sparsity>();
1929  }
1930  casadi_assert_in_range(rr, -numel()+ind1, numel()+ind1);
1931 
1932  // Handle index-1, negative indices
1933  if (ind1 || has_negative(rr)) {
1934  std::vector<casadi_int> rr_mod = rr;
1935  for (std::vector<casadi_int>::iterator i=rr_mod.begin(); i!=rr_mod.end(); ++i) {
1936  if (ind1) (*i)--;
1937  if (*i<0) *i += numel();
1938  }
1939  return _erase(rr_mod, false, mapping); // Call recursively
1940  }
1941 
1942  // Sort rr in non-deceasing order, if needed
1943  if (!is_nondecreasing(rr)) {
1944  std::vector<casadi_int> rr_sorted = rr;
1945  std::sort(rr_sorted.begin(), rr_sorted.end());
1946  return _erase(rr_sorted, false, mapping);
1947  }
1948 
1949  // Mapping
1950  mapping.resize(0);
1951 
1952  // Quick return if no elements
1953  if (numel()==0) return shared_from_this<Sparsity>();
1954 
1955  // Reserve memory
1956  mapping.reserve(nnz());
1957 
1958  // Number of non-zeros
1959  casadi_int nz=0;
1960 
1961  // Elements to be erased
1962  std::vector<casadi_int>::const_iterator next_rr = rr.begin();
1963 
1964  // Return value
1965  std::vector<casadi_int> ret_colind = get_colind(), ret_row = get_row();
1966 
1967  // First and last index for the column (note colind_ is being overwritten)
1968  casadi_int k_first, k_last=0;
1969 
1970  // Loop over columns
1971  for (casadi_int j=0; j<size2(); ++j) {
1972  // Update k range
1973  k_first = k_last;
1974  k_last = ret_colind[j+1];
1975 
1976  // Loop over nonzeros
1977  for (casadi_int k=k_first; k<k_last; ++k) {
1978  // Get row
1979  casadi_int i=ret_row[k];
1980 
1981  // Corresponding element
1982  casadi_int el = i+j*size1();
1983 
1984  // Continue to the next element to skip
1985  while (next_rr!=rr.end() && *next_rr<el) next_rr++;
1986 
1987  // Skip element if necessary
1988  if (next_rr!=rr.end() && *next_rr==el) {
1989  next_rr++;
1990  continue;
1991  }
1992 
1993  // Keep element
1994  mapping.push_back(k);
1995 
1996  // Update row
1997  ret_row[nz++] = i;
1998  }
1999 
2000  // Update colind
2001  ret_colind[j+1] = nz;
2002  }
2003 
2004  // Truncate row vector
2005  ret_row.resize(nz);
2006 
2007  return Sparsity(size1(), size2(), ret_colind, ret_row);
2008  }
casadi::SparsityInternal::numel
casadi_int numel() const
Number of elements.
Definition: sparsity_internal.cpp:477
casadi::SparsityInternal::_erase
Sparsity _erase(const std::vector< casadi_int > &rr, const std::vector< casadi_int > &cc, bool ind1, std::vector< casadi_int > &mapping) const
Erase rows and/or columns - does bounds checking.
Definition: sparsity_internal.cpp:2010
casadi::range
std::vector< casadi_int > range(casadi_int start, casadi_int stop, casadi_int step, casadi_int len)
Range function.
Definition: casadi_misc.cpp:104
casadi::is_nondecreasing
bool is_nondecreasing(const std::vector< T > &v)
Check if the vector is non-decreasing.
Definition: casadi_misc.hpp:755

References _erase(), get_colind(), get_row(), casadi::has_negative(), casadi::is_nondecreasing(), nnz(), numel(), casadi::range(), size1(), and size2().

◆ _erase() [2/2]

Sparsity casadi::SparsityInternal::_erase ( const std::vector< casadi_int > &  rr,
const std::vector< casadi_int > &  cc,
bool  ind1,
std::vector< casadi_int > &  mapping 
) const

Definition at line 2010 of file sparsity_internal.cpp.

2012  {
2013  casadi_assert_in_range(rr, -size1()+ind1, size1()+ind1);
2014  casadi_assert_in_range(cc, -size2()+ind1, size2()+ind1);
2015 
2016  // Handle index-1, negative indices, non-monotone rr and cc
2017  if (ind1 || has_negative(rr) || has_negative(cc)
2018  || !is_nondecreasing(rr) || !is_nondecreasing(cc)) {
2019  // Create substitute rr
2020  std::vector<casadi_int> rr_mod = rr;
2021  for (std::vector<casadi_int>::iterator i=rr_mod.begin(); i!=rr_mod.end(); ++i) {
2022  if (ind1) (*i)--;
2023  if (*i<0) *i += size1();
2024  }
2025  std::sort(rr_mod.begin(), rr_mod.end());
2026 
2027  // Create substitute cc
2028  std::vector<casadi_int> cc_mod = cc;
2029  for (std::vector<casadi_int>::iterator i=cc_mod.begin(); i!=cc_mod.end(); ++i) {
2030  if (ind1) (*i)--;
2031  if (*i<0) *i += size2();
2032  }
2033  std::sort(cc_mod.begin(), cc_mod.end());
2034 
2035  // Call recursively
2036  return _erase(rr_mod, cc_mod, false, mapping);
2037  }
2038 
2039  // Mapping
2040  mapping.resize(0);
2041 
2042  // Quick return if no elements
2043  if (numel()==0) return shared_from_this<Sparsity>();
2044 
2045  // Reserve memory
2046  mapping.reserve(nnz());
2047 
2048  // Return value
2049  std::vector<casadi_int> ret_colind = get_colind(), ret_row = get_row();
2050 
2051  // Number of non-zeros
2052  casadi_int nz=0;
2053 
2054  // Columns to be erased
2055  std::vector<casadi_int>::const_iterator ie = cc.begin();
2056 
2057  // First and last index for the col
2058  casadi_int el_first=0, el_last=0;
2059 
2060  // Loop over columns
2061  for (casadi_int i=0; i<size2(); ++i) {
2062  // Update beginning and end of non-zero indices
2063  el_first = el_last;
2064  el_last = ret_colind[i+1];
2065 
2066  // Is it a col that can be deleted
2067  bool deletable_col = ie!=cc.end() && *ie==i;
2068  if (deletable_col) {
2069  ie++;
2070 
2071  // Rows to be erased
2072  std::vector<casadi_int>::const_iterator je = rr.begin();
2073 
2074  // Loop over nonzero elements of the col
2075  for (casadi_int el=el_first; el<el_last; ++el) {
2076  // Row
2077  casadi_int j=ret_row[el];
2078 
2079  // Continue to the next row to skip
2080  for (; je!=rr.end() && *je<j; ++je) {}
2081 
2082  // Remove row if necessary
2083  if (je!=rr.end() && *je==j) {
2084  je++;
2085  continue;
2086  }
2087 
2088  // Save old nonzero for each new nonzero
2089  mapping.push_back(el);
2090 
2091  // Update row and increase nonzero counter
2092  ret_row[nz++] = j;
2093  }
2094  } else {
2095  // Loop over nonzero elements of the col
2096  for (casadi_int el=el_first; el<el_last; ++el) {
2097  // Row
2098  casadi_int j=ret_row[el];
2099 
2100  // Save old nonzero for each new nonzero
2101  mapping.push_back(el);
2102 
2103  // Update row and increase nonzero counter
2104  ret_row[nz++] = j;
2105  }
2106  }
2107 
2108  // Register last nonzero of the col
2109  ret_colind[i+1]=nz;
2110  }
2111 
2112  // Truncate row matrix
2113  ret_row.resize(nz);
2114 
2115  return Sparsity(size1(), size2(), ret_colind, ret_row);
2116  }

References get_colind(), get_row(), casadi::has_negative(), casadi::is_nondecreasing(), nnz(), numel(), size1(), and size2().

Referenced by _erase().

◆ _mtimes()

Sparsity casadi::SparsityInternal::_mtimes ( const Sparsity &  y) const

Definition at line 1722 of file sparsity_internal.cpp.

1722  {
1723  // Dimensions of the result
1724  casadi_int d1 = size1();
1725  casadi_int d2 = y.size2();
1726 
1727  // Elementwise multiplication if one factor is scalar
1728  if (is_scalar(false)) {
1729  return is_dense() ? y : Sparsity(y.size());
1730  } else if (y.is_scalar(false)) {
1731  return y.is_dense() ? shared_from_this<Sparsity>() : Sparsity(size());
1732  }
1733 
1734  // Quick return if both are dense
1735  if (is_dense() && y.is_dense()) {
1736  return !is_empty() && !y.is_empty() ? Sparsity::dense(d1, d2) :
1737  Sparsity(d1, d2);
1738  }
1739 
1740  // Quick return if first factor is diagonal
1741  if (is_diag()) return y;
1742 
1743  // Quick return if second factor is diagonal
1744  if (y.is_diag()) return shared_from_this<Sparsity>();
1745 
1746  // Direct access to the vectors
1747  const casadi_int* x_row = row();
1748  const casadi_int* x_colind = colind();
1749  const casadi_int* y_row = y.row();
1750  const casadi_int* y_colind = y.colind();
1751 
1752  // Sparsity pattern of the result
1753  std::vector<casadi_int> row, col;
1754 
1755  // Temporary vector for avoiding duplicate nonzeros
1756  std::vector<casadi_int> tmp(d1, -1);
1757 
1758  // Loop over the nonzeros of y
1759  for (casadi_int cc=0; cc<d2; ++cc) {
1760  for (casadi_int kk=y_colind[cc]; kk<y_colind[cc+1]; ++kk) {
1761  casadi_int rr = y_row[kk];
1762 
1763  // Loop over corresponding columns of x
1764  for (casadi_int kk1=x_colind[rr]; kk1<x_colind[rr+1]; ++kk1) {
1765  casadi_int rr1 = x_row[kk1];
1766 
1767  // Add to pattern if not already encountered
1768  if (tmp[rr1]!=cc) {
1769  tmp[rr1] = cc;
1770  row.push_back(rr1);
1771  col.push_back(cc);
1772  }
1773  }
1774  }
1775  }
1776 
1777  // Assemble sparsity pattern and return
1778  return Sparsity::triplet(d1, d2, row, col);
1779  }
casadi::SparsityInternal::is_empty
bool is_empty(bool both=false) const
Check if the sparsity is empty.
Definition: sparsity_internal.cpp:1801
casadi::SparsityInternal::size
std::pair< casadi_int, casadi_int > size() const
Shape.
Definition: sparsity_internal.cpp:1919
casadi::SparsityInternal::is_diag
bool is_diag() const
Is diagonal?
Definition: sparsity_internal.cpp:1805
casadi::SparsityInternal::is_scalar
bool is_scalar(bool scalar_and_dense) const
Is scalar?
Definition: sparsity_internal.cpp:1781
casadi::SparsityInternal::is_dense
bool is_dense() const
Is dense?
Definition: sparsity_internal.cpp:1785
casadi::Sparsity::dense
static Sparsity dense(casadi_int nrow, casadi_int ncol=1)
Create a dense rectangular sparsity pattern *.
Definition: sparsity.cpp:1012
casadi::Sparsity::triplet
static Sparsity triplet(casadi_int nrow, casadi_int ncol, const std::vector< casadi_int > &row, const std::vector< casadi_int > &col, std::vector< casadi_int > &mapping, bool invert_mapping)
Create a sparsity pattern given the nonzeros in sparse triplet form *.
Definition: sparsity.cpp:1127

References casadi::Sparsity::colind(), colind(), casadi::Sparsity::dense(), casadi::Sparsity::is_dense(), is_dense(), casadi::Sparsity::is_diag(), is_diag(), casadi::Sparsity::is_empty(), is_empty(), is_scalar(), casadi::Sparsity::is_scalar(), casadi::Sparsity::row(), row(), casadi::Sparsity::size(), size(), size1(), casadi::Sparsity::size2(), and casadi::Sparsity::triplet().

Referenced by casadi::Sparsity::mtimes().

◆ _removeDuplicates()

Sparsity casadi::SparsityInternal::_removeDuplicates ( std::vector< casadi_int > &  mapping) const

The same indices will be removed from the mapping vector, which must have the same length as the number of nonzeros

Extra doc: https://github.com/casadi/casadi/wiki/L_fm

Definition at line 2789 of file sparsity_internal.cpp.

2789  {
2790  casadi_assert_dev(mapping.size()==nnz());
2791 
2792  // Return value (to be hashed)
2793  std::vector<casadi_int> ret_colind = get_colind(), ret_row = get_row();
2794 
2795  // Nonzero counter without duplicates
2796  casadi_int k_strict=0;
2797 
2798  // Loop over columns
2799  for (casadi_int i=0; i<size2(); ++i) {
2800 
2801  // Last row encountered on the col so far
2802  casadi_int lastrow = -1;
2803 
2804  // Save new col offset (cannot set it yet, since we will need the old value below)
2805  casadi_int new_colind = k_strict;
2806 
2807  // Loop over nonzeros (including duplicates)
2808  for (casadi_int k=ret_colind[i]; k<ret_colind[i+1]; ++k) {
2809 
2810  // Make sure that the rows appear sequentially
2811  casadi_assert(ret_row[k] >= lastrow, "rows are not sequential");
2812 
2813  // Skip if duplicate
2814  if (ret_row[k] == lastrow) continue;
2815 
2816  // update last row encounterd on the col
2817  lastrow = ret_row[k];
2818 
2819  // Update mapping
2820  mapping[k_strict] = mapping[k];
2821 
2822  // Update row index
2823  ret_row[k_strict] = ret_row[k];
2824 
2825  // Increase the strict nonzero counter
2826  k_strict++;
2827  }
2828 
2829  // Update col offset
2830  ret_colind[i] = new_colind;
2831  }
2832 
2833  // Finalize the sparsity pattern
2834  ret_colind[size2()] = k_strict;
2835  ret_row.resize(k_strict);
2836  mapping.resize(k_strict);
2837  return Sparsity(size1(), size2(), ret_colind, ret_row);
2838  }

References get_colind(), get_row(), nnz(), size1(), and size2().

◆ _reshape()

Sparsity casadi::SparsityInternal::_reshape ( casadi_int  nrow,
casadi_int  ncol 
) const

Definition at line 2706 of file sparsity_internal.cpp.

2706  {
2707  // If a dimension is negative, call recursively
2708  if (nrow<0 && ncol>0) {
2709  return _reshape(numel()/ncol, ncol);
2710  } else if (nrow>0 && ncol<0) {
2711  return _reshape(nrow, numel()/nrow);
2712  }
2713 
2714  casadi_assert(numel() == nrow*ncol,
2715  "reshape: number of elements must remain the same. Old shape is "
2716  + dim() + ". New shape is " + str(nrow) + "x" + str(ncol)
2717  + "=" + str(nrow*ncol) + ".");
2718  std::vector<casadi_int> ret_col(nnz());
2719  std::vector<casadi_int> ret_row(nnz());
2720  const casadi_int* colind = this->colind();
2721  const casadi_int* row = this->row();
2722  for (casadi_int i=0; i<size2(); ++i) {
2723  for (casadi_int el=colind[i]; el<colind[i+1]; ++el) {
2724  casadi_int j = row[el];
2725 
2726  // Element number
2727  casadi_int k_ret = j+i*size1();
2728 
2729  // Col and row in the new matrix
2730  casadi_int i_ret = k_ret/nrow;
2731  casadi_int j_ret = k_ret%nrow;
2732  ret_col[el] = i_ret;
2733  ret_row[el] = j_ret;
2734  }
2735  }
2736  return Sparsity::triplet(nrow, ncol, ret_row, ret_col);
2737  }
casadi::SparsityInternal::_reshape
Sparsity _reshape(casadi_int nrow, casadi_int ncol) const
Reshape a sparsity, order of nonzeros remains the same.
Definition: sparsity_internal.cpp:2706
casadi::SparsityInternal::dim
std::string dim(bool with_nz=false) const
Get the dimension as a string.
Definition: sparsity_internal.cpp:1703

References colind(), dim(), nnz(), numel(), row(), size1(), size2(), casadi::str(), and casadi::Sparsity::triplet().

Referenced by casadi::Sparsity::reshape().

◆ _resize()

Sparsity casadi::SparsityInternal::_resize ( casadi_int  nrow,
casadi_int  ncol 
) const

Definition at line 2739 of file sparsity_internal.cpp.

2739  {
2740  // Col and row index of the new
2741  std::vector<casadi_int> row_new, colind_new(ncol+1, 0);
2742  const casadi_int* colind = this->colind();
2743  const casadi_int* row = this->row();
2744 
2745  // Loop over the columns which may contain nonzeros
2746  casadi_int i;
2747  for (i=0; i<size2() && i<ncol; ++i) {
2748  // First nonzero element of the col
2749  colind_new[i] = row_new.size();
2750 
2751  // Record rows of the nonzeros
2752  for (casadi_int el=colind[i]; el<colind[i+1] && row[el]<nrow; ++el) {
2753  row_new.push_back(row[el]);
2754  }
2755  }
2756 
2757  // Save col-indices for the rest of the columns
2758  for (; i<ncol+1; ++i) {
2759  colind_new[i] = row_new.size();
2760  }
2761 
2762  return Sparsity(nrow, ncol, colind_new, row_new);
2763  }

References colind(), row(), and size2().

◆ _tril()

Sparsity casadi::SparsityInternal::_tril ( bool  includeDiagonal) const

Definition at line 3747 of file sparsity_internal.cpp.

3747  {
3748  const casadi_int* colind = this->colind();
3749  const casadi_int* row = this->row();
3750  std::vector<casadi_int> ret_colind, ret_row;
3751  ret_colind.reserve(size2()+1);
3752  ret_colind.push_back(0);
3753  for (casadi_int cc=0; cc<size2(); ++cc) {
3754  for (casadi_int el=colind[cc]; el<colind[cc+1]; ++el) {
3755  casadi_int rr=row[el];
3756  if (rr>cc || (includeDiagonal && rr==cc)) {
3757  ret_row.push_back(rr);
3758  }
3759  }
3760  ret_colind.push_back(ret_row.size());
3761  }
3762  return Sparsity(size1(), size2(), ret_colind, ret_row);
3763  }

References colind(), row(), size1(), and size2().

Referenced by casadi::Sparsity::tril().

◆ _triu()

Sparsity casadi::SparsityInternal::_triu ( bool  includeDiagonal) const

Definition at line 3765 of file sparsity_internal.cpp.

3765  {
3766  const casadi_int* colind = this->colind();
3767  const casadi_int* row = this->row();
3768  std::vector<casadi_int> ret_colind, ret_row;
3769  ret_colind.reserve(size2()+1);
3770  ret_colind.push_back(0);
3771  for (casadi_int cc=0; cc<size2(); ++cc) {
3772  for (casadi_int el=colind[cc]; el<colind[cc+1]; ++el) {
3773  casadi_int rr=row[el];
3774  if (rr<cc || (includeDiagonal && rr==cc)) {
3775  ret_row.push_back(rr);
3776  }
3777  }
3778  ret_colind.push_back(ret_row.size());
3779  }
3780  return Sparsity(size1(), size2(), ret_colind, ret_row);
3781  }

References colind(), row(), size1(), and size2().

Referenced by casadi::Sparsity::triu().

◆ amd()

std::vector< casadi_int > casadi::SparsityInternal::amd ( ) const

The implementation is a modified version of cs_amd in CSparse Copyright(c) Timothy A. Davis, 2006-2009 Licensed as a derivative work under the GNU LGPL

Extra doc: https://github.com/casadi/casadi/wiki/L_en

Definition at line 660 of file sparsity_internal.cpp.

660  {
661  /*
662  Modified version of cs_amd in CSparse
663  Copyright(c) Timothy A. Davis, 2006-2009
664  Licensed as a derivative work under the GNU LGPL
665  */
666  casadi_assert(is_symmetric(), "AMD requires a symmetric matrix");
667  // Get sparsity
668  casadi_int n=size2();
669  std::vector<casadi_int> colind = get_colind();
670  std::vector<casadi_int> row = get_row();
671  // Drop diagonal entries
672  casadi_int nnz = 0; // number of nonzeros after pruning
673  casadi_int col_begin, col_end=0;
674  for (casadi_int c=0; c<n; ++c) {
675  // Get the range of nonzeros for the column, before pruning
676  col_begin = col_end;
677  col_end = colind[c+1];
678  // Loop over nonzeros
679  for (casadi_int k=col_begin; k<col_end; ++k) {
680  if (row[k]!=c) {
681  row[nnz++] = row[k];
682  }
683  }
684  colind[c+1] = nnz;
685  }
686  // dense threshold
687  casadi_int dense = static_cast<casadi_int>(10*sqrt(static_cast<double>(n)));
688  dense = std::max(casadi_int(16), dense);
689  dense = std::min(n-2, dense);
690  // Allocate result
691  std::vector<casadi_int> P(n+1);
692  // Work vectors
693  std::vector<casadi_int> len(n+1), nv(n+1), next(n+1), head(n+1), elen(n+1), degree(n+1),
694  w(n+1), hhead(n+1);
695  // Number of elements
696  casadi_int nel = 0;
697  // Minimal degree
698  casadi_int mindeg = 0;
699  // Maximum length of w
700  casadi_int lemax = 0;
701  // Degree
702  casadi_int d;
703  // ?
704  casadi_uint h;
705  // Flip
706  #define FLIP(i) (-(i)-2)
707  // Elbow room
708  //casadi_int t = nnz + nnz/5 + 2*n;
709  // Initialize quotient graph
710  for (casadi_int k = 0; k<n; ++k) len[k] = colind[k+1] - colind[k];
711  len[n] = 0;
712  casadi_int nzmax = row.size();
713  for (casadi_int i=0; i<=n; ++i) {
714  head[i] = -1; // degree list i is empty
715  P[i] = -1;
716  next[i] = -1;
717  hhead[i] = -1; // hash list i is empty
718  nv[i] = 1; // node i is just one node
719  w[i] = 1; // node i is alive
720  elen[i] = 0; // Ek of node i is empty
721  degree[i] = len[i]; // degree of node i
722  }
723  casadi_int mark = wclear(0, 0, get_ptr(w), n); // clear w
724  elen[n] = -2; // n is a dead element
725  colind[n] = -1; // n is a root of assembly tree
726  w[n] = 0; // n is a dead element
727  // Initialize degree lists
728  for (casadi_int i = 0; i < n; ++i) {
729  d = degree[i];
730  if (d == 0) { // node i is empty
731  elen[i] = -2; // element i is dead
732  nel++;
733  colind[i] = -1; // i is a root of assembly tree
734  w[i] = 0;
735  } else if (d > dense) { // node i is dense
736  nv[i] = 0; // absorb i into element n
737  elen[i] = -1; // node i is dead
738  nel++;
739  colind[i] = FLIP(n);
740  nv[n]++;
741  } else {
742  if (head[d] != -1) P[head[d]] = i;
743  next[i] = head[d]; // put node i in degree list d
744  head[d] = i;
745  }
746  }
747  while (nel < n) { // while (selecting pivots) do
748  // Select node of minimum approximate degree
749  casadi_int k;
750  for (k = -1; mindeg < n && (k = head[mindeg]) == -1; mindeg++) {}
751  if (next[k] != -1) P[next[k]] = -1;
752  head[mindeg] = next[k]; // remove k from degree list
753  casadi_int elenk = elen[k]; // elenk = |Ek|
754  casadi_int nvk = nv[k]; // # of nodes k represents
755  nel += nvk; // nv[k] nodes of A eliminated
756  // Garbage collection
757  if (elenk > 0 && nnz + mindeg >= nzmax) {
758  for (casadi_int j = 0; j < n; j++) {
759  casadi_int p;
760  if ((p = colind[j]) >= 0) { // j is a live node or element
761  colind[j] = row[p]; // save first entry of object
762  row[p] = FLIP(j); // first entry is now FLIP(j)
763  }
764  }
765  casadi_int q, p;
766  for (q = 0, p = 0; p < nnz; ) { // scan all of memory
767  casadi_int j;
768  if ((j = FLIP(row[p++])) >= 0) { // found object j
769  row[q] = colind[j]; // restore first entry of object
770  colind[j] = q++; // new pointer to object j
771  for (casadi_int k3 = 0; k3 < len[j]-1; k3++) row[q++] = row[p++];
772  }
773  }
774  nnz = q; // row[nnz...nzmax-1] now free
775  }
776  // Construct new element
777  casadi_int dk = 0;
778  nv[k] = -nvk; // flag k as in Lk
779  casadi_int p = colind[k];
780  casadi_int pk1 = (elenk == 0) ? p : nnz; // do in place if elen[k] == 0
781  casadi_int pk2 = pk1;
782  casadi_int e, pj, ln;
783  for (casadi_int k1 = 1; k1 <= elenk + 1; k1++) {
784  if (k1 > elenk) {
785  e = k; // search the nodes in k
786  pj = p; // list of nodes starts at row[pj]
787  ln = len[k] - elenk; // length of list of nodes in k
788  } else {
789  e = row[p++]; // search the nodes in e
790  pj = colind[e];
791  ln = len[e]; // length of list of nodes in e
792  }
793  for (casadi_int k2 = 1; k2 <= ln; k2++) {
794  casadi_int i = row[pj++];
795  casadi_int nvi;
796  if ((nvi = nv[i]) <= 0) continue; // node i dead, or seen
797  dk += nvi; // degree[Lk] += size of node i
798  nv[i] = -nvi; // negate nv[i] to denote i in Lk
799  row[pk2++] = i; // place i in Lk
800  if (next[i] != -1) P[next[i]] = P[i];
801  if (P[i] != -1) { // remove i from degree list
802  next[P[i]] = next[i];
803  } else {
804  head[degree[i]] = next[i];
805  }
806  }
807  if (e != k) {
808  colind[e] = FLIP(k); // absorb e into k
809  w[e] = 0; // e is now a dead element
810  }
811  }
812  if (elenk != 0) nnz = pk2; // row[nnz...nzmax] is free
813  degree[k] = dk; // external degree of k - |Lk\i|
814  colind[k] = pk1; // element k is in row[pk1..pk2-1]
815  len[k] = pk2 - pk1;
816  elen[k] = -2; // k is now an element
817  // Find set differences
818  mark = wclear(mark, lemax, get_ptr(w), n); // clear w if necessary
819  for (casadi_int pk = pk1; pk < pk2; pk++) { // scan 1: find |Le\Lk|
820  casadi_int i = row[pk];
821  casadi_int eln;
822  if ((eln = elen[i]) <= 0) continue; // skip if elen[i] empty
823  casadi_int nvi = -nv[i]; // nv[i] was negated
824  casadi_int wnvi = mark - nvi;
825  for (p = colind[i]; p <= colind[i] + eln - 1; p++) { // scan Ei
826  e = row[p];
827  if (w[e] >= mark) {
828  w[e] -= nvi; // decrement |Le\Lk|
829  } else if (w[e] != 0) { // ensure e is a live element
830  w[e] = degree[e] + wnvi; /* 1st time e seen in scan 1 */
831  }
832  }
833  }
834  // Degree update
835  for (casadi_int pk = pk1; pk < pk2; pk++) { // scan2: degree update
836  casadi_int i = row[pk]; // consider node i in Lk
837  casadi_int p1 = colind[i];
838  casadi_int p2 = p1 + elen[i] - 1;
839  casadi_int pn = p1;
840  for (h = 0, d = 0, p = p1; p <= p2; p++) { // scan Ei
841  e = row[p];
842  if (w[e] != 0) { // e is an unabsorbed element
843  casadi_int dext = w[e] - mark; // dext = |Le\Lk|
844  if (dext > 0) {
845  d += dext; // sum up the set differences
846  row[pn++] = e; // keep e in Ei
847  h += e; // compute the hash of node i
848  } else {
849  colind[e] = FLIP(k); // aggressive absorb. e->k
850  w[e] = 0; // e is a dead element
851  }
852  }
853  }
854  elen[i] = pn - p1 + 1; // elen[i] = |Ei|
855  casadi_int p3 = pn;
856  casadi_int p4 = p1 + len[i];
857  for (p = p2 + 1; p < p4; p++) { // prune edges in Ai
858  casadi_int j = row[p];
859  casadi_int nvj;
860  if ((nvj = nv[j]) <= 0) continue; // node j dead or in Lk
861  d += nvj; // degree(i) += |j|
862  row[pn++] = j; // place j in node list of i
863  h += j; // compute hash for node i
864  }
865  if (d == 0) { // check for mass elimination
866  colind[i] = FLIP(k); // absorb i into k
867  casadi_int nvi = -nv[i];
868  dk -= nvi; // |Lk| -= |i|
869  nvk += nvi; // |k| += nv[i]
870  nel += nvi;
871  nv[i] = 0;
872  elen[i] = -1; // node i is dead
873  } else {
874  degree[i] = std::min(degree[i], d); // update degree(i)
875  row[pn] = row[p3]; // move first node to end
876  row[p3] = row[p1]; // move 1st el. to end of Ei
877  row[p1] = k; // add k as 1st element in of Ei
878  len[i] = pn - p1 + 1; // new len of adj. list of node i
879  h %= n; // finalize hash of i
880  next[i] = hhead[h]; // place i in hash bucket
881  hhead[h] = i;
882  P[i] = h; // save hash of i in P[i]
883  }
884  } // scan2 is done
885  degree[k] = dk; // finalize |Lk|
886  lemax = std::max(lemax, dk);
887  mark = wclear(mark+lemax, lemax, get_ptr(w), n); // clear w
888  // Supernode detection
889  for (casadi_int pk = pk1; pk < pk2; pk++) {
890  casadi_int i = row[pk];
891  if (nv[i] >= 0) continue; // skip if i is dead
892  h = P[i]; // scan hash bucket of node i
893  i = hhead[h];
894  hhead[h] = -1; // hash bucket will be empty
895  for (; i != -1 && next[i] != -1; i = next[i], mark++) {
896  ln = len[i];
897  casadi_int eln = elen[i];
898  for (p = colind[i]+1; p <= colind[i] + ln-1; p++) w[row[p]] = mark;
899  casadi_int jlast = i;
900  for (casadi_int j = next[i]; j != -1; ) { // compare i with all j
901  casadi_int ok = (len[j] == ln) && (elen[j] == eln);
902  for (p = colind[j] + 1; ok && p <= colind[j] + ln - 1; p++) {
903  if (w[row[p]] != mark) ok = 0; // compare i and j
904  }
905  if (ok) { // i and j are identical
906  colind[j] = FLIP(i); // absorb j into i
907  nv[i] += nv[j];
908  nv[j] = 0;
909  elen[j] = -1; // node j is dead
910  j = next[j]; // delete j from hash bucket
911  next[jlast] = j;
912  } else {
913  jlast = j; // j and i are different
914  j = next[j];
915  }
916  }
917  }
918  }
919  // Finalize new element
920  casadi_int pk;
921  for (p = pk1, pk = pk1; pk < pk2; pk++) { // finalize Lk
922  casadi_int i = row[pk];
923  casadi_int nvi;
924  if ((nvi = -nv[i]) <= 0) continue; // skip if i is dead
925  nv[i] = nvi; // restore nv[i]
926  d = degree[i] + dk - nvi; // compute external degree(i)
927  d = std::min(d, n - nel - nvi);
928  if (head[d] != -1) P[head[d]] = i;
929  next[i] = head[d]; // put i back in degree list
930  P[i] = -1;
931  head[d] = i;
932  mindeg = std::min(mindeg, d); // find new minimum degree
933  degree[i] = d;
934  row[p++] = i; // place i in Lk
935  }
936  nv[k] = nvk; // # nodes absorbed into k
937  if ((len[k] = p-pk1) == 0) { // length of adj list of element k
938  colind[k] = -1; // k is a root of the tree
939  w[k] = 0; // k is now a dead element
940  }
941  if (elenk != 0) nnz = p; // free unused space in Lk
942  }
943  // Postordering
944  for (casadi_int i = 0; i < n; i++) colind[i] = FLIP(colind[i]); // fix assembly tree
945  for (casadi_int j = 0; j <= n; j++) head[j] = -1;
946  for (casadi_int j = n; j >= 0; j--) { // place unordered nodes in lists
947  if (nv[j] > 0) continue; // skip if j is an element
948  next[j] = head[colind[j]]; // place j in list of its parent
949  head[colind[j]] = j;
950  }
951  for (casadi_int e = n; e >= 0; e--) { // place elements in lists
952  if (nv[e] <= 0) continue; // skip unless e is an element
953  if (colind[e] != -1) {
954  next[e] = head[colind[e]]; // place e in list of its parent
955  head[colind[e]] = e;
956  }
957  }
958  for (casadi_int k = 0, i = 0; i <= n; i++) { // postorder the assembly tree
959  if (colind[i] == -1) k = postorder_dfs(i, k, get_ptr(head), get_ptr(next),
960  get_ptr(P), get_ptr(w));
961  }
962  P.resize(n);
963  return P;
964  #undef FLIP
965  }
casadi::SparsityInternal::wclear
static casadi_int wclear(casadi_int mark, casadi_int lemax, casadi_int *w, casadi_int n)
clear w
Definition: sparsity_internal.cpp:1525
casadi::SparsityInternal::postorder_dfs
static casadi_int postorder_dfs(casadi_int j, casadi_int k, casadi_int *head, const casadi_int *next, casadi_int *post, casadi_int *stack)
Traverse an elimination tree using depth first search.
Definition: sparsity_internal.cpp:74
casadi::SparsityInternal::is_symmetric
bool is_symmetric() const
Is symmetric?
Definition: sparsity_internal.cpp:1835
casadi::get_ptr
T * get_ptr(std::vector< T > &v)
Get a pointer to the data contained in the vector.
Definition: casadi_misc.hpp:857

References colind(), get_colind(), casadi::get_ptr(), get_row(), is_symmetric(), nnz(), casadi::P, postorder_dfs(), row(), size2(), and wclear().

◆ augment()

void casadi::SparsityInternal::augment ( casadi_int  k,
std::vector< casadi_int > &  jmatch,
casadi_int *  cheap,
std::vector< casadi_int > &  w,
casadi_int *  js,
casadi_int *  is,
casadi_int *  ps 
) const

The implementation is a modified version of cs_augment in CSparse Copyright(c) Timothy A. Davis, 2006-2009 Licensed as a derivative work under the GNU LGPL

Extra doc: https://github.com/casadi/casadi/wiki/L_f6

Definition at line 1080 of file sparsity_internal.cpp.

1081  {
1082  /*
1083  Modified version of cs_augment in CSparse
1084  Copyright(c) Timothy A. Davis, 2006-2009
1085  Licensed as a derivative work under the GNU LGPL
1086  */
1087  const casadi_int* colind = this->colind();
1088  const casadi_int* row = this->row();
1089 
1090  casadi_int found = 0, p, i = -1, head = 0, j ;
1091 
1092  // start with just node k in jstack
1093  js[0] = k ;
1094 
1095  while (head >= 0) {
1096  // --- Start (or continue) depth-first-search at node j -------------
1097 
1098  // get j from top of jstack
1099  j = js[head];
1100 
1101  // 1st time j visited for kth path
1102  if (w[j] != k) {
1103 
1104  // mark j as visited for kth path
1105  w[j] = k;
1106  for (p = cheap[j] ; p < colind[j+1] && !found; ++p) {
1107  i = row[p] ; /* try a cheap assignment (i, j) */
1108  found = (jmatch[i] == -1) ;
1109  }
1110 
1111  // start here next time j is traversed
1112  cheap[j] = p;
1113  if (found) {
1114  // row j matched with col i
1115  is[head] = i;
1116 
1117  // end of augmenting path
1118  break;
1119  }
1120 
1121  // no cheap match: start dfs for j
1122  ps[head] = colind[j];
1123  }
1124 
1125  // --- Depth-first-search of neighbors of j -------------------------
1126  for (p = ps[head]; p<colind[j+1]; ++p) {
1127 
1128  // consider col i
1129  i = row[p];
1130 
1131  // skip jmatch[i] if marked
1132  if (w[jmatch[i]] == k) continue;
1133 
1134  // pause dfs of node j
1135  ps[head] = p + 1;
1136 
1137  // i will be matched with j if found
1138  is[head] = i;
1139 
1140  // start dfs at row jmatch[i]
1141  js[++head] = jmatch[i];
1142  break ;
1143  }
1144 
1145  // node j is done; pop from stack
1146  if (p == colind[j+1]) head--;
1147  } // augment the match if path found:
1148 
1149  if (found)
1150  for (p = head; p>=0; --p)
1151  jmatch[is[p]] = js[p];
1152  }

References colind(), and row().

◆ bfs()

void casadi::SparsityInternal::bfs ( casadi_int  n,
std::vector< casadi_int > &  wi,
std::vector< casadi_int > &  wj,
std::vector< casadi_int > &  queue,
const std::vector< casadi_int > &  imatch,
const std::vector< casadi_int > &  jmatch,
casadi_int  mark 
) const

The implementation is a modified version of cs_bfs in CSparse Copyright(c) Timothy A. Davis, 2006-2009 Licensed as a derivative work under the GNU LGPL

Extra doc: https://github.com/casadi/casadi/wiki/L_ez

Definition at line 967 of file sparsity_internal.cpp.

969  {
970  /*
971  Modified version of cs_bfs in CSparse
972  Copyright(c) Timothy A. Davis, 2006-2009
973  Licensed as a derivative work under the GNU LGPL
974  */
975  casadi_int head = 0, tail = 0, j, i, p, j2 ;
976 
977  // place all unmatched nodes in queue
978  for (j=0; j<n; ++j) {
979  // skip j if matched
980  if (imatch[j] >= 0) continue;
981 
982  // j in set C0 (R0 if transpose)
983  wj[j] = 0;
984 
985  // place unmatched row j in queue
986  queue[tail++] = j;
987  }
988 
989  // quick return if no unmatched nodes
990  if (tail == 0) return;
991 
992  Sparsity trans;
993  const casadi_int *C_row, *C_colind;
994  if (mark == 1) {
995  C_row = row();
996  C_colind = colind();
997  } else {
998  trans = T();
999  C_row = trans.row();
1000  C_colind = trans.colind();
1001  }
1002 
1003  // while queue is not empty
1004  while (head < tail) {
1005 
1006  // get the head of the queue
1007  j = queue[head++];
1008  for (p = C_colind[j] ; p < C_colind[j+1] ; p++) {
1009  i = C_row[p] ;
1010 
1011  // skip if i is marked
1012  if (wi[i] >= 0) continue;
1013 
1014  // i in set R1 (C3 if transpose)
1015  wi[i] = mark;
1016 
1017  // traverse alternating path to j2
1018  j2 = jmatch[i];
1019 
1020  // skip j2 if it is marked
1021  if (wj[j2] >= 0) continue;
1022 
1023  // j2 in set C1 (R3 if transpose)
1024  wj[j2] = mark;
1025 
1026  // add j2 to queue
1027  queue[tail++] = j2;
1028  }
1029  }
1030  }
casadi::SparsityInternal::T
Sparsity T() const
Transpose the matrix.
Definition: sparsity_internal.cpp:501

References casadi::Sparsity::colind(), colind(), casadi::Sparsity::row(), row(), and T().

Referenced by dmperm().

◆ btf() [1/2]

const SparsityInternal::Btf & casadi::SparsityInternal::btf ( ) const

Definition at line 463 of file sparsity_internal.cpp.

463  {
464 #ifdef CASADI_WITH_THREADSAFE_SYMBOLICS
465  // Safe access to btf_
466  std::lock_guard<std::mutex> lock(btf_mtx_);
467 #endif // CASADI_WITH_THREADSAFE_SYMBOLICS
468  if (!btf_) {
469  btf_ = new SparsityInternal::Btf();
470  btf_->nb = btf(btf_->rowperm, btf_->colperm, btf_->rowblock, btf_->colblock,
471  btf_->coarse_rowblock, btf_->coarse_colblock);
472  }
473  return *btf_;
474  }
casadi::SparsityInternal::btf
const Btf & btf() const
Get cached block triangular form.
Definition: sparsity_internal.cpp:463

Referenced by spsolve().

◆ btf() [2/2]

casadi_int casadi::SparsityInternal::btf ( std::vector< casadi_int > &  rowperm,
std::vector< casadi_int > &  colperm,
std::vector< casadi_int > &  rowblock,
std::vector< casadi_int > &  colblock,
std::vector< casadi_int > &  coarse_rowblock,
std::vector< casadi_int > &  coarse_colblock 
) const
inline

Definition at line 362 of file sparsity_internal.hpp.

365  {
366  T()->dmperm(colperm, rowperm, colblock, rowblock,
367  coarse_colblock, coarse_rowblock);
368  return rowblock.size()-1;
369  }
casadi::SparsityInternal::dmperm
void dmperm(std::vector< casadi_int > &rowperm, std::vector< casadi_int > &colperm, std::vector< casadi_int > &rowblock, std::vector< casadi_int > &colblock, std::vector< casadi_int > &coarse_rowblock, std::vector< casadi_int > &coarse_colblock) const
Compute the Dulmage-Mendelsohn decomposition.
Definition: sparsity_internal.cpp:1244

References casadi::T.

◆ bw_lower()

casadi_int casadi::SparsityInternal::bw_lower ( ) const

Extra doc: https://github.com/casadi/casadi/wiki/L_ff

Definition at line 3822 of file sparsity_internal.cpp.

3822  {
3823  casadi_int bw = 0;
3824  const casadi_int* colind = this->colind();
3825  const casadi_int* row = this->row();
3826  for (casadi_int cc=0; cc<size2(); ++cc) {
3827  if (colind[cc] != colind[cc+1]) { // if there are any elements of the column
3828  casadi_int rr = row[colind[cc+1]-1];
3829  bw = std::max(bw, rr-cc);
3830  }
3831  }
3832  return bw;
3833  }

References colind(), row(), and size2().

◆ bw_upper()

casadi_int casadi::SparsityInternal::bw_upper ( ) const

Extra doc: https://github.com/casadi/casadi/wiki/L_fe

Definition at line 3809 of file sparsity_internal.cpp.

3809  {
3810  casadi_int bw = 0;
3811  const casadi_int* colind = this->colind();
3812  const casadi_int* row = this->row();
3813  for (casadi_int cc=0; cc<size2(); ++cc) {
3814  if (colind[cc] != colind[cc+1]) { // if there are any elements of the column
3815  casadi_int rr = row[colind[cc]];
3816  bw = std::max(bw, cc-rr);
3817  }
3818  }
3819  return bw;
3820  }

References colind(), row(), and size2().

◆ class_name()

std::string casadi::SparsityInternal::class_name ( ) const
inlineoverridevirtual

Implements casadi::SharedObjectInternal.

Definition at line 709 of file sparsity_internal.hpp.

709 {return "SparsityInternal";}

◆ colind()

const casadi_int* casadi::SparsityInternal::colind ( ) const
inline

Extra doc: https://github.com/casadi/casadi/wiki/L_ei

Definition at line 105 of file sparsity_internal.hpp.

105 { return &sp_.front()+2;}

Referenced by _mtimes(), _reshape(), _resize(), _tril(), _triu(), amd(), augment(), bfs(), bw_lower(), bw_upper(), combineGen(), dfs(), drop(), drop_diag(), etree(), export_code(), find(), get_col(), get_colind(), get_diag(), get_lower(), get_nz(), get_upper(), has_diag(), hash(), is_diag(), is_equal(), is_reshape(), is_stacked(), is_transpose(), is_tril(), is_triu(), ldl_colind(), ldl_row(), makeDense(), maxtrans(), nnz_diag(), nnz_lower(), nnz_upper(), pattern_inverse(), permute(), qr_nnz(), qr_sparsities(), rowsSequential(), scatter(), SparsityInternal(), spsolve(), spy(), star_coloring(), star_coloring2(), sub(), and uni_coloring().

◆ combine() [1/2]

Sparsity casadi::SparsityInternal::combine ( const Sparsity &  y,
bool  f0x_is_zero,
bool  function0_is_zero 
) const

Union of two sparsity patterns

Definition at line 2320 of file sparsity_internal.cpp.

2321  {
2322  static std::vector<unsigned char> mapping;
2323  return combineGen1<false>(y, f0x_is_zero, function0_is_zero, mapping);
2324  }

◆ combine() [2/2]

Sparsity casadi::SparsityInternal::combine ( const Sparsity &  y,
bool  f0x_is_zero,
bool  function0_is_zero,
std::vector< unsigned char > &  mapping 
) const

Union of two sparsity patterns

Definition at line 2326 of file sparsity_internal.cpp.

2328  {
2329  return combineGen1<true>(y, f0x_is_zero, function0_is_zero, mapping);
2330  }

Referenced by casadi::Sparsity::combine(), casadi::Sparsity::intersect(), and casadi::Sparsity::unite().

◆ combineGen()

template<bool with_mapping, bool f0x_is_zero, bool function0_is_zero>
Sparsity casadi::SparsityInternal::combineGen ( const Sparsity &  y,
std::vector< unsigned char > &  mapping 
) const

Union of two sparsity patterns

Definition at line 2370 of file sparsity_internal.cpp.

2371  {
2372 
2373  // Assert dimensions
2374  casadi_assert(size2()==y.size2() && size1()==y.size1(),
2375  "Dimension mismatch : " + str(size()) + " versus " + str(y.size()) + ".");
2376 
2377  // Sparsity pattern of the argument
2378  const casadi_int* y_colind = y.colind();
2379  const casadi_int* y_row = y.row();
2380  const casadi_int* colind = this->colind();
2381  const casadi_int* row = this->row();
2382 
2383  // Sparsity pattern of the result
2384  std::vector<casadi_int> ret_colind(size2()+1, 0);
2385  std::vector<casadi_int> ret_row;
2386 
2387  // Clear the mapping
2388  if (with_mapping) mapping.clear();
2389 
2390  // Loop over columns of both patterns
2391  for (casadi_int i=0; i<size2(); ++i) {
2392  // Non-zero element of the two matrices
2393  casadi_int el1 = colind[i];
2394  casadi_int el2 = y_colind[i];
2395 
2396  // End of the non-zero elements of the col for the two matrices
2397  casadi_int el1_last = colind[i+1];
2398  casadi_int el2_last = y_colind[i+1];
2399 
2400  // Loop over the non-zeros of both matrices
2401  while (el1<el1_last || el2<el2_last) {
2402  // Get the rows
2403  casadi_int row1 = el1<el1_last ? row[el1] : size1();
2404  casadi_int row2 = el2<el2_last ? y_row[el2] : size1();
2405 
2406  // Add to the return matrix
2407  if (row1==row2) { // both nonzero
2408  ret_row.push_back(row1);
2409  if (with_mapping) mapping.push_back( 1 | 2);
2410  el1++; el2++;
2411  } else if (row1<row2) { // only first argument is nonzero
2412  if (!function0_is_zero) {
2413  ret_row.push_back(row1);
2414  if (with_mapping) mapping.push_back(1);
2415  } else {
2416  if (with_mapping) mapping.push_back(1 | 4);
2417  }
2418  el1++;
2419  } else { // only second argument is nonzero
2420  if (!f0x_is_zero) {
2421  ret_row.push_back(row2);
2422  if (with_mapping) mapping.push_back(2);
2423  } else {
2424  if (with_mapping) mapping.push_back(2 | 4);
2425  }
2426  el2++;
2427  }
2428  }
2429 
2430  // Save the index of the last nonzero on the col
2431  ret_colind[i+1] = ret_row.size();
2432  }
2433 
2434  // Return cached object
2435  return Sparsity(size1(), size2(), ret_colind, ret_row);
2436  }

References casadi::Sparsity::colind(), colind(), casadi::Sparsity::row(), row(), casadi::Sparsity::size(), size(), casadi::Sparsity::size1(), size1(), casadi::Sparsity::size2(), size2(), and casadi::str().

◆ combineGen1()

template<bool with_mapping>
Sparsity casadi::SparsityInternal::combineGen1 ( const Sparsity &  y,
bool  f0x_is_zero,
bool  function0_is_zero,
std::vector< unsigned char > &  mapping 
) const

Union of two sparsity patterns

Definition at line 2343 of file sparsity_internal.cpp.

2345  {
2346 
2347  // Quick return if identical
2348  if (is_equal(y)) {
2349  if (with_mapping) {
2350  mapping.resize(y.nnz());
2351  std::fill(mapping.begin(), mapping.end(), 1 | 2);
2352  }
2353  return y;
2354  }
2355 
2356  if (f0x_is_zero) {
2357  if (function0_is_zero) {
2358  return combineGen<with_mapping, true, true>(y, mapping);
2359  } else {
2360  return combineGen<with_mapping, true, false>(y, mapping);
2361  }
2362  } else if (function0_is_zero) {
2363  return combineGen<with_mapping, false, true>(y, mapping);
2364  } else {
2365  return combineGen<with_mapping, false, false>(y, mapping);
2366  }
2367  }
casadi::SparsityInternal::is_equal
bool is_equal(const Sparsity &y) const
Check if two sparsity patterns are the same.
Definition: sparsity_internal.cpp:2474

References is_equal(), casadi::Sparsity::nnz(), and casadi::Sparsity::resize().

◆ debug_repr()

std::string casadi::GenericSharedInternal< SharedObject , SharedObjectInternal >::debug_repr ( const Internal *  i) const
inherited

Definition at line 62 of file generic_shared_internal.hpp.

162  {
163  // Note: i != this because of something something multiple inheritance
164  return str( (casadi_int)(i)) + "/" + static_cast<const Internal*>(this)->class_name();
165  }

◆ destroySingleton()

void casadi::GenericSharedInternal< SharedObject , SharedObjectInternal >::destroySingleton ( )
inlineprotectedinherited

Called in the destructor of singletons

Definition at line 77 of file generic_shared_internal.hpp.

77  {
78  static_cast<Internal*>(this)->count--;
79  }

◆ dfs()

casadi_int casadi::SparsityInternal::dfs ( casadi_int  j,
casadi_int  top,
std::vector< casadi_int > &  xi,
std::vector< casadi_int > &  pstack,
const std::vector< casadi_int > &  pinv,
std::vector< bool > &  marked 
) const

The implementation is a modified version of cs_dfs in CSparse Copyright(c) Timothy A. Davis, 2006-2009 Licensed as a derivative work under the GNU LGPL

Extra doc: https://github.com/casadi/casadi/wiki/L_el

Definition at line 518 of file sparsity_internal.cpp.

521  {
522  /*
523  Modified version of cs_dfs in CSparse
524  Copyright(c) Timothy A. Davis, 2006-2009
525  Licensed as a derivative work under the GNU LGPL
526  */
527  casadi_int head = 0;
528  const casadi_int* colind = this->colind();
529  const casadi_int* row = this->row();
530 
531  // initialize the recursion stack
532  xi[0] = j;
533  while (head >= 0) {
534 
535  // get j from the top of the recursion stack
536  j = xi[head];
537  casadi_int jnew = !pinv.empty() ? (pinv[j]) : j;
538  if (!marked[j]) {
539 
540  // mark node j as visited
541  marked[j]=true;
542  pstack[head] = (jnew < 0) ? 0 : colind[jnew];
543  }
544 
545  // node j done if no unvisited neighbors
546  casadi_int done = 1;
547  casadi_int p2 = (jnew < 0) ? 0 : colind[jnew+1];
548 
549  // examine all neighbors of j
550  for (casadi_int p = pstack[head]; p< p2; ++p) {
551 
552  // consider neighbor node i
553  casadi_int i = row[p];
554 
555  // skip visited node i
556  if (marked[i]) continue ;
557 
558  // pause depth-first search of node j
559  pstack[head] = p;
560 
561  // start dfs at node i
562  xi[++head] = i;
563 
564  // node j is not done
565  done = 0;
566 
567  // break, to start dfs (i)
568  break;
569  }
570 
571  //depth-first search at node j is done
572  if (done) {
573  // remove j from the recursion stack
574  head--;
575 
576  // and place in the output stack
577  xi[--top] = j ;
578  }
579  }
580  return (top) ;
581  }

References colind(), and row().

Referenced by scc().

◆ diag()

casadi_int casadi::SparsityInternal::diag ( casadi_int  i,
casadi_int  j,
double  aij,
void *  other 
)
static

The implementation is a modified version of cs_diag in CSparse Copyright(c) Timothy A. Davis, 2006-2009 Licensed as a derivative work under the GNU LGPL

Extra doc: https://github.com/casadi/casadi/wiki/L_fb

Definition at line 1540 of file sparsity_internal.cpp.

1540  {
1541  /*
1542  Modified version of cs_diag in CSparse
1543  Copyright(c) Timothy A. Davis, 2006-2009
1544  Licensed as a derivative work under the GNU LGPL
1545  */
1546  return (i != j) ;
1547  }

◆ dim()

std::string casadi::SparsityInternal::dim ( bool  with_nz = false) const

Definition at line 1703 of file sparsity_internal.cpp.

1703  {
1704  std::string ret = str(size1()) + "x" + str(size2());
1705  if (with_nz) ret += "," + str(nnz()) + "nz";
1706  return ret;
1707  }

References nnz(), size1(), size2(), and casadi::str().

Referenced by _reshape(), disp(), star_coloring(), and star_coloring2().

◆ disp()

void casadi::SparsityInternal::disp ( std::ostream &  stream,
bool  more 
) const
overridevirtual

Implements casadi::SharedObjectInternal.

Definition at line 481 of file sparsity_internal.cpp.

481  {
482  stream << dim(!is_dense());
483  if (more) {
484  stream << std::endl;
485  stream << "colind: " << get_colind() << std::endl;
486  stream << "row: " << get_row() << std::endl;
487  }
488  }

References dim(), get_colind(), get_row(), and is_dense().

◆ dmperm()

void casadi::SparsityInternal::dmperm ( std::vector< casadi_int > &  rowperm,
std::vector< casadi_int > &  colperm,
std::vector< casadi_int > &  rowblock,
std::vector< casadi_int > &  colblock,
std::vector< casadi_int > &  coarse_rowblock,
std::vector< casadi_int > &  coarse_colblock 
) const

The implementation is a modified version of cs_dmperm in CSparse Copyright(c) Timothy A. Davis, 2006-2009 Licensed as a derivative work under the GNU LGPL

Extra doc: https://github.com/casadi/casadi/wiki/L_f4

Definition at line 1244 of file sparsity_internal.cpp.

1249  {
1250  /*
1251  Modified version of cs_dmperm in CSparse
1252  Copyright(c) Timothy A. Davis, 2006-2009
1253  Licensed as a derivative work under the GNU LGPL
1254  */
1255  casadi_int seed = 0;
1256 
1257  // The transpose of the expression
1258  Sparsity trans;
1259 
1260  // Part 1: Maximum matching
1261 
1262  // col permutation
1263  rowperm.resize(size1());
1264 
1265  // row permutation
1266  colperm.resize(size2());
1267 
1268  // size nb+1, block k is columns r[k] to r[k+1]-1 in A(p, q)
1269  rowblock.resize(size1()+6);
1270 
1271  // size nb+1, block k is rows s[k] to s[k+1]-1 in A(p, q)
1272  colblock.resize(size2()+6);
1273 
1274  // coarse col decomposition
1275  coarse_rowblock.resize(5);
1276  std::fill(coarse_rowblock.begin(), coarse_rowblock.end(), 0);
1277 
1278  // coarse row decomposition
1279  coarse_colblock.resize(5);
1280  std::fill(coarse_colblock.begin(), coarse_colblock.end(), 0);
1281 
1282  // max transversal
1283  std::vector<casadi_int> imatch, jmatch;
1284  maxtrans(imatch, jmatch, trans, seed);
1285 
1286  // Coarse decomposition
1287 
1288  // use rowblock and colblock as workspace
1289  std::vector<casadi_int>& wi = rowblock;
1290  std::vector<casadi_int>& wj = colblock;
1291 
1292  // unmark all rows for bfs
1293  for (casadi_int j=0; j<size2(); ++j)
1294  wj[j] = -1;
1295 
1296  // unmark all columns for bfs
1297  for (casadi_int i=0; i<size1(); ++i)
1298  wi[i] = -1 ;
1299 
1300  // find C1, R1 from C0
1301  bfs(size2(), wi, wj, colperm, imatch, jmatch, 1);
1302 
1303  // find R3, C3 from R0
1304  bfs(size1(), wj, wi, rowperm, jmatch, imatch, 3);
1305 
1306  // unmatched set C0
1307  unmatched(size2(), wj, colperm, coarse_colblock, 0);
1308 
1309  // set R1 and C1
1310  matched(size2(), wj, imatch, rowperm, colperm, coarse_colblock, coarse_rowblock, 1, 1);
1311 
1312  // set R2 and C2
1313  matched(size2(), wj, imatch, rowperm, colperm, coarse_colblock, coarse_rowblock, 2, -1);
1314 
1315  // set R3 and C3
1316  matched(size2(), wj, imatch, rowperm, colperm, coarse_colblock, coarse_rowblock, 3, 3);
1317 
1318  // unmatched set R0
1319  unmatched(size1(), wi, rowperm, coarse_rowblock, 3);
1320 
1321  // --- Fine decomposition -----------------------------------------------
1322  // pinv=p'
1323  std::vector<casadi_int> pinv = invertPermutation(rowperm);
1324 
1325  // C=A(p, q) (it will hold A(R2, C2))
1326  std::vector<casadi_int> colind_C, row_C;
1327  permute(pinv, colperm, 0, colind_C, row_C);
1328 
1329  // delete rows C0, C1, and C3 from C
1330  casadi_int nc = coarse_colblock[3] - coarse_colblock[2];
1331  if (coarse_colblock[2] > 0) {
1332  for (casadi_int j = coarse_colblock[2]; j <= coarse_colblock[3]; ++j)
1333  colind_C[j-coarse_colblock[2]] = colind_C[j];
1334  }
1335  casadi_int ncol_C = nc;
1336 
1337  colind_C.resize(nc+1);
1338  // delete columns R0, R1, and R3 from C
1339  if (coarse_rowblock[2] - coarse_rowblock[1] < size1()) {
1340  drop(rprune, &coarse_rowblock, size1(), ncol_C, colind_C, row_C);
1341  casadi_int cnz = colind_C[nc];
1342  if (coarse_rowblock[1] > 0)
1343  for (casadi_int k=0; k<cnz; ++k)
1344  row_C[k] -= coarse_rowblock[1];
1345  }
1346  row_C.resize(colind_C.back());
1347  casadi_int nrow_C = nc ;
1348  Sparsity C(nrow_C, ncol_C, colind_C, row_C, true);
1349 
1350  // find strongly connected components of C
1351  std::vector<casadi_int> scc_p, scc_r;
1352  casadi_int scc_nb = C.scc(scc_p, scc_r);
1353 
1354  // --- Combine coarse and fine decompositions ---------------------------
1355 
1356  // C(ps, ps) is the permuted matrix
1357  std::vector<casadi_int> ps = scc_p;
1358 
1359  // kth block is rs[k]..rs[k+1]-1
1360  std::vector<casadi_int> rs = scc_r;
1361 
1362  // # of blocks of A(R2, C2)
1363  casadi_int nb1 = scc_nb;
1364 
1365  for (casadi_int k=0; k<nc; ++k)
1366  wj[k] = colperm[ps[k] + coarse_colblock[2]];
1367 
1368  for (casadi_int k=0; k<nc; ++k)
1369  colperm[k + coarse_colblock[2]] = wj[k];
1370 
1371  for (casadi_int k=0; k<nc; ++k)
1372  wi[k] = rowperm[ps[k] + coarse_rowblock[1]];
1373 
1374  for (casadi_int k=0; k<nc; ++k)
1375  rowperm[k + coarse_rowblock[1]] = wi[k];
1376 
1377  // create the fine block partitions
1378  casadi_int nb2 = 0;
1379  rowblock[0] = colblock[0] = 0;
1380 
1381  // leading coarse block A (R1, [C0 C1])
1382  if (coarse_colblock[2] > 0)
1383  nb2++ ;
1384 
1385  // coarse block A (R2, C2)
1386  for (casadi_int k=0; k<nb1; ++k) {
1387  // A (R2, C2) splits into nb1 fine blocks
1388  rowblock[nb2] = rs[k] + coarse_rowblock[1];
1389  colblock[nb2] = rs[k] + coarse_colblock[2] ;
1390  nb2++ ;
1391  }
1392 
1393  if (coarse_rowblock[2] < size1()) {
1394  // trailing coarse block A ([R3 R0], C3)
1395  rowblock[nb2] = coarse_rowblock[2];
1396  colblock[nb2] = coarse_colblock[3];
1397  nb2++ ;
1398  }
1399 
1400  rowblock[nb2] = size1();
1401  colblock[nb2] = size2() ;
1402 
1403  // Shrink rowblock and colblock
1404  rowblock.resize(nb2+1);
1405  colblock.resize(nb2+1);
1406  }
casadi::SparsityInternal::drop
static casadi_int drop(casadi_int(*fkeep)(casadi_int, casadi_int, double, void *), void *other, casadi_int nrow, casadi_int ncol, std::vector< casadi_int > &colind, std::vector< casadi_int > &row)
drop entries for which fkeep(A(i, j)) is false; return nz if OK, else -1
Definition: sparsity_internal.cpp:1496
casadi::SparsityInternal::bfs
void bfs(casadi_int n, std::vector< casadi_int > &wi, std::vector< casadi_int > &wj, std::vector< casadi_int > &queue, const std::vector< casadi_int > &imatch, const std::vector< casadi_int > &jmatch, casadi_int mark) const
Breadth-first search for coarse decomposition.
Definition: sparsity_internal.cpp:967
casadi::SparsityInternal::maxtrans
void maxtrans(std::vector< casadi_int > &imatch, std::vector< casadi_int > &jmatch, Sparsity &trans, casadi_int seed) const
Compute the maximum transversal (maximum matching)
Definition: sparsity_internal.cpp:1154
casadi::SparsityInternal::unmatched
static void unmatched(casadi_int m, const std::vector< casadi_int > &wi, std::vector< casadi_int > &p, std::vector< casadi_int > &rr, casadi_int set)
Collect unmatched columns into the permutation vector p.
Definition: sparsity_internal.cpp:1055
casadi::SparsityInternal::rprune
static casadi_int rprune(casadi_int i, casadi_int j, double aij, void *other)
return 1 if column i is in R2
Definition: sparsity_internal.cpp:1070
casadi::SparsityInternal::permute
Sparsity permute(const std::vector< casadi_int > &pinv, const std::vector< casadi_int > &q, casadi_int values) const
C = A(p, q) where p and q are permutations of 0..m-1 and 0..n-1.
Definition: sparsity_internal.cpp:1456
casadi::SparsityInternal::invertPermutation
static std::vector< casadi_int > invertPermutation(const std::vector< casadi_int > &p)
Invert a permutation vector.
Definition: sparsity_internal.cpp:1450
casadi::SparsityInternal::matched
static void matched(casadi_int n, const std::vector< casadi_int > &wj, const std::vector< casadi_int > &imatch, std::vector< casadi_int > &p, std::vector< casadi_int > &q, std::vector< casadi_int > &cc, std::vector< casadi_int > &rr, casadi_int set, casadi_int mark)
Collect matched columns and rows into p and q.
Definition: sparsity_internal.cpp:1032

References bfs(), casadi::C, drop(), invertPermutation(), matched(), maxtrans(), permute(), casadi::Sparsity::resize(), rprune(), size1(), size2(), and unmatched().

◆ drop()

casadi_int casadi::SparsityInternal::drop ( casadi_int(*)(casadi_int, casadi_int, double, void *)  fkeep,
void *  other,
casadi_int  nrow,
casadi_int  ncol,
std::vector< casadi_int > &  colind,
std::vector< casadi_int > &  row 
)
static

The implementation is a modified version of cs_fkeep in CSparse Copyright(c) Timothy A. Davis, 2006-2009 Licensed as a derivative work under the GNU LGPL

Extra doc: https://github.com/casadi/casadi/wiki/L_f3

Definition at line 1496 of file sparsity_internal.cpp.

1498  {
1499  /*
1500  Modified version of cs_fkeep in CSparse
1501  Copyright(c) Timothy A. Davis, 2006-2009
1502  Licensed as a derivative work under the GNU LGPL
1503  */
1504  casadi_int nz = 0;
1505 
1506  for (casadi_int j = 0; j<ncol; ++j) {
1507  // get current location of row j
1508  casadi_int p = colind[j];
1509 
1510  // record new location of row j
1511  colind[j] = nz;
1512  for ( ; p < colind[j+1] ; ++p) {
1513  if (fkeep(row[p], j, 1, other)) {
1514  // keep A(i, j)
1515  row[nz++] = row[p] ;
1516  }
1517  }
1518  }
1519 
1520  // finalize A
1521  colind[ncol] = nz;
1522  return nz ;
1523  }

References colind(), and row().

Referenced by dmperm().

◆ drop_diag()

Sparsity casadi::SparsityInternal::drop_diag ( ) const

Definition at line 1683 of file sparsity_internal.cpp.

1683  {
1684  casadi_int nrow = this->size1();
1685  casadi_int ncol = this->size2();
1686  const casadi_int* colind = this->colind();
1687  const casadi_int* row = this->row();
1688  // Return sparsity
1689  std::vector<casadi_int> ret_colind(ncol+1), ret_row;
1690  ret_colind[0] = 0;
1691  ret_row.reserve(nnz());
1692  for (casadi_int c=0; c<ncol; ++c) {
1693  for (casadi_int k=colind[c]; k<colind[c+1]; ++k) {
1694  if (row[k]!=c) {
1695  ret_row.push_back(row[k]);
1696  }
1697  }
1698  ret_colind[c+1] = ret_row.size();
1699  }
1700  return Sparsity(nrow, ncol, ret_colind, ret_row);
1701  }

References colind(), nnz(), row(), size1(), and size2().

◆ etree()

void casadi::SparsityInternal::etree ( const casadi_int *  sp,
casadi_int *  parent,
casadi_int *  w,
casadi_int  ata 
)
static

len[w] >= ata ? ncol + nrow : ncol len[parent] == ncol Ref: Chapter 4, Direct Methods for Sparse Linear Systems by Tim Davis Modified version of cs_etree in CSparse Copyright(c) Timothy A. Davis, 2006-2009 Licensed as a derivative work under the GNU LGPL

Extra doc: https://github.com/casadi/casadi/wiki/L_eo

Definition at line 35 of file sparsity_internal.cpp.

36  {
37  /*
38  Modified version of cs_etree in CSparse
39  Copyright(c) Timothy A. Davis, 2006-2009
40  Licensed as a derivative work under the GNU LGPL
41  */
42  casadi_int r, c, k, rnext;
43  // Extract sparsity
44  casadi_int nrow = *sp++, ncol = *sp++;
45  const casadi_int *colind = sp, *row = sp+ncol+1;
46  // Highest known ascestor of a node
47  casadi_int *ancestor=w;
48  // Path for A'A
49  casadi_int *prev;
50  if (ata) {
51  prev=w+ncol;
52  for (r=0; r<nrow; ++r) prev[r] = -1;
53  }
54  // Loop over columns
55  for (c=0; c<ncol; ++c) {
56  parent[c] = -1; // No parent yet
57  ancestor[c] = -1; // No ancestor
58  // Loop over nonzeros
59  for (k=colind[c]; k<colind[c+1]; ++k) {
60  r = row[k];
61  if (ata) r = prev[r];
62  // Traverse from r to c
63  while (r!=-1 && r<c) {
64  rnext = ancestor[r];
65  ancestor[r] = c;
66  if (rnext==-1) parent[r] = c;
67  r = rnext;
68  }
69  if (ata) prev[row[k]] = c;
70  }
71  }
72  }

References colind(), row(), and sp().

Referenced by casadi::Sparsity::etree(), and qr_init().

◆ export_code()

void casadi::SparsityInternal::export_code ( const std::string &  lang,
std::ostream &  stream,
const Dict &  options 
) const

Extra doc: https://github.com/casadi/casadi/wiki/L_fr

Definition at line 3607 of file sparsity_internal.cpp.

3608  {
3609  casadi_assert(lang=="matlab", "Only matlab language supported for now.");
3610 
3611  // Default values for options
3612  bool opt_inline = false;
3613  std::string name = "sp";
3614  bool as_matrix = true;
3615  casadi_int indent_level = 0;
3616  std::vector<std::string> nonzeros;
3617 
3618  // Read options
3619  for (auto&& op : options) {
3620  if (op.first=="inline") {
3621  opt_inline = op.second;
3622  } else if (op.first=="name") {
3623  name = op.second.to_string();
3624  } else if (op.first=="as_matrix") {
3625  as_matrix = op.second;
3626  } else if (op.first=="indent_level") {
3627  indent_level = op.second;
3628  } else if (op.first=="nonzeros") {
3629  nonzeros = op.second;
3630  } else {
3631  casadi_error("Unknown option '" + op.first + "'.");
3632  }
3633  }
3634 
3635  // Construct indent string
3636  std::string indent;
3637  for (casadi_int i=0; i < indent_level; ++i) indent += " ";
3638  casadi_assert(!opt_inline, "Inline not supported for now.");
3639 
3640  // Export dimensions
3641  stream << indent << name << "_m = " << size1() << ";\n";
3642  stream << indent << name << "_n = " << size2() << ";\n";
3643 
3644  // Matlab indices are one-based
3645  const casadi_int index_offset = 1;
3646 
3647  // Print columns
3648  const casadi_int* colind = this->colind();
3649  const casadi_int* row = this->row();
3650  stream << indent << name<< "_j = [";
3651  bool first = true;
3652  for (casadi_int i=0; i<size2(); ++i) {
3653  for (casadi_int el=colind[i]; el<colind[i+1]; ++el) {
3654  if (!first) stream << ", ";
3655  stream << (i+index_offset);
3656  first = false;
3657  }
3658  }
3659  stream << "];\n";
3660 
3661  // Print rows
3662  stream << indent << name << "_i = [";
3663  first = true;
3664  casadi_int nz = nnz();
3665  for (casadi_int i=0; i<nz; ++i) {
3666  if (!first) stream << ", ";
3667  stream << (row[i]+index_offset);
3668  first = false;
3669  }
3670  stream << "];\n";
3671 
3672  // Print nonzeros
3673  stream << indent << name << "_v = ";
3674  if (nonzeros.empty()) {
3675  stream << "ones(size(" << name << "_i));\n";
3676  } else {
3677  stream << "[";
3678  for (casadi_int i = 0; i < nonzeros.size(); ++i) {
3679  if (i > 0) stream << ", ";
3680  stream << nonzeros.at(i);
3681  }
3682  stream << "];\n";
3683  }
3684 
3685  if (as_matrix) {
3686  // Generate matrix
3687  stream << indent << name << " = sparse(" << name << "_i, " << name << "_j, ";
3688  stream << name << "_v, " << name << "_m, " << name << "_n);\n";
3689  }
3690 
3691  }

References colind(), nnz(), row(), size1(), and size2().

Referenced by spy_matlab().

◆ find()

void casadi::SparsityInternal::find ( std::vector< casadi_int > &  loc,
bool  ind1 
) const

Definition at line 2840 of file sparsity_internal.cpp.

2840  {
2841  casadi_assert(!std::mul_overflows(size1(), size2()), "Integer overflow detected");
2842  if (is_dense()) {
2843  loc = range(ind1, numel()+ind1);
2844  return;
2845  }
2846  const casadi_int* colind = this->colind();
2847  const casadi_int* row = this->row();
2848 
2849  // Element for each nonzero
2850  loc.resize(nnz());
2851 
2852  // Loop over columns
2853  for (casadi_int cc=0; cc<size2(); ++cc) {
2854 
2855  // Loop over the nonzeros
2856  for (casadi_int el=colind[cc]; el<colind[cc+1]; ++el) {
2857 
2858  // Get row
2859  casadi_int rr = row[el];
2860 
2861  // Get the element
2862  loc[el] = rr+cc*size1()+ind1;
2863  }
2864  }
2865  }
std::mul_overflows
bool mul_overflows(const T &a, const T &b)
Definition: casadi_misc.hpp:523

References colind(), is_dense(), std::mul_overflows(), nnz(), numel(), casadi::range(), row(), size1(), and size2().

◆ get_col()

std::vector< casadi_int > casadi::SparsityInternal::get_col ( ) const

Definition at line 490 of file sparsity_internal.cpp.

490  {
491  const casadi_int* colind = this->colind();
492  std::vector<casadi_int> col(nnz());
493  for (casadi_int r=0; r<size2(); ++r) {
494  for (casadi_int el = colind[r]; el < colind[r+1]; ++el) {
495  col[el] = r;
496  }
497  }
498  return col;
499  }

References colind(), nnz(), and size2().

Referenced by pmult(), repr_el(), and transpose().

◆ get_colind()

std::vector< casadi_int > casadi::SparsityInternal::get_colind ( ) const

Definition at line 3835 of file sparsity_internal.cpp.

3835  {
3836  const casadi_int* colind = this->colind();
3837  return std::vector<casadi_int>(colind, colind+size2()+1);
3838  }

References colind(), and size2().

Referenced by _appendColumns(), _enlargeColumns(), _enlargeRows(), _erase(), _removeDuplicates(), amd(), disp(), largest_first(), pattern_inverse(), and spy().

◆ get_diag()

Sparsity casadi::SparsityInternal::get_diag ( std::vector< casadi_int > &  mapping) const
Parameters
[out]mappingwill contain the nonzero mapping

Extra doc: https://github.com/casadi/casadi/wiki/L_ek

Definition at line 1619 of file sparsity_internal.cpp.

1619  {
1620  casadi_int nrow = this->size1();
1621  casadi_int ncol = this->size2();
1622  const casadi_int* colind = this->colind();
1623  const casadi_int* row = this->row();
1624 
1625  // Mapping
1626  mapping.clear();
1627 
1628  if (is_vector()) {
1629  // Sparsity pattern
1630  casadi_int n = nrow * ncol;
1631  std::vector<casadi_int> ret_colind(n+1, 0), ret_row;
1632 
1633  // Loop over all entries
1634  casadi_int ret_i=0;
1635  for (casadi_int cc=0; cc<ncol; ++cc) {
1636  for (casadi_int k = colind[cc]; k<colind[cc+1]; ++k) {
1637  casadi_int rr=row[k];
1638  casadi_int el=rr+nrow*cc; // Corresponding row in the return matrix
1639  while (ret_i<=el) ret_colind[ret_i++]=ret_row.size();
1640  ret_row.push_back(el);
1641  mapping.push_back(k);
1642  }
1643  }
1644  while (ret_i<=n) ret_colind[ret_i++]=ret_row.size();
1645 
1646  // Construct sparsity pattern
1647  return Sparsity(n, n, ret_colind, ret_row);
1648 
1649  } else {
1650  // Sparsity pattern
1651  casadi_int n = std::min(nrow, ncol);
1652  std::vector<casadi_int> ret_row, ret_colind(2, 0);
1653 
1654  // Loop over diagonal nonzeros
1655  for (casadi_int cc=0; cc<n; ++cc) {
1656  for (casadi_int el = colind[cc]; el<colind[cc+1]; ++el) {
1657  if (row[el]==cc) {
1658  ret_row.push_back(row[el]);
1659  ret_colind[1]++;
1660  mapping.push_back(el);
1661  }
1662  }
1663  }
1664 
1665  // Construct sparsity pattern
1666  return Sparsity(n, 1, ret_colind, ret_row);
1667  }
1668  }
casadi::SparsityInternal::is_vector
bool is_vector() const
Check if the pattern is a row or column vector.
Definition: sparsity_internal.cpp:1797

References colind(), is_vector(), row(), size1(), and size2().

Referenced by casadi::Sparsity::get_diag().

◆ get_lower()

std::vector< casadi_int > casadi::SparsityInternal::get_lower ( ) const

Definition at line 3783 of file sparsity_internal.cpp.

3783  {
3784  const casadi_int* colind = this->colind();
3785  const casadi_int* row = this->row();
3786  std::vector<casadi_int> ret;
3787  for (casadi_int cc=0; cc<size2(); ++cc) {
3788  for (casadi_int el = colind[cc]; el<colind[cc+1]; ++el) {
3789  if (row[el]>=cc) {
3790  ret.push_back(el);
3791  }
3792  }
3793  }
3794  return ret;
3795  }

References colind(), row(), and size2().

◆ get_nz() [1/3]

casadi_int casadi::SparsityInternal::get_nz ( casadi_int  rr,
casadi_int  cc 
) const

Definition at line 2676 of file sparsity_internal.cpp.

2676  {
2677  // If negative index, count from the back
2678  if (rr<0) rr += size1();
2679  if (cc<0) cc += size2();
2680  const casadi_int* colind = this->colind();
2681  const casadi_int* row = this->row();
2682 
2683  // Check consistency
2684  casadi_assert(rr>=0 && rr<size1(), "Row index " + str(rr)
2685  + " out of bounds [0, " + str(size1()) + ")");
2686  casadi_assert(cc>=0 && cc<size2(), "Column index " + str(cc)
2687  + " out of bounds [0, " + str(size2()) + ")");
2688 
2689  // Quick return if matrix is dense
2690  if (is_dense()) return rr+cc*size1();
2691 
2692  // Quick return if past the end
2693  if (colind[cc]==nnz() || (colind[cc+1]==nnz() && row[nnz()-1]<rr)) return -1;
2694 
2695  // Find sparse element
2696  for (casadi_int ind=colind[cc]; ind<colind[cc+1]; ++ind) {
2697  if (row[ind] == rr) {
2698  return ind; // element exists
2699  } else if (row[ind] > rr) {
2700  break; // break at the place where the element should be added
2701  }
2702  }
2703  return -1;
2704  }

References colind(), is_dense(), nnz(), row(), size1(), size2(), and casadi::str().

Referenced by get_nz(), and sub().

◆ get_nz() [2/3]

std::vector< casadi_int > casadi::SparsityInternal::get_nz ( const std::vector< casadi_int > &  rr,
const std::vector< casadi_int > &  cc 
) const

Definition at line 2118 of file sparsity_internal.cpp.

2119  {
2120  casadi_assert_bounded(rr, size1());
2121  casadi_assert_bounded(cc, size2());
2122 
2123  std::vector<casadi_int> rr_sorted;
2124  std::vector<casadi_int> rr_sorted_index;
2125 
2126  sort(rr, rr_sorted, rr_sorted_index);
2127 
2128  std::vector<casadi_int> ret(cc.size()*rr.size());
2129 
2130  casadi_int stride = rr.size();
2131  const casadi_int* colind = this->colind();
2132  const casadi_int* row = this->row();
2133 
2134  for (casadi_int i=0;i<cc.size();++i) {
2135  casadi_int it = cc[i];
2136  casadi_int el=colind[it];
2137  for (casadi_int j=0;j<rr_sorted.size();++j) {
2138  casadi_int jt=rr_sorted[j];
2139  // Continue to the non-zero element
2140  for (; el<colind[it+1] && row[el]<jt; ++el) {}
2141  // Add the non-zero element, if there was an element in the location exists
2142  if (el<colind[it+1] && row[el]== jt) {
2143  ret[i*stride+rr_sorted_index[j]] = el;
2144  } else {
2145  ret[i*stride+rr_sorted_index[j]] = -1;
2146  }
2147  }
2148  }
2149  return ret;
2150  }
casadi::sort
void sort(const std::vector< T > &values, std::vector< T > &sorted_values, std::vector< casadi_int > &indices, bool invert_indices=false)
Sort the data in a vector.
Definition: casadi_misc.hpp:881

References colind(), row(), size1(), size2(), and casadi::sort().

◆ get_nz() [3/3]

void casadi::SparsityInternal::get_nz ( std::vector< casadi_int > &  indices) const

Definition at line 2867 of file sparsity_internal.cpp.

2867  {
2868  // Quick return if no elements
2869  if (indices.empty()) return;
2870  // In a dense matrix, all indices can be found
2871  if (is_dense()) return;
2872  const casadi_int* colind = this->colind();
2873  const casadi_int* row = this->row();
2874 
2875  // Make a sanity check
2876  casadi_int last=-1;
2877  for (std::vector<casadi_int>::iterator it=indices.begin(); it!=indices.end(); ++it) {
2878  if (*it>=0) {
2879  casadi_int el = *it;
2880  if (el<last) {
2881  // Sort rr in nondecreasing order, if needed
2882  std::vector<casadi_int> indices_sorted, mapping;
2883  sort(indices, indices_sorted, mapping, false);
2884  get_nz(indices_sorted);
2885  for (size_t i=0; i<indices.size(); ++i) {
2886  indices[mapping[i]] = indices_sorted[i];
2887  }
2888  return;
2889  }
2890  last = el;
2891  }
2892  }
2893 
2894  // Quick return if no elements
2895  if (last<0) return;
2896 
2897  // Iterator to input/output
2898  std::vector<casadi_int>::iterator it=indices.begin();
2899  while (*it<0) it++; // first non-ignored
2900 
2901  // Position in flattened matrix
2902  casadi_int cur_pos = -1;
2903 
2904  casadi_int col_pos = 0;
2905 
2906  // Loop over columns
2907  for (casadi_int i=0; i<size2(); ++i, col_pos+=size1()) {
2908 
2909  // Last position in flattened matrix for current column
2910  casadi_int last_pos = -1;
2911 
2912  // Early skip to next column
2913  if (colind[i+1]>colind[i]) {
2914  casadi_int el = colind[i+1] - 1;
2915  casadi_int j = row[el];
2916  last_pos = col_pos + j;
2917  } else {
2918  continue;
2919  }
2920 
2921  // Loop over the nonzeros
2922  for (casadi_int el=colind[i]; el<colind[i+1] && last_pos >= *it; ++el) {
2923  // Get row
2924  casadi_int j = row[el];
2925 
2926  cur_pos = col_pos + j;
2927 
2928  // Add leading elements not in pattern
2929  while (*it < cur_pos) {
2930  // Mark as not found
2931  *it = -1;
2932  if (++it==indices.end()) return;
2933  }
2934 
2935  while (cur_pos == *it) {
2936  // Save element index
2937  *it = el;
2938 
2939  // Increase index and terminate if end of vector reached
2940  do {
2941  if (++it==indices.end()) return;
2942  } while (*it<0);
2943  }
2944  }
2945  }
2946 
2947  // Add trailing elements not in pattern
2948  std::fill(it, indices.end(), -1);
2949  }
casadi::SparsityInternal::get_nz
casadi_int get_nz(casadi_int rr, casadi_int cc) const
Get the index of an existing non-zero element.
Definition: sparsity_internal.cpp:2676

References colind(), get_nz(), is_dense(), row(), size1(), size2(), and casadi::sort().

◆ get_row()

std::vector< casadi_int > casadi::SparsityInternal::get_row ( ) const

Definition at line 3840 of file sparsity_internal.cpp.

3840  {
3841  const casadi_int* row = this->row();
3842  return std::vector<casadi_int>(row, row+nnz());
3843  }

References nnz(), and row().

Referenced by _appendColumns(), _appendVector(), _enlargeColumns(), _enlargeRows(), _erase(), _removeDuplicates(), amd(), disp(), and transpose().

◆ get_upper()

std::vector< casadi_int > casadi::SparsityInternal::get_upper ( ) const

Definition at line 3797 of file sparsity_internal.cpp.

3797  {
3798  const casadi_int* colind = this->colind();
3799  const casadi_int* row = this->row();
3800  std::vector<casadi_int> ret;
3801  for (casadi_int cc=0; cc<size2(); ++cc) {
3802  for (casadi_int el = colind[cc]; el<colind[cc+1] && row[el]<=cc; ++el) {
3803  ret.push_back(el);
3804  }
3805  }
3806  return ret;
3807  }

References colind(), row(), and size2().

◆ getCount()

casadi_int casadi::GenericSharedInternal< SharedObject , SharedObjectInternal >::getCount
inherited

Definition at line 60 of file generic_shared_internal.hpp.

186  {
187  return static_cast<const Internal*>(this)->count;
188  }

◆ has_diag()

bool casadi::SparsityInternal::has_diag ( ) const

Definition at line 1670 of file sparsity_internal.cpp.

1670  {
1671  casadi_int nrow = this->size1();
1672  casadi_int ncol = this->size2();
1673  const casadi_int* colind = this->colind();
1674  const casadi_int* row = this->row();
1675  for (casadi_int c=0; c<ncol && c<nrow; ++c) {
1676  for (casadi_int k=colind[c]; k<colind[c+1]; ++k) {
1677  if (row[k]==c) return true;
1678  }
1679  }
1680  return false;
1681  }

References colind(), row(), size1(), and size2().

◆ hash()

std::size_t casadi::SparsityInternal::hash ( ) const

Definition at line 3711 of file sparsity_internal.cpp.

3711  {
3712  return hash_sparsity(size1(), size2(), colind(), row());
3713  }
casadi::hash_sparsity
std::size_t hash_sparsity(casadi_int nrow, casadi_int ncol, const std::vector< casadi_int > &colind, const std::vector< casadi_int > &row)
Hash a sparsity pattern.
Definition: sparsity.cpp:996

References colind(), casadi::hash_sparsity(), row(), size1(), and size2().

◆ initSingleton()

void casadi::GenericSharedInternal< SharedObject , SharedObjectInternal >::initSingleton ( )
inlineprotectedinherited

Called in the constructor of singletons to avoid that the counter reaches zero

Definition at line 71 of file generic_shared_internal.hpp.

71  {
72  casadi_assert_dev(static_cast<Internal*>(this)->count==0);
73  static_cast<Internal*>(this)->count++;
74  }

◆ invertPermutation()

std::vector< casadi_int > casadi::SparsityInternal::invertPermutation ( const std::vector< casadi_int > &  p)
static

Extra doc: https://github.com/casadi/casadi/wiki/L_f8

Definition at line 1450 of file sparsity_internal.cpp.

1450  {
1451  std::vector<casadi_int> pinv(p.size());
1452  for (casadi_int k=0; k<p.size(); ++k) pinv[p[k]] = k;
1453  return pinv;
1454  }

Referenced by dmperm().

◆ is_column()

bool casadi::SparsityInternal::is_column ( ) const

Extra doc: https://github.com/casadi/casadi/wiki/L_fi

Definition at line 1793 of file sparsity_internal.cpp.

1793  {
1794  return size2()==1;
1795  }

References size2().

Referenced by is_vector().

◆ is_dense()

bool casadi::SparsityInternal::is_dense ( ) const

Definition at line 1785 of file sparsity_internal.cpp.

1785  {
1786  return nnz() == numel();
1787  }

References nnz(), and numel().

Referenced by _mtimes(), disp(), find(), get_nz(), is_reshape(), is_transpose(), and pattern_inverse().

◆ is_diag()

bool casadi::SparsityInternal::is_diag ( ) const

Definition at line 1805 of file sparsity_internal.cpp.

1805  {
1806  const casadi_int* colind = this->colind();
1807  const casadi_int* row = this->row();
1808 
1809  // Check if matrix is square
1810  if (size2() != size1()) return false;
1811 
1812  // Check if correct number of non-zeros (one per column)
1813  if (nnz() != size2()) return false;
1814 
1815  // Check that the row indices are correct
1816  for (casadi_int i=0; i<nnz(); ++i) {
1817  if (row[i]!=i)
1818  return false;
1819  }
1820 
1821  // Make sure that the col indices are correct
1822  for (casadi_int i=0; i<size2(); ++i) {
1823  if (colind[i]!=i)
1824  return false;
1825  }
1826 
1827  // Diagonal if reached this point
1828  return true;
1829  }

References colind(), nnz(), row(), size1(), and size2().

Referenced by _mtimes().

◆ is_empty()

bool casadi::SparsityInternal::is_empty ( bool  both = false) const

A sparsity is considered empty if one of the dimensions is zero (or optionally both dimensions)

Extra doc: https://github.com/casadi/casadi/wiki/L_fg

Definition at line 1801 of file sparsity_internal.cpp.

1801  {
1802  return both ? size2()==0 && size1()==0 : size2()==0 || size1()==0;
1803  }

References size1(), and size2().

Referenced by _mtimes(), and pattern_inverse().

◆ is_equal() [1/3]

bool casadi::SparsityInternal::is_equal ( casadi_int  y_nrow,
casadi_int  y_ncol,
const casadi_int *  y_colind,
const casadi_int *  y_row 
) const

Definition at line 2532 of file sparsity_internal.cpp.

2533  {
2534  const casadi_int* colind = this->colind();
2535  const casadi_int* row = this->row();
2536 
2537  // Get number of nonzeros
2538  casadi_int nz = y_colind[y_ncol];
2539 
2540  // First check dimensions and number of non-zeros
2541  if (nnz()!=nz || size2()!=y_ncol || size1()!=y_nrow) return false;
2542 
2543  // Check if dense
2544  if (nnz()==numel()) return true;
2545 
2546  // Check the number of non-zeros per col
2547  if (!std::equal(colind, colind+size2()+1, y_colind)) return false;
2548 
2549  // Finally check the row indices
2550  if (!std::equal(row, row+nz, y_row)) return false;
2551 
2552  // Equal if reached this point
2553  return true;
2554  }

References colind(), nnz(), numel(), row(), size1(), and size2().

◆ is_equal() [2/3]

bool casadi::SparsityInternal::is_equal ( casadi_int  y_nrow,
casadi_int  y_ncol,
const std::vector< casadi_int > &  y_colind,
const std::vector< casadi_int > &  y_row 
) const

Definition at line 2524 of file sparsity_internal.cpp.

2526  {
2527  casadi_assert_dev(y_colind.size()==y_ncol+1);
2528  casadi_assert_dev(y_row.size()==y_colind.back());
2529  return is_equal(y_nrow, y_ncol, get_ptr(y_colind), get_ptr(y_row));
2530  }

References casadi::get_ptr(), and is_equal().

◆ is_equal() [3/3]

bool casadi::SparsityInternal::is_equal ( const Sparsity &  y) const

Definition at line 2474 of file sparsity_internal.cpp.

2474  {
2475  // Quick true if the objects are the same
2476  if (this == y.get()) return true;
2477 
2478  // Otherwise, compare the patterns
2479  return is_equal(y.size1(), y.size2(), y.colind(), y.row());
2480  }

References casadi::Sparsity::colind(), casadi::Sparsity::get(), casadi::Sparsity::row(), casadi::Sparsity::size1(), and casadi::Sparsity::size2().

Referenced by combineGen1(), is_equal(), is_reshape(), is_stacked(), and is_subset().

◆ is_orthonormal()

bool casadi::SparsityInternal::is_orthonormal ( bool  allow_empty = false) const

Extra doc: https://github.com/casadi/casadi/wiki/L_25x

Definition at line 1843 of file sparsity_internal.cpp.

1843  {
1844  if (!allow_empty) {
1845  if (!is_square()) return false;
1846  if (nnz()!=size1()) return false;
1847  }
1848 
1849  Sparsity sp = shared_from_this<Sparsity>();
1850  if (sum2(sp).nnz()!=nnz()) return false;
1851  if (sum1(sp).nnz()!=nnz()) return false;
1852  return true;
1853  }
casadi::SparsityInternal::is_square
bool is_square() const
Is square?
Definition: sparsity_internal.cpp:1831

References is_square(), nnz(), size1(), and sp().

Referenced by is_permutation().

◆ is_orthonormal_columns()

bool casadi::SparsityInternal::is_orthonormal_columns ( bool  allow_empty = false) const

Extra doc: https://github.com/casadi/casadi/wiki/L_25z

Definition at line 1867 of file sparsity_internal.cpp.

1867  {
1868  if (!allow_empty) {
1869  if (size2()>size1()) return false;
1870  if (nnz()!=size2()) return false;
1871  }
1872 
1873  Sparsity sp = shared_from_this<Sparsity>();
1874  if (sum2(sp).nnz()!=nnz()) return false;
1875  if (sum1(sp).nnz()!=nnz()) return false;
1876  return true;
1877  }

References nnz(), size1(), size2(), and sp().

◆ is_orthonormal_rows()

bool casadi::SparsityInternal::is_orthonormal_rows ( bool  allow_empty = false) const

Extra doc: https://github.com/casadi/casadi/wiki/L_25y

Definition at line 1855 of file sparsity_internal.cpp.

1855  {
1856  if (!allow_empty) {
1857  if (size1()>size2()) return false;
1858  if (nnz()!=size1()) return false;
1859  }
1860 
1861  Sparsity sp = shared_from_this<Sparsity>();
1862  if (sum2(sp).nnz()!=nnz()) return false;
1863  if (sum1(sp).nnz()!=nnz()) return false;
1864  return true;
1865  }

References nnz(), size1(), size2(), and sp().

Referenced by is_selection().

◆ is_permutation()

bool casadi::SparsityInternal::is_permutation ( ) const

Extra doc: https://github.com/casadi/casadi/wiki/L_25v

Definition at line 1839 of file sparsity_internal.cpp.

1839  {
1840  return is_orthonormal();
1841  }
casadi::SparsityInternal::is_orthonormal
bool is_orthonormal(bool allow_empty=false) const
Are the rows and columns of the pattern orthonormal ?
Definition: sparsity_internal.cpp:1843

References is_orthonormal().

◆ is_reshape()

bool casadi::SparsityInternal::is_reshape ( const SparsityInternal &  y) const

Definition at line 3542 of file sparsity_internal.cpp.

3542  {
3543  // Quick return if the objects are the same
3544  if (this==&y) return true;
3545 
3546  // Check if same number of entries and nonzeros
3547  if (numel()!=y.numel() || nnz()!=y.nnz()) return false;
3548 
3549  // Quick return if empty interior or dense
3550  if (nnz()==0 || is_dense()) return true;
3551 
3552  // Get Pattern
3553  const casadi_int* colind = this->colind();
3554  const casadi_int* row = this->row();
3555  const casadi_int* y_colind = y.colind();
3556  const casadi_int* y_row = y.row();
3557 
3558  // If same number of rows, check if patterns are identical
3559  if (size1()==y.size1()) return is_equal(y.size1(), y.size2(), y_colind, y_row);
3560 
3561  // Loop over the elements
3562  for (casadi_int cc=0; cc<size2(); ++cc) {
3563  for (casadi_int el=colind[cc]; el<colind[cc+1]; ++el) {
3564  casadi_int rr=row[el];
3565 
3566  // Get row and column of y
3567  casadi_int loc = rr+size1()*cc;
3568  casadi_int rr_y = loc % y.size1();
3569  casadi_int cc_y = loc / y.size1();
3570 
3571  // Make sure matching
3572  if (rr_y != y_row[el] || el<y_colind[cc_y] || el>=y_colind[cc_y+1])
3573  return false;
3574  }
3575  }
3576 
3577  // Reshape if reached this point
3578  return true;
3579  }

References colind(), is_dense(), is_equal(), nnz(), numel(), row(), size1(), and size2().

◆ is_row()

bool casadi::SparsityInternal::is_row ( ) const

Extra doc: https://github.com/casadi/casadi/wiki/L_fh

Definition at line 1789 of file sparsity_internal.cpp.

1789  {
1790  return size1()==1;
1791  }

References size1().

Referenced by is_vector().

◆ is_scalar()

bool casadi::SparsityInternal::is_scalar ( bool  scalar_and_dense) const

Definition at line 1781 of file sparsity_internal.cpp.

1781  {
1782  return size2()==1 && size1()==1 && (!scalar_and_dense || nnz()==1);
1783  }

References nnz(), size1(), and size2().

Referenced by _mtimes().

◆ is_selection()

bool casadi::SparsityInternal::is_selection ( bool  allow_empty = false) const

Extra doc: https://github.com/casadi/casadi/wiki/L_25w

Definition at line 1879 of file sparsity_internal.cpp.

1879  {
1880  return is_orthonormal_rows(allow_empty);
1881  }
casadi::SparsityInternal::is_orthonormal_rows
bool is_orthonormal_rows(bool allow_empty=false) const
Are the rows of the pattern orthonormal ?
Definition: sparsity_internal.cpp:1855

References is_orthonormal_rows().

◆ is_square()

bool casadi::SparsityInternal::is_square ( ) const

Definition at line 1831 of file sparsity_internal.cpp.

1831  {
1832  return size2() == size1();
1833  }

References size1(), and size2().

Referenced by is_orthonormal(), star_coloring(), and star_coloring2().

◆ is_stacked()

bool casadi::SparsityInternal::is_stacked ( const Sparsity &  y,
casadi_int  n 
) const

Definition at line 2438 of file sparsity_internal.cpp.

2438  {
2439  // Quick true if the objects are equal
2440  if (n==1 && is_equal(y)) return true;
2441  // Get sparsity patterns
2442  casadi_int size1 = this->size1();
2443  casadi_int size2 = this->size2();
2444  const casadi_int* colind = this->colind();
2445  const casadi_int* row = this->row();
2446  casadi_int y_size1 = y.size1();
2447  casadi_int y_size2 = y.size2();
2448  const casadi_int* y_colind = y.colind();
2449  const casadi_int* y_row = y.row();
2450  // Make sure dimensions are consistent
2451  if (size1!=y_size1 || size2!=n*y_size2) return false;
2452  // Make sure number of nonzeros are consistent
2453  casadi_int nnz = colind[size2], y_nnz = y_colind[y_size2];
2454  if (nnz!=n*y_nnz) return false;
2455  // Quick return if dense
2456  if (y_nnz==y_size1*y_size2) return true;
2457  // Offset
2458  casadi_int offset = 0;
2459  // Skip the initial zero
2460  colind++;
2461  // For all repeats
2462  for (int i=0; i<n; ++i) {
2463  // Compare column offsets
2464  for (int c=0; c<y_size2; ++c) if (y_colind[c+1]+offset != *colind++) return false;
2465  // Compare row indices
2466  for (int k=0; k<y_nnz; ++k) if (y_row[k] != *row++) return false;
2467  // Update nonzero offset
2468  offset += y_nnz;
2469  }
2470  // Equal if reached this point
2471  return true;
2472  }

References casadi::Sparsity::colind(), colind(), is_equal(), nnz(), casadi::Sparsity::row(), row(), casadi::Sparsity::size1(), size1(), casadi::Sparsity::size2(), and size2().

◆ is_subset()

bool casadi::SparsityInternal::is_subset ( const Sparsity &  rhs) const

Definition at line 2332 of file sparsity_internal.cpp.

2332  {
2333  if (is_equal(rhs)) return true;
2334  std::vector<unsigned char> mapping;
2335  shared_from_this<Sparsity>().unite(rhs, mapping);
2336  for (auto e : mapping) {
2337  if (e==1) return false;
2338  }
2339  return true;
2340  }

References is_equal().

◆ is_symmetric()

bool casadi::SparsityInternal::is_symmetric ( ) const

Definition at line 1835 of file sparsity_internal.cpp.

1835  {
1836  return is_transpose(*this);
1837  }
casadi::SparsityInternal::is_transpose
bool is_transpose(const SparsityInternal &y) const
Check if the sparsity is the transpose of another.
Definition: sparsity_internal.cpp:3496

References is_transpose().

Referenced by amd().

◆ is_transpose()

bool casadi::SparsityInternal::is_transpose ( const SparsityInternal &  y) const

Definition at line 3496 of file sparsity_internal.cpp.

3496  {
3497  // Assert dimensions and number of nonzeros
3498  if (size2()!=y.size1() || size1()!=y.size2() || nnz()!=y.nnz())
3499  return false;
3500 
3501  // Quick return if empty interior or dense
3502  if (nnz()==0 || is_dense())
3503  return true;
3504 
3505  // Run algorithm on the pattern with the least number of rows
3506  if (size1()>size2()) return y.is_transpose(*this);
3507 
3508  // Index counter for columns of the possible transpose
3509  std::vector<casadi_int> y_col_count(y.size2(), 0);
3510  const casadi_int* colind = this->colind();
3511  const casadi_int* row = this->row();
3512  const casadi_int* y_colind = y.colind();
3513  const casadi_int* y_row = y.row();
3514 
3515  // Loop over the columns
3516  for (casadi_int i=0; i<size2(); ++i) {
3517 
3518  // Loop over the nonzeros
3519  for (casadi_int el=colind[i]; el<colind[i+1]; ++el) {
3520 
3521  // Get the row
3522  casadi_int j=row[el];
3523 
3524  // Get the element of the possible transpose
3525  casadi_int el_y = y_colind[j] + y_col_count[j]++;
3526 
3527  // Quick return if element doesn't exist
3528  if (el_y>=y_colind[j+1]) return false;
3529 
3530  // Get the row of the possible transpose
3531  casadi_int j_y = y_row[el_y];
3532 
3533  // Quick return if mismatch
3534  if (j_y != i) return false;
3535  }
3536  }
3537 
3538  // Transpose if reached this point
3539  return true;
3540  }

References colind(), is_dense(), is_transpose(), nnz(), row(), size1(), and size2().

Referenced by is_symmetric(), and is_transpose().

◆ is_tril()

bool casadi::SparsityInternal::is_tril ( bool  strictly) const

Definition at line 3715 of file sparsity_internal.cpp.

3715  {
3716  const casadi_int* colind = this->colind();
3717  const casadi_int* row = this->row();
3718  // Loop over columns
3719  for (casadi_int i=0; i<size2(); ++i) {
3720  if (colind[i] != colind[i+1]) { // if there are any elements of the column
3721  // check row of the top-most element of the column
3722  casadi_int rr = row[colind[i]];
3723  // Not lower triangular if row > i
3724  if (strictly ? rr <= i : rr < i) return false;
3725  }
3726  }
3727  // all columns ok
3728  return true;
3729  }

References colind(), row(), and size2().

◆ is_triu()

bool casadi::SparsityInternal::is_triu ( bool  strictly) const

Definition at line 3731 of file sparsity_internal.cpp.

3731  {
3732  const casadi_int* colind = this->colind();
3733  const casadi_int* row = this->row();
3734  // Loop over columns
3735  for (casadi_int i=0; i<size2(); ++i) {
3736  if (colind[i] != colind[i+1]) { // if there are any elements of the column
3737  // check row of the bottom-most element of the column
3738  casadi_int rr = row[colind[i+1]-1];
3739  // Not upper triangular if row>i
3740  if (strictly ? rr >= i : rr > i) return false;
3741  }
3742  }
3743  // all columns ok
3744  return true;
3745  }

References colind(), row(), and size2().

◆ is_vector()

bool casadi::SparsityInternal::is_vector ( ) const

Extra doc: https://github.com/casadi/casadi/wiki/L_fj

Definition at line 1797 of file sparsity_internal.cpp.

1797  {
1798  return is_row() || is_column();
1799  }
casadi::SparsityInternal::is_column
bool is_column() const
Check if the pattern is a column vector (i.e. size2()==1)
Definition: sparsity_internal.cpp:1793
casadi::SparsityInternal::is_row
bool is_row() const
Check if the pattern is a row vector (i.e. size1()==1)
Definition: sparsity_internal.cpp:1789

References is_column(), and is_row().

Referenced by get_diag().

◆ largest_first()

std::vector< casadi_int > casadi::SparsityInternal::largest_first ( ) const

Definition at line 3406 of file sparsity_internal.cpp.

3406  {
3407  std::vector<casadi_int> degree = get_colind();
3408  casadi_int max_degree = 0;
3409  for (casadi_int k=0; k<size2(); ++k) {
3410  degree[k] = degree[k+1]-degree[k];
3411  max_degree = std::max(max_degree, 1+degree[k]);
3412  }
3413  degree.resize(size2());
3414 
3415  // Vector for binary sort
3416  std::vector<casadi_int> degree_count(max_degree+1, 0);
3417  for (std::vector<casadi_int>::const_iterator it=degree.begin(); it!=degree.end(); ++it) {
3418  degree_count.at(*it+1)++;
3419  }
3420 
3421  // Cumsum to get the offset for each degree
3422  for (casadi_int d=0; d<max_degree; ++d) {
3423  degree_count[d+1] += degree_count[d];
3424  }
3425 
3426  // Now a bucket sort
3427  std::vector<casadi_int> ordering(size2());
3428  for (casadi_int k=size2()-1; k>=0; --k) {
3429  ordering[degree_count[degree[k]]++] = k;
3430  }
3431 
3432  // Invert the ordering
3433  std::vector<casadi_int>& reverse_ordering = degree_count; // reuse memory
3434  reverse_ordering.resize(ordering.size());
3435  std::copy(ordering.begin(), ordering.end(), reverse_ordering.rbegin());
3436 
3437  // Return the ordering
3438  return reverse_ordering;
3439  }

References get_colind(), and size2().

Referenced by star_coloring(), and star_coloring2().

◆ ldl_colind()

void casadi::SparsityInternal::ldl_colind ( const casadi_int *  sp,
casadi_int *  parent,
casadi_int *  l_colind,
casadi_int *  w 
)
static

Strictly lower entries only. Ref: User Guide for LDL by Tim Davis len[colind] = ncol+1 len[parent] = ncol len[w] >= ncol Modified version of LDL Copyright(c) Timothy A. Davis, 2005-2013 Licensed as a derivative work under the GNU LGPL

Extra doc: https://github.com/casadi/casadi/wiki/L_ew

Definition at line 377 of file sparsity_internal.cpp.

378  {
379  /* Modified version of LDL
380  Copyright(c) Timothy A. Davis, 2005-2013
381  Licensed as a derivative work under the GNU LGPL
382  */
383  casadi_int n = sp[0];
384  const casadi_int *colind=sp+2, *row=sp+2+n+1;
385  // Local variables
386  casadi_int r, c, k;
387  // Work vectors
388  casadi_int* visited=w; w+=n;
389  // Loop over columns
390  for (c=0; c<n; ++c) {
391  // L(c,:) pattern: all nodes reachable in etree from nz in A(0:c-1,c)
392  parent[c] = -1; // parent of c is not yet known
393  visited[c] = c; // mark node c as visited
394  l_colind[1+c] = 0; // count of nonzeros in column c of L
395  // Loop over strictly upper triangular entries A
396  for (k=colind[c]; k<colind[c+1] && (r=row[k])<c; ++k) {
397  // Follow path from r to root of etree, stop at visited node
398  while (visited[r]!=c) {
399  // Find parent of r if not yet determined
400  if (parent[r]==-1) parent[r]=c;
401  l_colind[1+r]++; // L(c,r) is nonzero
402  visited[r] = c; // mark r as visited
403  r=parent[r]; // proceed to parent row
404  }
405  }
406  }
407  // Cumsum
408  l_colind[0] = 0;
409  for (c=0; c<n; ++c) l_colind[c+1] += l_colind[c];
410  }

References colind(), row(), and sp().

Referenced by casadi::Sparsity::ldl().

◆ ldl_row()

void casadi::SparsityInternal::ldl_row ( const casadi_int *  sp,
const casadi_int *  parent,
casadi_int *  l_colind,
casadi_int *  l_row,
casadi_int *  w 
)
static

Strictly lower entries only. Ref: User Guide for LDL by Tim Davis len[w] >= n Modified version of LDL Copyright(c) Timothy A. Davis, 2005-2013 Licensed as a derivative work under the GNU LGPL

Extra doc: https://github.com/casadi/casadi/wiki/L_ex

Definition at line 412 of file sparsity_internal.cpp.

414  {
415  /* Modified version of LDL
416  Copyright(c) Timothy A. Davis, 2005-2013
417  Licensed as a derivative work under the GNU LGPL
418  */
419  // Extract sparsity
420  casadi_int n = sp[0];
421  const casadi_int *colind = sp+2, *row = sp+n+3;
422  // Work vectors
423  casadi_int *visited=w; w+=n;
424  // Local variables
425  casadi_int r, c, k;
426  // Compute nonzero pattern of kth row of L
427  for (c=0; c<n; ++c) {
428  // Not yet visited
429  visited[c] = c;
430  // Loop over nonzeros in upper triangular half
431  for (k=colind[c]; k<colind[c+1] && (r=row[k])<c; ++k) {
432  // Loop over dependent rows
433  while (visited[r]!=c) {
434  l_row[l_colind[r]++] = c; // L(c,r) is nonzero
435  visited[r] = c; // mark r as visited
436  r=parent[r]; // proceed to parent row
437  }
438  }
439  }
440  // Restore l_colind by shifting it forward
441  k=0;
442  for (c=0; c<n; ++c) {
443  r=l_colind[c];
444  l_colind[c]=k;
445  k=r;
446  }
447  }

References colind(), row(), and sp().

Referenced by casadi::Sparsity::ldl().

◆ leaf()

casadi_int casadi::SparsityInternal::leaf ( casadi_int  i,
casadi_int  j,
const casadi_int *  first,
casadi_int *  maxfirst,
casadi_int *  prevleaf,
casadi_int *  ancestor,
casadi_int *  jleaf 
)
static

Ref: Chapter 4, Direct Methods for Sparse Linear Systems by Tim Davis Modified version of cs_leaf in CSparse Copyright(c) Timothy A. Davis, 2006-2009 Licensed as a derivative work under the GNU LGPL

Extra doc: https://github.com/casadi/casadi/wiki/L_er

Definition at line 127 of file sparsity_internal.cpp.

129  {
130  /* Modified version of cs_leaf in CSparse
131  Copyright(c) Timothy A. Davis, 2006-2009
132  Licensed as a derivative work under the GNU LGPL
133  */
134  casadi_int q, s, sparent, jprev;
135  *jleaf = 0;
136  // Quick return if j is not a leaf
137  if (i<=j || first[j]<=maxfirst[i]) return -1;
138  // Update max first[j] seen so far
139  maxfirst[i] = first[j];
140  // Previous leaf of ith subtree
141  jprev = prevleaf[i];
142  prevleaf[i] = j;
143  // j is first or subsequent leaf
144  *jleaf = (jprev == -1) ? 1 : 2;
145  // if first leaf, q is root of ith subtree
146  if (*jleaf==1) return i;
147  // Path compression
148  for (q=jprev; q!=ancestor[q]; q=ancestor[q]) {}
149  for (s=jprev; s!=q; s=sparent) {
150  sparent = ancestor[s];
151  ancestor[s] = q;
152  }
153  // Return least common ancestor
154  return q;
155  }

Referenced by qr_counts().

◆ makeDense()

Sparsity casadi::SparsityInternal::makeDense ( std::vector< casadi_int > &  mapping) const

Definition at line 2662 of file sparsity_internal.cpp.

2662  {
2663  const casadi_int* colind = this->colind();
2664  const casadi_int* row = this->row();
2665  mapping.resize(nnz());
2666  for (casadi_int i=0; i<size2(); ++i) {
2667  for (casadi_int el=colind[i]; el<colind[i+1]; ++el) {
2668  casadi_int j = row[el];
2669  mapping[el] = j + i*size1();
2670  }
2671  }
2672 
2673  return Sparsity::dense(size1(), size2());
2674  }

References colind(), casadi::Sparsity::dense(), nnz(), row(), size1(), and size2().

Referenced by casadi::Sparsity::makeDense().

◆ matched()

void casadi::SparsityInternal::matched ( casadi_int  n,
const std::vector< casadi_int > &  wj,
const std::vector< casadi_int > &  imatch,
std::vector< casadi_int > &  p,
std::vector< casadi_int > &  q,
std::vector< casadi_int > &  cc,
std::vector< casadi_int > &  rr,
casadi_int  set,
casadi_int  mark 
)
static

The implementation is a modified version of cs_matched in CSparse Copyright(c) Timothy A. Davis, 2006-2009 Licensed as a derivative work under the GNU LGPL

Extra doc: https://github.com/casadi/casadi/wiki/L_f0

Definition at line 1032 of file sparsity_internal.cpp.

1035  {
1036  /*
1037  Modified version of cs_matched in CSparse
1038  Copyright(c) Timothy A. Davis, 2006-2009
1039  Licensed as a derivative work under the GNU LGPL
1040  */
1041  casadi_int kc = cc[set];
1042  casadi_int kr = rr[set-1] ;
1043  for (casadi_int j=0; j<n; ++j) {
1044  // skip if j is not in C set
1045  if (wj[j] != mark) continue;
1046 
1047  p[kr++] = imatch[j] ;
1048  q[kc++] = j ;
1049  }
1050 
1051  cc[set+1] = kc ;
1052  rr[set] = kr ;
1053  }

Referenced by dmperm().

◆ maxtrans()

void casadi::SparsityInternal::maxtrans ( std::vector< casadi_int > &  imatch,
std::vector< casadi_int > &  jmatch,
Sparsity &  trans,
casadi_int  seed 
) const

The implementation is a modified version of cs_maxtrans in CSparse Copyright(c) Timothy A. Davis, 2006-2009 Licensed as a derivative work under the GNU LGPL

Extra doc: https://github.com/casadi/casadi/wiki/L_f5

Definition at line 1154 of file sparsity_internal.cpp.

1155  {
1156  /*
1157  Modified version of cs_maxtrans in CSparse
1158  Copyright(c) Timothy A. Davis, 2006-2009
1159  Licensed as a derivative work under the GNU LGPL
1160  */
1161  const casadi_int* colind = this->colind();
1162  const casadi_int* row = this->row();
1163 
1164  casadi_int n2 = 0, m2 = 0;
1165 
1166  // allocate result
1167  jmatch.resize(size1());
1168  imatch.resize(size2());
1169  std::vector<casadi_int> w(size1()+size2());
1170 
1171  // count nonempty columns and rows
1172  casadi_int k=0;
1173  for (casadi_int j=0; j<size2(); ++j) {
1174  n2 += (colind[j] < colind[j+1]);
1175  for (casadi_int p=colind[j]; p < colind[j+1]; ++p) {
1176  w[row[p]] = 1;
1177 
1178  // count entries already on diagonal
1179  k += (j == row[p]);
1180  }
1181  }
1182 
1183  // quick return if diagonal zero-free
1184  if (k == std::min(size1(), size2())) {
1185  casadi_int i;
1186  for (i=0; i<k; ++i) jmatch[i] = i;
1187  for (; i<size1(); ++i) jmatch[i] = -1;
1188 
1189  casadi_int j;
1190  for (j=0; j<k; ++j) imatch[j] = j;
1191  for (; j<size2(); ++j) imatch[j] = -1;
1192  }
1193 
1194  for (casadi_int i=0; i<size1(); ++i) m2 += w[i];
1195 
1196  // transpose if needed
1197  if (m2 < n2 && trans.is_null())
1198  trans = T();
1199 
1200  // Get pointer to sparsity
1201  const SparsityInternal* C = m2 < n2 ? static_cast<const SparsityInternal*>(trans.get()) : this;
1202  const casadi_int* C_colind = C->colind();
1203 
1204  std::vector<casadi_int>& Cjmatch = m2 < n2 ? imatch : jmatch;
1205  std::vector<casadi_int>& Cimatch = m2 < n2 ? jmatch : imatch;
1206 
1207  // get workspace
1208  w.resize(5 * C->size2());
1209 
1210  casadi_int *cheap = &w.front() + C->size2();
1211  casadi_int *js = &w.front() + 2*C->size2();
1212  casadi_int *is = &w.front() + 3*C->size2();
1213  casadi_int *ps = &w.front() + 4*C->size2();
1214 
1215  // for cheap assignment
1216  for (casadi_int j=0; j<C->size2(); ++j)
1217  cheap[j] = C_colind[j];
1218 
1219  // all rows unflagged
1220  for (casadi_int j=0; j<C->size2(); ++j)
1221  w[j] = -1;
1222 
1223  // nothing matched yet
1224  for (casadi_int i=0; i<C->size1(); ++i)
1225  Cjmatch[i] = -1;
1226 
1227  // q = random permutation
1228  std::vector<casadi_int> q = randperm(C->size2(), seed);
1229 
1230  // augment, starting at row q[k]
1231  for (k=0; k<C->size2(); ++k) {
1232  C->augment(!q.empty() ? q[k]: k, Cjmatch, cheap, w, js, is, ps);
1233  }
1234 
1235  // find col match
1236  for (casadi_int j=0; j<C->size2(); ++j)
1237  Cimatch[j] = -1;
1238 
1239  for (casadi_int i = 0; i<C->size1(); ++i)
1240  if (Cjmatch[i] >= 0)
1241  Cimatch[Cjmatch[i]] = i;
1242  }
casadi::SparsityInternal::SparsityInternal
SparsityInternal(casadi_int nrow, casadi_int ncol, const casadi_int *colind, const casadi_int *row)
Construct a sparsity pattern from arrays.
Definition: sparsity_internal.cpp:450
casadi::SparsityInternal::randperm
static std::vector< casadi_int > randperm(casadi_int n, casadi_int seed)
return a random permutation vector
Definition: sparsity_internal.cpp:1408

References casadi::C, colind(), casadi::Sparsity::get(), casadi::GenericShared< Shared, Internal >::is_null(), randperm(), row(), size1(), size2(), and T().

Referenced by dmperm().

◆ multiply()

Sparsity casadi::SparsityInternal::multiply ( const Sparsity &  B) const

The implementation is a modified version of cs_multiply in CSparse Copyright(c) Timothy A. Davis, 2006-2009 Licensed as a derivative work under the GNU LGPL

Extra doc: https://github.com/casadi/casadi/wiki/L_fc

Definition at line 1575 of file sparsity_internal.cpp.

1575  {
1576  /*
1577  Modified version of cs_multiply in CSparse
1578  Copyright(c) Timothy A. Davis, 2006-2009
1579  Licensed as a derivative work under the GNU LGPL
1580  */
1581  casadi_int nz = 0;
1582  casadi_assert(size2() == B.size1(), "Dimension mismatch.");
1583  casadi_int m = size1();
1584  casadi_int anz = nnz();
1585  casadi_int n = B.size2();
1586  const casadi_int* Bp = B.colind();
1587  const casadi_int* Bi = B.row();
1588  casadi_int bnz = Bp[n];
1589 
1590  // get workspace
1591  std::vector<casadi_int> w(m);
1592 
1593  // allocate result
1594  std::vector<casadi_int> C_colind(n+1, 0), C_row;
1595 
1596  C_colind.resize(anz + bnz);
1597 
1598  casadi_int* Cp = &C_colind.front();
1599  for (casadi_int j=0; j<n; ++j) {
1600  if (nz+m > C_row.size()) {
1601  C_row.resize(2*(C_row.size())+m);
1602  }
1603 
1604  // row j of C starts here
1605  Cp[j] = nz;
1606  for (casadi_int p = Bp[j] ; p<Bp[j+1] ; ++p) {
1607  nz = scatter(Bi[p], w, j+1, &C_row.front(), nz);
1608  }
1609  }
1610 
1611  // finalize the last row of C
1612  Cp[n] = nz;
1613  C_row.resize(nz);
1614 
1615  // Success
1616  return Sparsity(m, n, C_colind, C_row);
1617  }
casadi::SparsityInternal::scatter
casadi_int scatter(casadi_int j, std::vector< casadi_int > &w, casadi_int mark, casadi_int *Ci, casadi_int nz) const
x = x + beta * A(:, j), where x is a dense vector and A(:, j) is sparse
Definition: sparsity_internal.cpp:1549

References casadi::Sparsity::colind(), nnz(), casadi::Sparsity::row(), scatter(), casadi::Sparsity::size1(), size1(), casadi::Sparsity::size2(), and size2().

◆ nnz()

casadi_int casadi::SparsityInternal::nnz ( ) const
inline

Definition at line 113 of file sparsity_internal.hpp.

113 { return colind()[size2()];}

Referenced by _appendColumns(), _appendVector(), _enlargeColumns(), _enlargeRows(), _erase(), _removeDuplicates(), _reshape(), amd(), dim(), drop_diag(), export_code(), find(), get_col(), get_nz(), get_row(), is_dense(), is_diag(), is_equal(), is_orthonormal(), is_orthonormal_columns(), is_orthonormal_rows(), is_reshape(), is_scalar(), is_stacked(), is_transpose(), makeDense(), multiply(), nnz_diag(), nnz_lower(), nnz_upper(), permute(), pmult(), repr_el(), star_coloring2(), and sub().

◆ nnz_diag()

casadi_int casadi::SparsityInternal::nnz_diag ( ) const

Definition at line 1895 of file sparsity_internal.cpp.

1895  {
1896  const casadi_int* colind = this->colind();
1897  const casadi_int* row = this->row();
1898  casadi_int nnz = 0;
1899  for (casadi_int cc=0; cc<size2(); ++cc) {
1900  for (casadi_int el = colind[cc]; el < colind[cc+1]; ++el) {
1901  nnz += row[el]==cc;
1902  }
1903  }
1904  return nnz;
1905  }

References colind(), nnz(), row(), and size2().

◆ nnz_lower()

casadi_int casadi::SparsityInternal::nnz_lower ( bool  strictly = false) const

Definition at line 1883 of file sparsity_internal.cpp.

1883  {
1884  const casadi_int* colind = this->colind();
1885  const casadi_int* row = this->row();
1886  casadi_int nnz = 0;
1887  for (casadi_int cc=0; cc<size2(); ++cc) {
1888  for (casadi_int el = colind[cc]; el<colind[cc+1]; ++el) {
1889  if (cc<row[el] || (!strictly && cc==row[el])) nnz++;
1890  }
1891  }
1892  return nnz;
1893  }

References colind(), nnz(), row(), and size2().

◆ nnz_upper()

casadi_int casadi::SparsityInternal::nnz_upper ( bool  strictly = false) const

Definition at line 1907 of file sparsity_internal.cpp.

1907  {
1908  const casadi_int* colind = this->colind();
1909  const casadi_int* row = this->row();
1910  casadi_int nnz = 0;
1911  for (casadi_int cc=0; cc<size2(); ++cc) {
1912  for (casadi_int el = colind[cc]; el<colind[cc+1]; ++el) {
1913  if (cc>row[el] || (!strictly && cc==row[el])) nnz++;
1914  }
1915  }
1916  return nnz;
1917  }

References colind(), nnz(), row(), and size2().

Referenced by star_coloring2().

◆ numel()

casadi_int casadi::SparsityInternal::numel ( ) const

Definition at line 477 of file sparsity_internal.cpp.

477  {
478  return size1()*size2();
479  }

References size1(), and size2().

Referenced by _erase(), _reshape(), find(), is_dense(), is_equal(), is_reshape(), repr_el(), and sub().

◆ pattern_inverse()

Sparsity casadi::SparsityInternal::pattern_inverse ( ) const

Definition at line 2482 of file sparsity_internal.cpp.

2482  {
2483  // Quick return clauses
2484  if (is_empty()) return Sparsity::dense(size1(), size2());
2485  if (is_dense()) return Sparsity(size1(), size2());
2486 
2487  // Sparsity of the result
2488  std::vector<casadi_int> row_ret;
2489  std::vector<casadi_int> colind_ret=get_colind();
2490  const casadi_int* colind = this->colind();
2491  const casadi_int* row = this->row();
2492 
2493  // Loop over columns
2494  for (casadi_int i=0;i<size2();++i) {
2495  // Update colind vector of the result
2496  colind_ret[i+1]=colind_ret[i]+size1()-(colind[i+1]-colind[i]);
2497 
2498  // Counter of new row indices
2499  casadi_int j=0;
2500 
2501  // Loop over all nonzeros
2502  for (casadi_int k=colind[i];k<colind[i+1];++k) {
2503 
2504  // Try to reach current nonzero
2505  while (j<row[k]) {
2506  // And meanwhile, add nonzeros to the result
2507  row_ret.push_back(j);
2508  j++;
2509  }
2510  j++;
2511  }
2512  // Process the remainder up to the row size
2513  while (j < size1()) {
2514  row_ret.push_back(j);
2515  j++;
2516  }
2517  }
2518 
2519  // Return result
2520  return Sparsity(size1(), size2(), colind_ret, row_ret);
2521  }

References colind(), casadi::Sparsity::dense(), get_colind(), is_dense(), is_empty(), row(), size1(), and size2().

Referenced by casadi::Sparsity::pattern_inverse().

◆ permute() [1/2]

Sparsity casadi::SparsityInternal::permute ( const std::vector< casadi_int > &  pinv,
const std::vector< casadi_int > &  q,
casadi_int  values 
) const

Definition at line 1456 of file sparsity_internal.cpp.

1457  {
1458  std::vector<casadi_int> colind_C, row_C;
1459  permute(pinv, q, values, colind_C, row_C);
1460  return Sparsity(size1(), size2(), colind_C, row_C);
1461  }

References size1(), and size2().

Referenced by dmperm().

◆ permute() [2/2]

void casadi::SparsityInternal::permute ( const std::vector< casadi_int > &  pinv,
const std::vector< casadi_int > &  q,
casadi_int  values,
std::vector< casadi_int > &  colind_C,
std::vector< casadi_int > &  row_C 
) const

The implementation is a modified version of cs_permute in CSparse Copyright(c) Timothy A. Davis, 2006-2009 Licensed as a derivative work under the GNU LGPL

Extra doc: https://github.com/casadi/casadi/wiki/L_f9

Definition at line 1463 of file sparsity_internal.cpp.

1466  {
1467  /*
1468  Modified version of cs_permute in CSparse
1469  Copyright(c) Timothy A. Davis, 2006-2009
1470  Licensed as a derivative work under the GNU LGPL
1471  */
1472  const casadi_int* colind = this->colind();
1473  const casadi_int* row = this->row();
1474 
1475  // alloc column offsets
1476  colind_C.resize(size2()+1);
1477 
1478  // Row for each nonzero
1479  row_C.resize(nnz());
1480  casadi_int nz = 0;
1481  for (casadi_int k = 0; k<size2(); ++k) {
1482  // row k of C is row q[k] of A
1483  colind_C[k] = nz;
1484 
1485  casadi_int j = !q.empty() ? (q[k]) : k;
1486 
1487  for (casadi_int t = colind[j]; t<colind[j+1]; ++t) {
1488  row_C[nz++] = !pinv.empty() ? (pinv[row[t]]) : row[t] ;
1489  }
1490  }
1491 
1492  // finalize the last row of C
1493  colind_C[size2()] = nz;
1494  }

References colind(), nnz(), row(), and size2().

◆ pmult()

Sparsity casadi::SparsityInternal::pmult ( const std::vector< casadi_int > &  p,
bool  permute_rows = true,
bool  permute_cols = true,
bool  invert_permutation = false 
) const

Definition at line 3441 of file sparsity_internal.cpp.

3442  {
3443  // Invert p, possibly
3444  std::vector<casadi_int> p_inv;
3445  if (invert_permutation) {
3446  p_inv.resize(p.size());
3447  for (casadi_int k=0; k<p.size(); ++k) {
3448  p_inv[p[k]] = k;
3449  }
3450  }
3451  const std::vector<casadi_int>& pp = invert_permutation ? p_inv : p;
3452 
3453  // Get columns
3454  std::vector<casadi_int> col = get_col();
3455 
3456  // Get rows
3457  const casadi_int* row = this->row();
3458 
3459  // Sparsity of the return matrix
3460  std::vector<casadi_int> new_row(col.size()), new_col(col.size());
3461 
3462  // Possibly permute columns
3463  if (permute_columns) {
3464  // Assert dimensions
3465  casadi_assert_dev(p.size()==size2());
3466 
3467  // Permute
3468  for (casadi_int k=0; k<col.size(); ++k) {
3469  new_col[k] = pp[col[k]];
3470  }
3471 
3472  } else {
3473  // No permutation of columns
3474  std::copy(col.begin(), col.end(), new_col.begin());
3475  }
3476 
3477  // Possibly permute rows
3478  if (permute_rows) {
3479  // Assert dimensions
3480  casadi_assert_dev(p.size()==size1());
3481 
3482  // Permute
3483  for (casadi_int k=0; k<nnz(); ++k) {
3484  new_row[k] = pp[row[k]];
3485  }
3486 
3487  } else {
3488  // No permutation of rows
3489  std::copy(row, row+nnz(), new_row.begin());
3490  }
3491 
3492  // Return permuted matrix
3493  return Sparsity::triplet(size1(), size2(), new_row, new_col);
3494  }
casadi::SparsityInternal::get_col
std::vector< casadi_int > get_col() const
Get the column for each nonzero.
Definition: sparsity_internal.cpp:490
casadi::invert_permutation
std::vector< casadi_int > invert_permutation(const std::vector< casadi_int > &a)
inverse a permutation vector
Definition: casadi_misc.cpp:178

References get_col(), casadi::invert_permutation(), nnz(), casadi::Sparsity::resize(), row(), size1(), size2(), and casadi::Sparsity::triplet().

Referenced by casadi::Sparsity::pmult(), star_coloring(), and star_coloring2().

◆ postorder()

void casadi::SparsityInternal::postorder ( const casadi_int *  parent,
casadi_int  n,
casadi_int *  post,
casadi_int *  w 
)
static

Ref: Chapter 4, Direct Methods for Sparse Linear Systems by Tim Davis len[w] >= 3*n len[post] == n Modified version of cs_post in CSparse Copyright(c) Timothy A. Davis, 2006-2009 Licensed as a derivative work under the GNU LGPL

Extra doc: https://github.com/casadi/casadi/wiki/L_eq

Definition at line 99 of file sparsity_internal.cpp.

100  {
101  /* Modified version of cs_post in CSparse
102  Copyright(c) Timothy A. Davis, 2006-2009
103  Licensed as a derivative work under the GNU LGPL
104  */
105  casadi_int j, k=0;
106  // Work vectors
107  casadi_int *head, *next, *stack;
108  head=w; w+=n;
109  next=w; w+=n;
110  stack=w; w+=n;
111  // Empty linked lists
112  for (j=0; j<n; ++j) head[j] = -1;
113  // Traverse nodes in reverse order
114  for (j=n-1; j>=0; --j) {
115  if (parent[j]!=-1) {
116  next[j] = head[parent[j]];
117  head[parent[j]] = j;
118  }
119  }
120  for (j=0; j<n; j++) {
121  if (parent[j]==-1) {
122  k = postorder_dfs(j, k, head, next, post, stack);
123  }
124  }
125  }

References postorder_dfs().

Referenced by qr_init().

◆ postorder_dfs()

casadi_int casadi::SparsityInternal::postorder_dfs ( casadi_int  j,
casadi_int  k,
casadi_int *  head,
const casadi_int *  next,
casadi_int *  post,
casadi_int *  stack 
)
static

Ref: Chapter 4, Direct Methods for Sparse Linear Systems by Tim Davis Modified version of cs_tdfs in CSparse Copyright(c) Timothy A. Davis, 2006-2009 Licensed as a derivative work under the GNU LGPL

Extra doc: https://github.com/casadi/casadi/wiki/L_ep

Definition at line 74 of file sparsity_internal.cpp.

76  {
77  /* Modified version of cs_tdfs in CSparse
78  Copyright(c) Timothy A. Davis, 2006-2009
79  Licensed as a derivative work under the GNU LGPL
80  */
81  casadi_int i, p, top=0;
82  stack[0] = j;
83  while (top>=0) {
84  p = stack[top];
85  i = head[p];
86  if (i==-1) {
87  // No children
88  top--;
89  post[k++] = p;
90  } else {
91  // Add to stack
92  head[p] = next[i];
93  stack[++top] = i;
94  }
95  }
96  return k;
97  }

Referenced by amd(), and postorder().

◆ qr_counts()

casadi_int casadi::SparsityInternal::qr_counts ( const casadi_int *  tr_sp,
const casadi_int *  parent,
const casadi_int *  post,
casadi_int *  counts,
casadi_int *  w 
)
static

Ref: Chapter 4, Direct Methods for Sparse Linear Systems by Tim Davis len[counts] = ncol len[w] >= 5*ncol + nrow + 1 Modified version of cs_counts in CSparse Copyright(c) Timothy A. Davis, 2006-2009 Licensed as a derivative work under the GNU LGPL

Extra doc: https://github.com/casadi/casadi/wiki/L_es

Definition at line 157 of file sparsity_internal.cpp.

159  {
160  /* Modified version of cs_counts in CSparse
161  Copyright(c) Timothy A. Davis, 2006-2009
162  Licensed as a derivative work under the GNU LGPL
163  */
164  casadi_int ncol = *tr_sp++, nrow = *tr_sp++;
165  const casadi_int *rowind=tr_sp, *col=tr_sp+nrow+1;
166  casadi_int i, j, k, J, p, q, jleaf, *maxfirst, *prevleaf,
167  *ancestor, *head=nullptr, *next=nullptr, *first;
168  // Work vectors
169  ancestor=w; w+=ncol;
170  maxfirst=w; w+=ncol;
171  prevleaf=w; w+=ncol;
172  first=w; w+=ncol;
173  head=w; w+=ncol+1;
174  next=w; w+=nrow;
175  // Find first [j]
176  for (k=0; k<ncol; ++k) first[k]=-1;
177  for (k=0; k<ncol; ++k) {
178  j=post[k];
179  // counts[j]=1 if j is a leaf
180  counts[j] = (first[j]==-1) ? 1 : 0;
181  for (; j!=-1 && first[j]==-1; j=parent[j]) first[j]=k;
182  }
183  // Invert post (use ancestor as work vector)
184  for (k=0; k<ncol; ++k) ancestor[post[k]] = k;
185  for (k=0; k<ncol+1; ++k) head[k]=-1;
186  for (i=0; i<nrow; ++i) {
187  for (k=ncol, p=rowind[i]; p<rowind[i+1]; ++p) {
188  k = std::min(k, ancestor[col[p]]);
189  }
190  // Place row i in linked list k
191  next[i] = head[k];
192  head[k] = i;
193  }
194 
195  // Clear workspace
196  for (k=0; k<ncol; ++k) maxfirst[k]=-1;
197  for (k=0; k<ncol; ++k) prevleaf[k]=-1;
198  // Each node in its own set
199  for (i=0; i<ncol; ++i) ancestor[i]=i;
200  for (k=0; k<ncol; ++k) {
201  // j is the kth node in the postordered etree
202  j=post[k];
203  if (parent[j]!=-1) counts[parent[j]]--; // j is not a root
204  J=head[k];
205  while (J!=-1) { // J=j for LL' = A case
206  for (p=rowind[J]; p<rowind[J+1]; ++p) {
207  i=col[p];
208  q = leaf(i, j, first, maxfirst, prevleaf, ancestor, &jleaf);
209  if (jleaf>=1) counts[j]++; // A(i,j) is in skeleton
210  if (jleaf==2) counts[q]--; // account for overlap in q
211  }
212  J = next[J];
213  }
214  if (parent[j]!=-1) ancestor[j]=parent[j];
215  }
216  // Sum up counts of each child
217  for (j=0; j<ncol; ++j) {
218  if (parent[j]!=-1) counts[parent[j]] += counts[j];
219  }
220 
221  // Sum of counts
222  casadi_int sum_counts = 0;
223  for (j=0; j<ncol; ++j) sum_counts += counts[j];
224  return sum_counts;
225  }
casadi::SparsityInternal::leaf
static casadi_int leaf(casadi_int i, casadi_int j, const casadi_int *first, casadi_int *maxfirst, casadi_int *prevleaf, casadi_int *ancestor, casadi_int *jleaf)
Needed by casadi_qr_colind.
Definition: sparsity_internal.cpp:128

References leaf().

Referenced by qr_init().

◆ qr_init()

void casadi::SparsityInternal::qr_init ( const casadi_int *  sp,
const casadi_int *  sp_tr,
casadi_int *  leftmost,
casadi_int *  parent,
casadi_int *  pinv,
casadi_int *  nrow_ext,
casadi_int *  v_nnz,
casadi_int *  r_nnz,
casadi_int *  w 
)
static

Ref: Chapter 5, Direct Methods for Sparse Linear Systems by Tim Davis len[w] >= nrow + 7*ncol + 1 len[pinv] == nrow + ncol len[leftmost] == nrow

Extra doc: https://github.com/casadi/casadi/wiki/L_eu

Definition at line 288 of file sparsity_internal.cpp.

291  {
292  // Extract sparsity
293  casadi_int ncol = sp[1];
294  // Calculate elimination tree for A'A
295  etree(sp, parent, w, 1); // len[w] >= nrow+ncol
296  // Calculate postorder
297  casadi_int* post = w; w += ncol;
298  postorder(parent, ncol, post, w); // len[w] >= 3*ncol
299  // Calculate nnz in R
300  *r_nnz = qr_counts(sp_tr, parent, post, w, w+ncol);
301  // Calculate nnz in V
302  *v_nnz = qr_nnz(sp, pinv, leftmost, parent, nrow_ext, w);
303  }
casadi::SparsityInternal::postorder
static void postorder(const casadi_int *parent, casadi_int n, casadi_int *post, casadi_int *w)
Calculate the postorder permuation.
Definition: sparsity_internal.cpp:99
casadi::SparsityInternal::qr_counts
static casadi_int qr_counts(const casadi_int *tr_sp, const casadi_int *parent, const casadi_int *post, casadi_int *counts, casadi_int *w)
Calculate the column offsets for the QR R matrix.
Definition: sparsity_internal.cpp:158
casadi::SparsityInternal::qr_nnz
static casadi_int qr_nnz(const casadi_int *sp, casadi_int *pinv, casadi_int *leftmost, const casadi_int *parent, casadi_int *nrow_ext, casadi_int *w)
Calculate the number of nonzeros in the QR V matrix.
Definition: sparsity_internal.cpp:228
casadi::SparsityInternal::etree
static void etree(const casadi_int *sp, casadi_int *parent, casadi_int *w, casadi_int ata)
Calculate the elimination tree for a matrix.
Definition: sparsity_internal.cpp:35

References etree(), postorder(), qr_counts(), qr_nnz(), and sp().

Referenced by casadi::Sparsity::qr_sparse().

◆ qr_nnz()

casadi_int casadi::SparsityInternal::qr_nnz ( const casadi_int *  sp,
casadi_int *  pinv,
casadi_int *  leftmost,
const casadi_int *  parent,
casadi_int *  nrow_ext,
casadi_int *  w 
)
static

Ref: Chapter 5, Direct Methods for Sparse Linear Systems by Tim Davis len[w] >= nrow + 3*ncol len[pinv] == nrow + ncol len[leftmost] == nrow Modified version of cs_sqr in CSparse Copyright(c) Timothy A. Davis, 2006-2009 Licensed as a derivative work under the GNU LGPL

Extra doc: https://github.com/casadi/casadi/wiki/L_et

Definition at line 227 of file sparsity_internal.cpp.

229  {
230  /* Modified version of cs_sqr in CSparse
231  Copyright(c) Timothy A. Davis, 2006-2009
232  Licensed as a derivative work under the GNU LGPL
233  */
234  // Extract sparsity
235  casadi_int nrow = sp[0], ncol = sp[1];
236  const casadi_int *colind=sp+2, *row=sp+2+ncol+1;
237  // Work vectors
238  casadi_int *next=w; w+=nrow;
239  casadi_int *head=w; w+=ncol;
240  casadi_int *tail=w; w+=ncol;
241  casadi_int *nque=w; w+=ncol;
242  // Local variables
243  casadi_int r, c, k, pa;
244  // Clear queue
245  for (c=0; c<ncol; ++c) head[c] = -1;
246  for (c=0; c<ncol; ++c) tail[c] = -1;
247  for (c=0; c<ncol; ++c) nque[c] = 0;
248  for (r=0; r<nrow; ++r) leftmost[r] = -1;
249  // leftmost[r] = min(find(A(r,:)))
250  for (c=ncol-1; c>=0; --c) {
251  for (k=colind[c]; k<colind[c+1]; ++k) {
252  leftmost[row[k]] = c;
253  }
254  }
255  // Scan rows in reverse order
256  for (r=nrow-1; r>=0; --r) {
257  pinv[r] = -1; // row r not yet ordered
258  c=leftmost[r];
259 
260  if (c==-1) continue; // row r is empty
261  if (nque[c]++ == 0) tail[c]=r; // first row in queue c
262  next[r] = head[c]; // put r at head of queue c
263  head[c] = r;
264  }
265  // Find row permutation and nnz(V)
266  casadi_int v_nnz = 0;
267  casadi_int nrow_new = nrow;
268  for (c=0; c<ncol; ++c) {
269  r = head[c]; // remove r from queue c
270  v_nnz++; // count V(c,c) as nonzero
271  if (r<0) r=nrow_new++; // add a fictitious row
272  pinv[r] = c; // associate row r with V(:,c)
273  if (--nque[c]<=0) continue; // skip if V(c+1,nrow,c) is empty
274  v_nnz += nque[c]; // nque[c] is nnz(V(c+1:nrow, c))
275  if ((pa=parent[c]) != -1) {
276  // Move all rows to parent of c
277  if (nque[pa]==0) tail[pa] = tail[c];
278  next[tail[c]] = head[pa];
279  head[pa] = next[r];
280  nque[pa] += nque[c];
281  }
282  }
283  for (r=0; r<nrow; ++r) if (pinv[r]<0) pinv[r] = c++;
284  if (nrow_ext) *nrow_ext = nrow_new;
285  return v_nnz;
286  }

References colind(), row(), and sp().

Referenced by qr_init().

◆ qr_sparsities()

void casadi::SparsityInternal::qr_sparsities ( const casadi_int *  sp_a,
casadi_int  nrow_ext,
casadi_int *  sp_v,
casadi_int *  sp_r,
const casadi_int *  leftmost,
const casadi_int *  parent,
const casadi_int *  pinv,
casadi_int *  iw 
)
static

Ref: Chapter 5, Direct Methods for Sparse Linear Systems by Tim Davis Note: nrow <= nrow_ext <= nrow+ncol len[iw] = nrow_ext + ncol len[x] = nrow_ext sp_v = [nrow_ext, ncol, 0, 0, ...] len[3 + ncol + nnz_v] len[v] nnz_v sp_r = [nrow_ext, ncol, 0, 0, ...] len[3 + ncol + nnz_r] len[r] nnz_r len[beta] ncol Modified version of cs_qr in CSparse Copyright(c) Timothy A. Davis, 2006-2009 Licensed as a derivative work under the GNU LGPL

Extra doc: https://github.com/casadi/casadi/wiki/L_ev

Definition at line 305 of file sparsity_internal.cpp.

308  {
309  /* Modified version of cs_qr in CSparse
310  Copyright(c) Timothy A. Davis, 2006-2009
311  Licensed as a derivative work under the GNU LGPL
312  */
313  // Extract sparsities
314  casadi_int ncol = sp_a[1];
315  const casadi_int *colind=sp_a+2, *row=sp_a+2+ncol+1;
316  casadi_int *v_colind=sp_v+2, *v_row=sp_v+2+ncol+1;
317  casadi_int *r_colind=sp_r+2, *r_row=sp_r+2+ncol+1;
318  // Specify dimensions of V and R
319  sp_v[0] = sp_r[0] = nrow_ext;
320  sp_v[1] = sp_r[1] = ncol;
321  // Work vectors
322  casadi_int* s = iw; iw += ncol;
323  // Local variables
324  casadi_int r, c, k, k1, top, len, k2, r2;
325  // Clear w to mark nodes
326  for (r=0; r<nrow_ext; ++r) iw[r] = -1;
327  // Number of nonzeros in v and r
328  casadi_int nnz_r=0, nnz_v=0;
329  // Compute V and R
330  for (c=0; c<ncol; ++c) {
331  // R(:,c) starts here
332  r_colind[c] = nnz_r;
333  // V(:, c) starts here
334  v_colind[c] = k1 = nnz_v;
335  // Add V(c,c) to pattern of V
336  iw[c] = c;
337  v_row[nnz_v++] = c;
338  top = ncol;
339  for (k=colind[c]; k<colind[c+1]; ++k) {
340  r = leftmost[row[k]]; // r = min(find(A(r,:))
341  // Traverse up c
342  for (len=0; iw[r]!=c; r=parent[r]) {
343  s[len++] = r;
344  iw[r] = c;
345  }
346  while (len>0) s[--top] = s[--len]; // push path on stack
347  r = pinv[row[k]]; // r = permuted row of A(:,c)
348  if (r>c && iw[r]<c) {
349  v_row[nnz_v++] = r; // add r to pattern of V(:,c)
350  iw[r] = c;
351  }
352  }
353  // For each r in pattern of R(:,c)
354  for (k = top; k<ncol; ++k) {
355  // R(r,c) is nonzero
356  r = s[k];
357  // Apply (V(r), beta(r)) to x: x -= v*beta*v'*x
358  r_row[nnz_r++] = r;
359  if (parent[r]==c) {
360  for (k2=v_colind[r]; k2<v_colind[r+1]; ++k2) {
361  r2 = v_row[k2];
362  if (iw[r2]<c) {
363  iw[r2] = c;
364  v_row[nnz_v++] = r2;
365  }
366  }
367  }
368  }
369  // R(c,c) = norm(x)
370  r_row[nnz_r++] = c;
371  }
372  // Finalize R, V
373  r_colind[ncol] = nnz_r;
374  v_colind[ncol] = nnz_v;
375  }

References colind(), and row().

Referenced by casadi::Sparsity::qr_sparse().

◆ randperm()

std::vector< casadi_int > casadi::SparsityInternal::randperm ( casadi_int  n,
casadi_int  seed 
)
static

return a random permutation vector, the identity perm, or p = n-1:-1:0. seed = -1 means p = n-1:-1:0. seed = 0 means p = identity. The implementation is a modified version of cs_randperm in CSparse Copyright(c) Timothy A. Davis, 2006-2009 Licensed as a derivative work under the GNU LGPL

Extra doc: https://github.com/casadi/casadi/wiki/L_f7

Definition at line 1408 of file sparsity_internal.cpp.

1408  {
1409  /*
1410  Modified version of cs_randperm in CSparse
1411  Copyright(c) Timothy A. Davis, 2006-2009
1412  Licensed as a derivative work under the GNU LGPL
1413  */
1414  // Return object
1415  std::vector<casadi_int> p;
1416 
1417  // return p = empty (identity)
1418  if (seed==0) return p;
1419 
1420  // allocate result
1421  p.resize(n);
1422 
1423  for (casadi_int k=0; k<n; ++k)
1424  p[k] = n-k-1;
1425 
1426  // return reverse permutation
1427  if (seed==-1) return p;
1428  #if defined(_WIN32)
1429  srand(seed);
1430  #else
1431  unsigned int seedu = static_cast<unsigned int>(seed);
1432  #endif
1433 
1434  for (casadi_int k=0; k<n; ++k) {
1435  // j = rand casadi_int in range k to n-1
1436  #if defined(_WIN32)
1437  casadi_int j = k + (rand() % (n-k)); // NOLINT(runtime/threadsafe_fn)
1438  #else
1439  casadi_int j = k + (rand_r(&seedu) % (n-k));
1440  #endif
1441  // swap p[k] and p[j]
1442  casadi_int t = p[j];
1443  p[j] = p[k];
1444  p[k] = t;
1445  }
1446 
1447  return p;
1448  }

Referenced by maxtrans().

◆ repr_el()

std::string casadi::SparsityInternal::repr_el ( casadi_int  k) const

Definition at line 1709 of file sparsity_internal.cpp.

1709  {
1710  casadi_int start_index = GlobalOptions::start_index;
1711  std::stringstream ss;
1712  if (numel()!=nnz()) {
1713  ss << "nonzero index " << k+start_index << " ";
1714  }
1715  casadi_int r = row()[k];
1716  casadi_int c = get_col()[k];
1717  ss << "(row " << r+start_index << ", col " << c+start_index << ")";
1718 
1719  return ss.str();
1720  }
casadi::GlobalOptions::start_index
static casadi_int start_index
Definition: global_options.hpp:71

References get_col(), nnz(), numel(), row(), and casadi::GlobalOptions::start_index.

◆ row()

const casadi_int* casadi::SparsityInternal::row ( ) const
inline

Extra doc: https://github.com/casadi/casadi/wiki/L_ej

Definition at line 110 of file sparsity_internal.hpp.

110 { return colind()+size2()+1;}

Referenced by _mtimes(), _reshape(), _resize(), _tril(), _triu(), amd(), augment(), bfs(), bw_lower(), bw_upper(), combineGen(), dfs(), drop(), drop_diag(), etree(), export_code(), find(), get_diag(), get_lower(), get_nz(), get_row(), get_upper(), has_diag(), hash(), is_diag(), is_equal(), is_reshape(), is_stacked(), is_transpose(), is_tril(), is_triu(), ldl_colind(), ldl_row(), makeDense(), maxtrans(), nnz_diag(), nnz_lower(), nnz_upper(), pattern_inverse(), permute(), pmult(), qr_nnz(), qr_sparsities(), repr_el(), rowsSequential(), scatter(), SparsityInternal(), spsolve(), spy(), star_coloring(), star_coloring2(), sub(), and uni_coloring().

◆ rowsSequential()

bool casadi::SparsityInternal::rowsSequential ( bool  strictly) const

Definition at line 2765 of file sparsity_internal.cpp.

2765  {
2766  const casadi_int* colind = this->colind();
2767  const casadi_int* row = this->row();
2768  for (casadi_int i=0; i<size2(); ++i) {
2769  casadi_int lastrow = -1;
2770  for (casadi_int k=colind[i]; k<colind[i+1]; ++k) {
2771 
2772  // check if not in sequence
2773  if (row[k] < lastrow)
2774  return false;
2775 
2776  // Check if duplicate
2777  if (strictly && row[k] == lastrow)
2778  return false;
2779 
2780  // update last row of the col
2781  lastrow = row[k];
2782  }
2783  }
2784 
2785  // sequential if reached this point
2786  return true;
2787  }

References colind(), row(), and size2().

◆ rprune()

casadi_int casadi::SparsityInternal::rprune ( casadi_int  i,
casadi_int  j,
double  aij,
void *  other 
)
static

The implementation is a modified version of cs_rprune in CSparse Copyright(c) Timothy A. Davis, 2006-2009 Licensed as a derivative work under the GNU LGPL

Extra doc: https://github.com/casadi/casadi/wiki/L_f2

Definition at line 1070 of file sparsity_internal.cpp.

1070  {
1071  /*
1072  Modified version of cs_rprune in CSparse
1073  Copyright(c) Timothy A. Davis, 2006-2009
1074  Licensed as a derivative work under the GNU LGPL
1075  */
1076  std::vector<casadi_int> &rr = *static_cast<std::vector<casadi_int> *>(other);
1077  return (i >= rr[1] && i < rr[2]) ;
1078  }

Referenced by dmperm().

◆ scatter()

casadi_int casadi::SparsityInternal::scatter ( casadi_int  j,
std::vector< casadi_int > &  w,
casadi_int  mark,
casadi_int *  Ci,
casadi_int  nz 
) const

The implementation is a modified version of cs_scatter in CSparse Copyright(c) Timothy A. Davis, 2006-2009 Licensed as a derivative work under the GNU LGPL

Extra doc: https://github.com/casadi/casadi/wiki/L_fd

Definition at line 1549 of file sparsity_internal.cpp.

1550  {
1551  /*
1552  Modified version of cs_scatter in CSparse
1553  Copyright(c) Timothy A. Davis, 2006-2009
1554  Licensed as a derivative work under the GNU LGPL
1555  */
1556  casadi_int i, p;
1557  const casadi_int *Ap = colind();
1558  const casadi_int *Ai = row();
1559 
1560  for (p = Ap[j]; p<Ap[j+1]; ++p) {
1561  // A(i, j) is nonzero
1562  i = Ai[p];
1563 
1564  if (w[i] < mark) {
1565  // i is new entry in row j
1566  w[i] = mark;
1567 
1568  // add i to pattern of C(:, j)
1569  Ci[nz++] = i;
1570  }
1571  }
1572  return nz;
1573  }

References colind(), and row().

Referenced by multiply().

◆ scc()

casadi_int casadi::SparsityInternal::scc ( std::vector< casadi_int > &  p,
std::vector< casadi_int > &  r 
) const

The implementation is a modified version of cs_scc in CSparse Copyright(c) Timothy A. Davis, 2006-2009 Licensed as a derivative work under the GNU LGPL

Extra doc: https://github.com/casadi/casadi/wiki/L_em

Definition at line 583 of file sparsity_internal.cpp.

584  {
585  /*
586  Modified version of cs_scc in CSparse
587  Copyright(c) Timothy A. Davis, 2006-2009
588  Licensed as a derivative work under the GNU LGPL
589  */
590  std::vector<casadi_int> tmp;
591 
592  Sparsity AT = T();
593 
594  std::vector<casadi_int> xi(2*size2()+1);
595  std::vector<casadi_int>& Blk = xi;
596 
597  std::vector<casadi_int> pstack(size2()+1);
598 
599  p.resize(size2());
600  r.resize(size2()+6);
601 
602  std::vector<bool> marked(size2(), false);
603 
604  casadi_int top = size2();
605 
606  //first dfs(A) to find finish times (xi)
607  for (casadi_int i = 0; i<size2(); ++i) {
608  if (!marked[i])
609  top = dfs(i, top, xi, pstack, tmp, marked);
610  }
611 
612  //restore A; unmark all nodes
613  std::fill(marked.begin(), marked.end(), false);
614 
615  top = size2();
616  casadi_int nb = size2();
617 
618  // dfs(A') to find strongly connnected comp
619  for (casadi_int k=0 ; k < size2() ; ++k) {
620  // get i in reverse order of finish times
621  casadi_int i = xi[k];
622 
623  // skip node i if already ordered
624  if (marked[i]) continue;
625 
626  // node i is the start of a component in p
627  r[nb--] = top;
628  top = AT.dfs(i, top, p, pstack, tmp, marked);
629  }
630 
631  // first block starts at zero; shift r up
632  r[nb] = 0;
633  for (casadi_int k = nb ; k <= size2() ; ++k)
634  r[k-nb] = r[k] ;
635 
636  // nb = # of strongly connected components
637  nb = size2()-nb;
638 
639  // sort each block in natural order
640  for (casadi_int b = 0 ; b < nb ; b++) {
641  for (casadi_int k = r[b]; k<r[b+1] ; ++k)
642  Blk[p[k]] = b ;
643  }
644 
645  // Get p; shift r down (side effect)
646  for (casadi_int i=0; i<size2(); ++i) {
647  p[r[Blk[i]]++] = i;
648  }
649 
650  // Shift up r
651  r.resize(nb+1);
652  for (casadi_int i=nb; i>0; --i) {
653  r[i]=r[i-1];
654  }
655  r[0]=0;
656 
657  return nb;
658  }
casadi::SparsityInternal::dfs
casadi_int dfs(casadi_int j, casadi_int top, std::vector< casadi_int > &xi, std::vector< casadi_int > &pstack, const std::vector< casadi_int > &pinv, std::vector< bool > &marked) const
Depth-first search.
Definition: sparsity_internal.cpp:518

References casadi::Sparsity::dfs(), dfs(), size2(), and T().

◆ shared_from_this() [1/2]

B casadi::GenericSharedInternal< SharedObject , SharedObjectInternal >::shared_from_this ( )
inlineprotectedinherited

Definition at line 83 of file generic_shared_internal.hpp.

83  {
84  casadi_assert_dev(B::test_cast(static_cast<Internal*>(this)));
85  B ret;
86  ret.own(static_cast<Internal*>(this));
87  return ret;
88  }

◆ shared_from_this() [2/2]

const B casadi::GenericSharedInternal< SharedObject , SharedObjectInternal >::shared_from_this ( ) const
inlineprotectedinherited

Definition at line 92 of file generic_shared_internal.hpp.

92  {
93  casadi_assert_dev(B::test_cast(static_cast<const Internal*>(this)));
94  B ret;
95  ret.own(const_cast<Internal*>(static_cast<const Internal*>(this)));
96  return ret;
97  }

◆ size()

std::pair< casadi_int, casadi_int > casadi::SparsityInternal::size ( ) const

Definition at line 1919 of file sparsity_internal.cpp.

1919  {
1920  return std::pair<casadi_int, casadi_int>(size1(), size2());
1921  }

References size1(), and size2().

Referenced by _mtimes(), and combineGen().

◆ size1()

casadi_int casadi::SparsityInternal::size1 ( ) const
inline

Extra doc: https://github.com/casadi/casadi/wiki/L_eg

Definition at line 95 of file sparsity_internal.hpp.

95 { return sp_[0];}

Referenced by _appendColumns(), _appendVector(), _enlargeColumns(), _enlargeRows(), _erase(), _mtimes(), _removeDuplicates(), _reshape(), _tril(), _triu(), combineGen(), dim(), dmperm(), drop_diag(), export_code(), find(), get_diag(), get_nz(), has_diag(), hash(), is_diag(), is_empty(), is_equal(), is_orthonormal(), is_orthonormal_columns(), is_orthonormal_rows(), is_reshape(), is_row(), is_scalar(), is_square(), is_stacked(), is_transpose(), makeDense(), maxtrans(), multiply(), numel(), pattern_inverse(), permute(), pmult(), size(), spy(), sub(), and transpose().

◆ size2()

casadi_int casadi::SparsityInternal::size2 ( ) const
inline

Extra doc: https://github.com/casadi/casadi/wiki/L_eh

Definition at line 100 of file sparsity_internal.hpp.

100 { return sp_[1];}

Referenced by _appendColumns(), _appendVector(), _enlargeRows(), _erase(), _removeDuplicates(), _reshape(), _resize(), _tril(), _triu(), amd(), bw_lower(), bw_upper(), combineGen(), dim(), dmperm(), drop_diag(), export_code(), find(), get_col(), get_colind(), get_diag(), get_lower(), get_nz(), get_upper(), has_diag(), hash(), is_column(), is_diag(), is_empty(), is_equal(), is_orthonormal_columns(), is_orthonormal_rows(), is_reshape(), is_scalar(), is_square(), is_stacked(), is_transpose(), is_tril(), is_triu(), largest_first(), makeDense(), maxtrans(), multiply(), nnz_diag(), nnz_lower(), nnz_upper(), numel(), pattern_inverse(), permute(), pmult(), rowsSequential(), scc(), size(), spy(), star_coloring(), star_coloring2(), sub(), transpose(), and uni_coloring().

◆ sp()

const std::vector<casadi_int>& casadi::SparsityInternal::sp ( ) const
inline

Extra doc: https://github.com/casadi/casadi/wiki/L_ef

Definition at line 90 of file sparsity_internal.hpp.

90 { return sp_;}

Referenced by _appendColumns(), _appendVector(), etree(), is_orthonormal(), is_orthonormal_columns(), is_orthonormal_rows(), ldl_colind(), ldl_row(), qr_init(), qr_nnz(), and sub().

◆ spsolve()

void casadi::SparsityInternal::spsolve ( bvec_t *  X,
bvec_t *  B,
bool  tr 
) const

Definition at line 3845 of file sparsity_internal.cpp.

3846  {
3847  const Btf& btf = this->btf();
3848  const casadi_int* colind = this->colind();
3849  const casadi_int* row = this->row();
3850 
3851  if (!tr) {
3852  for (casadi_int b = 0; b < btf.nb; ++b) { // loop over the blocks forward
3853 
3854  // Get dependencies from all right-hand-sides in the block ...
3855  bvec_t block_dep = 0;
3856  for (casadi_int el=btf.rowblock[b]; el<btf.rowblock[b+1]; ++el) {
3857  casadi_int rr = btf.rowperm[el];
3858  block_dep |= B[rr];
3859  }
3860 
3861  // ... as well as all other variables in the block
3862  for (casadi_int el=btf.colblock[b]; el<btf.colblock[b+1]; ++el) {
3863  casadi_int cc = btf.colperm[el];
3864  block_dep |= X[cc];
3865  }
3866 
3867  // Propagate ...
3868  for (casadi_int el=btf.colblock[b]; el<btf.colblock[b+1]; ++el) {
3869  casadi_int cc = btf.colperm[el];
3870 
3871  // ... to all variables in the block ...
3872  X[cc] |= block_dep;
3873 
3874  // ... as well as to all right-hand-sides that depend on variables in the block
3875  for (casadi_int k=colind[cc]; k<colind[cc+1]; ++k) {
3876  casadi_int rr=row[k];
3877  B[rr] |= block_dep;
3878  }
3879  }
3880  }
3881 
3882  } else { // transpose
3883  for (casadi_int b = btf.nb; b-- > 0; ) { // loop over the blocks backward
3884 
3885  // Get dependencies ...
3886  bvec_t block_dep = 0;
3887  for (casadi_int el=btf.colblock[b]; el<btf.colblock[b+1]; ++el) {
3888  casadi_int cc = btf.colperm[el];
3889 
3890  // .. from all right-hand-sides in the block ...
3891  block_dep |= B[cc];
3892 
3893  // ... as well as from all depending variables ...
3894  for (casadi_int k=colind[cc]; k<colind[cc+1]; ++k) {
3895  casadi_int rr=row[k];
3896  block_dep |= X[rr];
3897  }
3898  }
3899 
3900  // Propagate to all right-hand-sides in the block ...
3901  for (casadi_int el=btf.colblock[b]; el<btf.colblock[b+1]; ++el) {
3902  casadi_int cc = btf.colperm[el];
3903  B[cc] |= block_dep;
3904  }
3905 
3906  // ... as well as all other variables in the block
3907  for (casadi_int el=btf.rowblock[b]; el<btf.rowblock[b+1]; ++el) {
3908  casadi_int rr = btf.rowperm[el];
3909  X[rr] |= block_dep;
3910  }
3911  }
3912  }
3913  }
casadi::bvec_t
unsigned long long bvec_t
Definition: casadi_common.hpp:76

References btf(), colind(), row(), and casadi::X.

◆ spy()

void casadi::SparsityInternal::spy ( std::ostream &  stream) const

Extra doc: https://github.com/casadi/casadi/wiki/L_fq

Definition at line 3581 of file sparsity_internal.cpp.

3581  {
3582 
3583  // Index counter for each column
3584  std::vector<casadi_int> cind = get_colind();
3585  const casadi_int* colind = this->colind();
3586  const casadi_int* row = this->row();
3587 
3588  // Loop over rows
3589  for (casadi_int rr=0; rr<size1(); ++rr) {
3590 
3591  // Loop over columns
3592  for (casadi_int cc=0; cc<size2(); ++cc) {
3593  // Check if nonzero
3594  if (cind[cc]<colind[cc+1] && row[cind[cc]]==rr) {
3595  stream << "*";
3596  cind[cc]++;
3597  } else {
3598  stream << ".";
3599  }
3600  }
3601 
3602  // End of row
3603  stream << std::endl;
3604  }
3605  }

References colind(), get_colind(), row(), size1(), and size2().

◆ spy_matlab()

void casadi::SparsityInternal::spy_matlab ( const std::string &  mfile) const

Definition at line 3693 of file sparsity_internal.cpp.

3693  {
3694  // Create the .m file
3695  std::ofstream mfile;
3696  mfile.open(mfile_name.c_str());
3697 
3698  // Header
3699  mfile << "% This function was automatically generated by CasADi" << std::endl;
3700 
3701  Dict opts;
3702  opts["name"] = "A";
3703  export_code("matlab", mfile, opts);
3704 
3705  // Issue spy command
3706  mfile << "spy(A);" << std::endl;
3707 
3708  mfile.close();
3709  }
casadi::SparsityInternal::export_code
void export_code(const std::string &lang, std::ostream &stream, const Dict &options) const
Export sparsity in Matlab format.
Definition: sparsity_internal.cpp:3607
casadi::Dict
GenericType::Dict Dict
C++ equivalent of Python's dict or MATLAB's struct.
Definition: generic_type.hpp:258

References export_code().

◆ star_coloring()

Sparsity casadi::SparsityInternal::star_coloring ( casadi_int  ordering,
casadi_int  cutoff 
) const

See description in public class.

Extra doc: https://github.com/casadi/casadi/wiki/L_fo

Definition at line 3293 of file sparsity_internal.cpp.

3293  {
3294  if (!is_square()) {
3295  // NOTE(@jaeandersson) Why warning and not error?
3296  casadi_message("StarColoring requires a square matrix, got " + dim() + ".");
3297  }
3298 
3299  // Reorder, if necessary
3300  if (ordering!=0) {
3301  casadi_assert_dev(ordering==1);
3302 
3303  // Ordering
3304  std::vector<casadi_int> ord = largest_first();
3305 
3306  // Create a new sparsity pattern
3307  Sparsity sp_permuted = pmult(ord, true, true, true);
3308 
3309  // Star coloring for the permuted matrix
3310  Sparsity ret_permuted = sp_permuted.star_coloring(0);
3311 
3312  // Permute result back
3313  return ret_permuted.pmult(ord, true, false, false);
3314  }
3315 
3316  // Allocate temporary vectors
3317  const casadi_int* colind = this->colind();
3318  const casadi_int* row = this->row();
3319  std::vector<casadi_int> forbiddenColors;
3320  forbiddenColors.reserve(size2());
3321  std::vector<casadi_int> color(size2(), -1);
3322 
3323  // 4: for i <- 1 to |V | do
3324  for (casadi_int i=0; i<size2(); ++i) {
3325 
3326  // 5: for each w \in N1 (vi) do
3327  for (casadi_int w_el=colind[i]; w_el<colind[i+1]; ++w_el) {
3328  casadi_int w = row[w_el];
3329 
3330  // 6: if w is colored then
3331  if (color[w]!=-1) {
3332 
3333  // 7: forbiddenColors[color[w]] <- v
3334  forbiddenColors[color[w]] = i;
3335 
3336  } // 8: end if
3337 
3338  // 9: for each colored vertex x \in N1 (w) do
3339  for (casadi_int x_el=colind[w]; x_el<colind[w+1]; ++x_el) {
3340  casadi_int x = row[x_el];
3341  if (color[x]==-1) continue;
3342 
3343  // 10: if w is not colored then
3344  if (color[w]==-1) {
3345 
3346  //11: forbiddenColors[color[x]] <- vi
3347  forbiddenColors[color[x]] = i;
3348 
3349  } else { // 12: else
3350 
3351  // 13: for each colored vertex y \in N1 (x), y != w do
3352  for (casadi_int y_el=colind[x]; y_el<colind[x+1]; ++y_el) {
3353  casadi_int y = row[y_el];
3354  if (color[y]==-1 || y==w) continue;
3355 
3356  // 14: if color[y] = color[w] then
3357  if (color[y]==color[w]) {
3358 
3359  // 15: forbiddenColors[color[x]] <- vi
3360  forbiddenColors[color[x]] = i;
3361 
3362  // 16: break
3363  break;
3364 
3365  } // 17: end if
3366 
3367  } // 18: end for
3368 
3369  } // 19: end if
3370 
3371  } // 20 end for
3372 
3373  } // 21 end for
3374 
3375  // 22: color[v] <- min {c > 0 : forbiddenColors[c] = v}
3376  bool new_color = true;
3377  for (casadi_int color_i=0; color_i<forbiddenColors.size(); ++color_i) {
3378  // Break if color is ok
3379  if (forbiddenColors[color_i]!=i) {
3380  color[i] = color_i;
3381  new_color = false;
3382  break;
3383  }
3384  }
3385 
3386  // New color if reached end
3387  if (new_color) {
3388  color[i] = forbiddenColors.size();
3389  forbiddenColors.push_back(-1);
3390 
3391  // Cutoff if too many colors
3392  if (forbiddenColors.size()>cutoff) {
3393  return Sparsity();
3394  }
3395  }
3396 
3397  } // 23 end for
3398 
3399  // Number of colors used
3400  casadi_int num_colors = forbiddenColors.size();
3401 
3402  // Return sparsity in sparse triplet format
3403  return Sparsity::triplet(size2(), num_colors, range(color.size()), color);
3404  }
casadi::SparsityInternal::pmult
Sparsity pmult(const std::vector< casadi_int > &p, bool permute_rows=true, bool permute_cols=true, bool invert_permutation=false) const
Permute rows and/or columns.
Definition: sparsity_internal.cpp:3441
casadi::SparsityInternal::largest_first
std::vector< casadi_int > largest_first() const
Order the columns by decreasing degree.
Definition: sparsity_internal.cpp:3406

References colind(), dim(), is_square(), largest_first(), pmult(), casadi::Sparsity::pmult(), casadi::range(), row(), size2(), casadi::Sparsity::star_coloring(), and casadi::Sparsity::triplet().

Referenced by casadi::Sparsity::star_coloring().

◆ star_coloring2()

Sparsity casadi::SparsityInternal::star_coloring2 ( casadi_int  ordering,
casadi_int  cutoff 
) const

See description in public class.

Extra doc: https://github.com/casadi/casadi/wiki/L_fp

Definition at line 3040 of file sparsity_internal.cpp.

3040  {
3041  if (!is_square()) {
3042  // NOTE(@jaeandersson) Why warning and not error?
3043  casadi_message("StarColoring requires a square matrix, got " + dim() + ".");
3044  }
3045 
3046  // TODO(Joel): What we need here, is a distance-2 smallest last ordering
3047  // Reorder, if necessary
3048  const casadi_int* colind = this->colind();
3049  const casadi_int* row = this->row();
3050  if (ordering!=0) {
3051  casadi_assert_dev(ordering==1);
3052 
3053  // Ordering
3054  std::vector<casadi_int> ord = largest_first();
3055 
3056  // Create a new sparsity pattern
3057  Sparsity sp_permuted = pmult(ord, true, true, true);
3058 
3059  // Star coloring for the permuted matrix
3060  Sparsity ret_permuted = sp_permuted.star_coloring2(0);
3061 
3062  // Permute result back
3063  return ret_permuted.pmult(ord, true, false, false);
3064  }
3065 
3066  // Allocate temporary vectors
3067  std::vector<casadi_int> forbiddenColors;
3068  forbiddenColors.reserve(size2());
3069  std::vector<casadi_int> color(size2(), -1);
3070 
3071  std::vector<casadi_int> firstNeighborP(size2(), -1);
3072  std::vector<casadi_int> firstNeighborQ(size2(), -1);
3073  std::vector<casadi_int> firstNeighborQ_el(size2(), -1);
3074 
3075  std::vector<casadi_int> treated(size2(), -1);
3076  std::vector<casadi_int> hub(nnz_upper(), -1);
3077 
3078  std::vector<casadi_int> Tmapping;
3079  transpose(Tmapping);
3080 
3081  std::vector<casadi_int> star(nnz());
3082  casadi_int k = 0;
3083  for (casadi_int i=0; i<size2(); ++i) {
3084  for (casadi_int j_el=colind[i]; j_el<colind[i+1]; ++j_el) {
3085  casadi_int j = row[j_el];
3086  if (i<j) {
3087  star[j_el] = k;
3088  star[Tmapping[j]] = k;
3089  k++;
3090  }
3091  }
3092  }
3093 
3094 
3095 
3096  casadi_int starID = 0;
3097 
3098  // 3: for each v \in V do
3099  for (casadi_int v=0; v<size2(); ++v) {
3100 
3101  // 4: for each colored w \in N1(v) do
3102  for (casadi_int w_el=colind[v]; w_el<colind[v+1]; ++w_el) {
3103  casadi_int w = row[w_el];
3104  casadi_int colorW = color[w];
3105  if (colorW==-1) continue;
3106 
3107  // 5: forbiddenColors[color[w]] <- v
3108  forbiddenColors[colorW] = v;
3109 
3110  // 6: (p, q) <- firstNeighbor[color[w]]
3111  casadi_int p = firstNeighborP[colorW];
3112  casadi_int q = firstNeighborQ[colorW];
3113 
3114  // 7: if p = v then < Case 1
3115  if (v==p) {
3116 
3117  // 8: if treated[q] != v then
3118  if (treated[q]!=v) {
3119 
3120  // 9: treat(v, q) < forbid colors of neighbors of q
3121 
3122  // treat@2: for each colored x \in N1 (q) do
3123  for (casadi_int x_el=colind[q]; x_el<colind[q+1]; ++x_el) {
3124  casadi_int x = row[x_el];
3125  if (color[x]==-1) continue;
3126 
3127  // treat@3: forbiddenColors[color[x]] <- v
3128  forbiddenColors[color[x]] = v;
3129  }
3130 
3131  // treat@4: treated[q] <- v
3132  treated[q] = v;
3133 
3134  }
3135  // 10: treat(v, w) < forbid colors of neighbors of w
3136 
3137  // treat@2: for each colored x \in N1 (w) do
3138  for (casadi_int x_el=colind[w]; x_el<colind[w+1]; ++x_el) {
3139  casadi_int x = row[x_el];
3140  if (color[x]==-1) continue;
3141 
3142  // treat@3: forbiddenColors[color[x]] <- v
3143  forbiddenColors[color[x]] = v;
3144  }
3145 
3146  // treat@4: treated[w] <- v
3147  treated[w] = v;
3148 
3149  // 11: else
3150  } else {
3151 
3152  // 12: firstNeighbor[color[w]] <- (v, w)
3153  firstNeighborP[colorW] = v;
3154  firstNeighborQ[colorW] = w;
3155  firstNeighborQ_el[colorW] = w_el;
3156 
3157  // 13: for each colored vertex x \in N1 (w) do
3158  casadi_int x_el_end = colind[w+1];
3159  casadi_int x, colorx;
3160  for (casadi_int x_el=colind[w]; x_el < x_el_end; ++x_el) {
3161  x = row[x_el];
3162  colorx = color[x];
3163  if (colorx==-1 || x==v) continue;
3164 
3165  // 14: if x = hub[star[wx]] then potential Case 2
3166  if (hub[star[x_el]]==x) {
3167 
3168  // 15: forbiddenColors[color[x]] <- v
3169  forbiddenColors[colorx] = v;
3170 
3171  }
3172  }
3173  }
3174 
3175  }
3176 
3177  // 16: color[v] <- min {c > 0 : forbiddenColors[c] != v}
3178  bool new_color = true;
3179  for (casadi_int color_i=0; color_i<forbiddenColors.size(); ++color_i) {
3180  // Break if color is ok
3181  if (forbiddenColors[color_i]!=v) {
3182  color[v] = color_i;
3183  new_color = false;
3184  break;
3185  }
3186  }
3187 
3188  // New color if reached end
3189  if (new_color) {
3190  color[v] = forbiddenColors.size();
3191  forbiddenColors.push_back(-1);
3192 
3193  // Cutoff if too many colors
3194  if (forbiddenColors.size()>cutoff) {
3195  return Sparsity();
3196  }
3197  }
3198 
3199  // 17: updateStars(v)
3200 
3201  // updateStars@2: for each colored w \in N1 (v) do
3202  for (casadi_int w_el=colind[v]; w_el<colind[v+1]; ++w_el) {
3203  casadi_int w = row[w_el];
3204  casadi_int colorW = color[w];
3205  if (colorW==-1) continue;
3206 
3207  // updateStars@3: if exits x \in N1 (w) where x = v and color[x] = color[v] then
3208  bool check = false;
3209  casadi_int x;
3210  casadi_int x_el;
3211  for (x_el=colind[w]; x_el<colind[w+1]; ++x_el) {
3212  x = row[x_el];
3213  if (x==v || color[x]!=color[v]) continue;
3214  check = true;
3215  break;
3216  }
3217  if (check) {
3218 
3219  // updateStars@4: hub[star[wx]] <- w
3220  casadi_int starwx = star[x_el];
3221  hub[starwx] = w;
3222 
3223  // updateStars@5: star[vw] <- star[wx]
3224  star[w_el] = starwx;
3225  star[Tmapping[w_el]] = starwx;
3226 
3227  // updateStars@6: else
3228  } else {
3229 
3230  // updateStars@7: (p, q) <- firstNeighbor[color[w]]
3231  casadi_int p = firstNeighborP[colorW];
3232  casadi_int q = firstNeighborQ[colorW];
3233  casadi_int q_el = firstNeighborQ_el[colorW];
3234 
3235  // updateStars@8: if (p = v) and (q = w) then
3236  if (p==v && q!=w) {
3237 
3238  // updateStars@9: hub[star[vq]] <- v
3239  casadi_int starvq = star[q_el];
3240  hub[starvq] = v;
3241 
3242  // updateStars@10: star[vw] <- star[vq]
3243  star[w_el] = starvq;
3244  star[Tmapping[w_el]] = starvq;
3245 
3246  // updateStars@11: else
3247  } else {
3248 
3249  // updateStars@12: starID <- starID + 1
3250  starID+= 1;
3251 
3252  // updateStars@13: star[vw] <- starID
3253  star[w_el] = starID;
3254  star[Tmapping[w_el]]= starID;
3255 
3256  }
3257 
3258  }
3259 
3260  }
3261 
3262  }
3263 
3264  // Create return sparsity containing the coloring
3265  std::vector<casadi_int> ret_colind(forbiddenColors.size()+1, 0), ret_row;
3266 
3267  // Get the number of rows for each col
3268  for (casadi_int i=0; i<color.size(); ++i) {
3269  ret_colind[color[i]+1]++;
3270  }
3271 
3272  // Cumsum
3273  for (casadi_int j=0; j<forbiddenColors.size(); ++j) {
3274  ret_colind[j+1] += ret_colind[j];
3275  }
3276 
3277  // Get row for each col
3278  ret_row.resize(color.size());
3279  for (casadi_int j=0; j<ret_row.size(); ++j) {
3280  ret_row[ret_colind[color[j]]++] = j;
3281  }
3282 
3283  // Swap index back one step
3284  for (casadi_int j=ret_colind.size()-2; j>=0; --j) {
3285  ret_colind[j+1] = ret_colind[j];
3286  }
3287  ret_colind[0] = 0;
3288 
3289  // Return the coloring
3290  return Sparsity(size2(), forbiddenColors.size(), ret_colind, ret_row);
3291  }
casadi::SparsityInternal::transpose
Sparsity transpose(std::vector< casadi_int > &mapping, bool invert_mapping=false) const
Transpose the matrix and get the reordering of the non-zero entries,.
Definition: sparsity_internal.cpp:508
casadi::SparsityInternal::nnz_upper
casadi_int nnz_upper(bool strictly=false) const
Number of non-zeros in the upper triangular half.
Definition: sparsity_internal.cpp:1907

References colind(), dim(), is_square(), largest_first(), nnz(), nnz_upper(), pmult(), casadi::Sparsity::pmult(), row(), size2(), casadi::Sparsity::star_coloring2(), and transpose().

Referenced by casadi::Sparsity::star_coloring2().

◆ sub() [1/2]

Sparsity casadi::SparsityInternal::sub ( const std::vector< casadi_int > &  rr,
const SparsityInternal &  sp,
std::vector< casadi_int > &  mapping,
bool  ind1 
) const

Does bounds checking rr is not required to be monotonous

Extra doc: https://github.com/casadi/casadi/wiki/L_fl

Definition at line 2152 of file sparsity_internal.cpp.

2153  {
2154  casadi_assert_dev(rr.size()==sp.nnz());
2155 
2156  // Check bounds
2157  casadi_assert_in_range(rr, -numel()+ind1, numel()+ind1);
2158 
2159  // Handle index-1, negative indices
2160  if (ind1 || has_negative(rr)) {
2161  std::vector<casadi_int> rr_mod = rr;
2162  for (std::vector<casadi_int>::iterator i=rr_mod.begin(); i!=rr_mod.end(); ++i) {
2163  casadi_assert(!(ind1 && (*i)<=0),
2164  "Matlab is 1-based, but requested index " + str(*i) + ". "
2165  "Note that negative slices are disabled in the Matlab interface. "
2166  "Possibly you may want to use 'end'.");
2167  if (ind1) (*i)--;
2168  if (*i<0) *i += numel();
2169  }
2170  return sub(rr_mod, sp, mapping, false); // Call recursively
2171  }
2172 
2173  // Find the nonzeros corresponding to rr
2174  mapping.resize(rr.size());
2175  std::copy(rr.begin(), rr.end(), mapping.begin());
2176  get_nz(mapping);
2177 
2178  // Construct new pattern of the corresponding elements
2179  std::vector<casadi_int> ret_colind(sp.size2()+1), ret_row;
2180  ret_colind[0] = 0;
2181  const casadi_int* sp_colind = sp.colind();
2182  const casadi_int* sp_row = sp.row();
2183  for (casadi_int c=0; c<sp.size2(); ++c) {
2184  for (casadi_int el=sp_colind[c]; el<sp_colind[c+1]; ++el) {
2185  if (mapping[el]>=0) {
2186  mapping[ret_row.size()] = mapping[el];
2187  ret_row.push_back(sp_row[el]);
2188  }
2189  }
2190  ret_colind[c+1] = ret_row.size();
2191  }
2192  mapping.resize(ret_row.size());
2193  return Sparsity(sp.size1(), sp.size2(), ret_colind, ret_row);
2194  }
casadi::SparsityInternal::sub
Sparsity sub(const std::vector< casadi_int > &rr, const std::vector< casadi_int > &cc, std::vector< casadi_int > &mapping, bool ind1) const
Get a submatrix.
Definition: sparsity_internal.cpp:2196

References get_nz(), casadi::has_negative(), numel(), sp(), casadi::str(), and sub().

◆ sub() [2/2]

Sparsity casadi::SparsityInternal::sub ( const std::vector< casadi_int > &  rr,
const std::vector< casadi_int > &  cc,
std::vector< casadi_int > &  mapping,
bool  ind1 
) const

Does bounds checking rr and rr are not required to be monotonous

Extra doc: https://github.com/casadi/casadi/wiki/L_fk

Definition at line 2196 of file sparsity_internal.cpp.

2198  {
2199  casadi_assert_in_range(rr, -size1()+ind1, size1()+ind1);
2200  casadi_assert_in_range(cc, -size2()+ind1, size2()+ind1);
2201 
2202  // Handle index-1, negative indices in rr
2203  std::vector<casadi_int> tmp = rr;
2204  for (std::vector<casadi_int>::iterator i=tmp.begin(); i!=tmp.end(); ++i) {
2205  if (ind1) (*i)--;
2206  if (*i<0) *i += size1();
2207  }
2208  std::vector<casadi_int> rr_sorted, rr_sorted_index;
2209  sort(tmp, rr_sorted, rr_sorted_index, false);
2210 
2211  // Handle index-1, negative indices in cc
2212  tmp = cc;
2213  for (std::vector<casadi_int>::iterator i=tmp.begin(); i!=tmp.end(); ++i) {
2214  if (ind1) (*i)--;
2215  if (*i<0) *i += size2();
2216  }
2217  std::vector<casadi_int> cc_sorted, cc_sorted_index;
2218  sort(tmp, cc_sorted, cc_sorted_index, false);
2219  std::vector<casadi_int> columns, rows;
2220 
2221  // With lookup vector
2222  bool with_lookup = static_cast<double>(cc.size())*static_cast<double>(rr.size()) > nnz();
2223  std::vector<casadi_int> rrlookup;
2224  if (with_lookup) {
2225  // Time complexity: O(ii.size()*(nnz per column))
2226  // Typical use case:
2227  // a = SX::sym("a", sp_diag(50000))
2228  // a[:, :]
2229  rrlookup = lookupvector(rr_sorted, size1());
2230  // Else: Time complexity: O(ii.size()*jj.size())
2231  // Typical use case:
2232  // a = DM.ones(1000, 1000)
2233  // a[[0, 1],[0, 1]]
2234  }
2235 
2236  // count the number of non-zeros
2237  casadi_int nnz = 0;
2238 
2239  // loop over the columns of the slice
2240  const casadi_int* colind = this->colind();
2241  const casadi_int* row = this->row();
2242  for (casadi_int i=0; i<cc.size(); ++i) {
2243  casadi_int it = cc_sorted[i];
2244  if (with_lookup) {
2245  // loop over the non-zeros of the matrix
2246  for (casadi_int el=colind[it]; el<colind[it+1]; ++el) {
2247  casadi_int j = row[el];
2248  casadi_int ji = rrlookup[j];
2249  if (ji!=-1) {
2250  casadi_int jv = rr_sorted[ji];
2251  while (ji>=0 && jv == rr_sorted[ji--]) nnz++;
2252  }
2253  }
2254  } else {
2255  // Loop over rr
2256  casadi_int el = colind[it];
2257  for (casadi_int j=0; j<rr_sorted.size(); ++j) {
2258  casadi_int jt=rr_sorted[j];
2259  // Continue to the non-zero element
2260  while (el<colind[it+1] && row[el]<jt) el++;
2261  // Add the non-zero element, if there was an element in the location exists
2262  if (el<colind[it+1] && row[el]== jt) nnz++;
2263  }
2264  }
2265  }
2266 
2267  mapping.resize(nnz);
2268  columns.resize(nnz);
2269  rows.resize(nnz);
2270 
2271  casadi_int k = 0;
2272  // loop over the columns of the slice
2273  for (casadi_int i=0; i<cc.size(); ++i) {
2274  casadi_int it = cc_sorted[i];
2275  if (with_lookup) {
2276  // loop over the non-zeros of the matrix
2277  for (casadi_int el=colind[it]; el<colind[it+1]; ++el) {
2278  casadi_int jt = row[el];
2279  casadi_int ji = rrlookup[jt];
2280  if (ji!=-1) {
2281  casadi_int jv = rr_sorted[ji];
2282  while (ji>=0 && jv == rr_sorted[ji]) {
2283  rows[k] = rr_sorted_index[ji];
2284  columns[k] = cc_sorted_index[i];
2285  mapping[k] = el;
2286  k++;
2287  ji--;
2288  }
2289  }
2290  }
2291  } else {
2292  // Loop over rr
2293  casadi_int el = colind[it];
2294  for (casadi_int j=0; j<rr_sorted.size(); ++j) {
2295  casadi_int jt=rr_sorted[j];
2296  // Continue to the non-zero element
2297  while (el<colind[it+1] && row[el]<jt) el++;
2298  // Add the non-zero element, if there was an element in the location exists
2299  if (el<colind[it+1] && row[el]== jt) {
2300  rows[k] = rr_sorted_index[j];
2301  columns[k] = cc_sorted_index[i];
2302  mapping[k] = el;
2303  k++;
2304  }
2305  }
2306  }
2307  }
2308 
2309  std::vector<casadi_int> sp_mapping;
2310  std::vector<casadi_int> mapping_ = mapping;
2311  Sparsity ret = Sparsity::triplet(rr.size(), cc.size(), rows, columns, sp_mapping, false);
2312 
2313  for (casadi_int i=0; i<mapping.size(); ++i)
2314  mapping[i] = mapping_[sp_mapping[i]];
2315 
2316  // Create sparsity pattern
2317  return ret;
2318  }
casadi::lookupvector
std::vector< casadi_int > lookupvector(const std::vector< casadi_int > &v, casadi_int size)
Returns a vector for quickly looking up entries of supplied list.
Definition: casadi_misc.cpp:155

References colind(), casadi::lookupvector(), nnz(), row(), size1(), size2(), casadi::sort(), and casadi::Sparsity::triplet().

Referenced by casadi::Sparsity::sub(), and sub().

◆ T()

Sparsity casadi::SparsityInternal::T ( ) const

Definition at line 501 of file sparsity_internal.cpp.

501  {
502  // Dummy mapping
503  std::vector<casadi_int> mapping;
504 
505  return transpose(mapping);
506  }

References transpose().

Referenced by bfs(), maxtrans(), scc(), and casadi::Sparsity::T().

◆ transpose()

Sparsity casadi::SparsityInternal::transpose ( std::vector< casadi_int > &  mapping,
bool  invert_mapping = false 
) const
Parameters
[out]mappingthe non-zeros of the original matrix for each non-zero of the new matrix

Extra doc: https://github.com/casadi/casadi/wiki/L_ey

Definition at line 508 of file sparsity_internal.cpp.

509  {
510  // Get the sparsity of the transpose in sparse triplet form
511  std::vector<casadi_int> trans_col = get_row();
512  std::vector<casadi_int> trans_row = get_col();
513 
514  // Create the sparsity pattern
515  return Sparsity::triplet(size2(), size1(), trans_row, trans_col, mapping, invert_mapping);
516  }

References get_col(), get_row(), size1(), size2(), and casadi::Sparsity::triplet().

Referenced by star_coloring2(), T(), and casadi::Sparsity::transpose().

◆ uni_coloring()

Sparsity casadi::SparsityInternal::uni_coloring ( const Sparsity &  AT,
casadi_int  cutoff 
) const

A greedy distance-2 coloring algorithm (Algorithm 3.1 in A. H. GEBREMEDHIN, F. MANNE, A. POTHEN)

Extra doc: https://github.com/casadi/casadi/wiki/L_fn

Definition at line 2951 of file sparsity_internal.cpp.

2951  {
2952 
2953  // Allocate temporary vectors
2954  std::vector<casadi_int> forbiddenColors;
2955  forbiddenColors.reserve(size2());
2956  std::vector<casadi_int> color(size2(), 0);
2957 
2958  // Access the sparsity of the transpose
2959  const casadi_int* AT_colind = AT.colind();
2960  const casadi_int* AT_row = AT.row();
2961  const casadi_int* colind = this->colind();
2962  const casadi_int* row = this->row();
2963 
2964  // Loop over columns
2965  for (casadi_int i=0; i<size2(); ++i) {
2966 
2967  // Loop over nonzero elements
2968  for (casadi_int el=colind[i]; el<colind[i+1]; ++el) {
2969 
2970  // Get row
2971  casadi_int c = row[el];
2972 
2973  // Loop over previous columns that have an element in row c
2974  for (casadi_int el_prev=AT_colind[c]; el_prev<AT_colind[c+1]; ++el_prev) {
2975 
2976  // Get the col
2977  casadi_int i_prev = AT_row[el_prev];
2978 
2979  // Escape loop if we have arrived at the current col
2980  if (i_prev>=i)
2981  break;
2982 
2983  // Get the color of the col
2984  casadi_int color_prev = color[i_prev];
2985 
2986  // Mark the color as forbidden for the current col
2987  forbiddenColors[color_prev] = i;
2988  }
2989  }
2990 
2991  // Get the first nonforbidden color
2992  casadi_int color_i;
2993  for (color_i=0; color_i<forbiddenColors.size(); ++color_i) {
2994  // Break if color is ok
2995  if (forbiddenColors[color_i]!=i) break;
2996  }
2997  color[i] = color_i;
2998 
2999  // Add color if reached end
3000  if (color_i==forbiddenColors.size()) {
3001  forbiddenColors.push_back(0);
3002 
3003  // Cutoff if too many colors
3004  if (forbiddenColors.size()>cutoff) {
3005  return Sparsity();
3006  }
3007  }
3008  }
3009 
3010  // Create return sparsity containing the coloring
3011  std::vector<casadi_int> ret_colind(forbiddenColors.size()+1, 0), ret_row;
3012 
3013  // Get the number of rows for each col
3014  for (casadi_int i=0; i<color.size(); ++i) {
3015  ret_colind[color[i]+1]++;
3016  }
3017 
3018  // Cumsum
3019  for (casadi_int j=0; j<forbiddenColors.size(); ++j) {
3020  ret_colind[j+1] += ret_colind[j];
3021  }
3022 
3023  // Get row for each col
3024  ret_row.resize(color.size());
3025  for (casadi_int j=0; j<ret_row.size(); ++j) {
3026  ret_row[ret_colind[color[j]]++] = j;
3027  }
3028 
3029  // Swap index back one step
3030  for (casadi_int j=ret_colind.size()-2; j>=0; --j) {
3031  ret_colind[j+1] = ret_colind[j];
3032  }
3033  ret_colind[0] = 0;
3034 
3035  // Return the coloring
3036  return Sparsity(size2(), forbiddenColors.size(), ret_colind, ret_row);
3037 ;
3038  }

References casadi::Sparsity::colind(), colind(), casadi::Sparsity::row(), row(), and size2().

Referenced by casadi::Sparsity::uni_coloring().

◆ unmatched()

void casadi::SparsityInternal::unmatched ( casadi_int  m,
const std::vector< casadi_int > &  wi,
std::vector< casadi_int > &  p,
std::vector< casadi_int > &  rr,
casadi_int  set 
)
static

The implementation is a modified version of cs_unmatched in CSparse Copyright(c) Timothy A. Davis, 2006-2009 Licensed as a derivative work under the GNU LGPL

Extra doc: https://github.com/casadi/casadi/wiki/L_f1

Definition at line 1055 of file sparsity_internal.cpp.

1056  {
1057  /*
1058  Modified version of cs_unmatched in CSparse
1059  Copyright(c) Timothy A. Davis, 2006-2009
1060  Licensed as a derivative work under the GNU LGPL
1061  */
1062  casadi_int i, kr = rr[set] ;
1063  for (i=0; i<m; i++)
1064  if (wi[i] == 0)
1065  p[kr++] = i;
1066 
1067  rr[set+1] = kr;
1068  }

Referenced by dmperm().

◆ wclear()

casadi_int casadi::SparsityInternal::wclear ( casadi_int  mark,
casadi_int  lemax,
casadi_int *  w,
casadi_int  n 
)
static

The implementation is a modified version of cs_wclear in CSparse Copyright(c) Timothy A. Davis, 2006-2009 Licensed as a derivative work under the GNU LGPL

Extra doc: https://github.com/casadi/casadi/wiki/L_fa

Definition at line 1525 of file sparsity_internal.cpp.

1526  {
1527  /*
1528  Modified version of cs_wclear in CSparse
1529  Copyright(c) Timothy A. Davis, 2006-2009
1530  Licensed as a derivative work under the GNU LGPL
1531  */
1532  if (mark < 2 || (mark + lemax < 0)) {
1533  for (casadi_int k = 0; k<n; ++k) if (w[k] != 0) w[k] = 1;
1534  mark = 2 ;
1535  }
1536  // at this point, w [0..n-1] < mark holds
1537  return mark;
1538  }

Referenced by amd().

◆ weak()

GenericWeakRef< SharedObject , SharedObjectInternal > * casadi::GenericSharedInternal< SharedObject , SharedObjectInternal >::weak
inherited

Extra doc: https://github.com/casadi/casadi/wiki/L_1ai

Definition at line 67 of file generic_shared_internal.hpp.

191  {
192  if (weak_ref_==nullptr) {
193  weak_ref_ = new GenericWeakRef<Shared, Internal>(static_cast<Internal*>(this));
194  }
195  return weak_ref_;
196  }
The documentation for this class was generated from the following files: