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 46 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...
 
bool is_compactible (std::vector< casadi_int > &row, std::vector< casadi_int > &col) const
 Check if the nonzero pattern is a Cartesian product (row x col) 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_new (std::vector< casadi_int > &which_color, const Dict &opts) const
 Perform a star 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 450 of file sparsity_internal.cpp.

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

References colind(), and row().

◆ ~SparsityInternal()

casadi::SparsityInternal::~SparsityInternal ( )
override

Definition at line 460 of file sparsity_internal.cpp.

460  {
461  delete btf_;
462  }

Member Function Documentation

◆ _appendColumns()

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

Definition at line 2585 of file sparsity_internal.cpp.

2585  {
2586  casadi_assert(size1()== sp.size1(),
2587  "_appendColumns(sp): row sizes must match but got " + str(size1())
2588  + " for lhs, and " + str(sp.size1()) + " for rhs.");
2589 
2590  // Append rows
2591  std::vector<casadi_int> new_row = get_row();
2592  const casadi_int* sp_row = sp.row();
2593  new_row.insert(new_row.end(), sp_row, sp_row+sp.nnz());
2594 
2595  // Get column indices
2596  std::vector<casadi_int> new_colind = get_colind();
2597  const casadi_int* sp_colind = sp.colind();
2598  new_colind.resize(size2() + sp.size2() + 1);
2599  for (casadi_int i = size2()+1; i<new_colind.size(); ++i)
2600  new_colind[i] = sp_colind[i-size2()] + nnz();
2601 
2602  return Sparsity(size1(), size2()+sp.size2(), new_colind, new_row);
2603  }
casadi::SparsityInternal::sp
const std::vector< casadi_int > & sp() const
Get number of rows (see public class)
Definition: sparsity_internal.hpp:92
casadi::SparsityInternal::size1
casadi_int size1() const
Get number of rows (see public class)
Definition: sparsity_internal.hpp:97
casadi::SparsityInternal::size2
casadi_int size2() const
Get number of columns (see public class)
Definition: sparsity_internal.hpp:102
casadi::SparsityInternal::get_colind
std::vector< casadi_int > get_colind() const
Get colind() as a vector.
Definition: sparsity_internal.cpp:3949
casadi::SparsityInternal::get_row
std::vector< casadi_int > get_row() const
Get row() as a vector.
Definition: sparsity_internal.cpp:3954
casadi::SparsityInternal::nnz
casadi_int nnz() const
Number of structural non-zeros.
Definition: sparsity_internal.hpp:115
casadi::str
std::string str(const T &v)
String representation, any type.
Definition: casadi_common.hpp:248

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 2564 of file sparsity_internal.cpp.

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

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 2605 of file sparsity_internal.cpp.

2606  {
2607  casadi_assert_in_range(cc, -ncol+ind1, ncol+ind1);
2608 
2609  // Handle index-1, negative indices
2610  if (ind1 || has_negative(cc)) {
2611  std::vector<casadi_int> cc_mod = cc;
2612  for (std::vector<casadi_int>::iterator i=cc_mod.begin(); i!=cc_mod.end(); ++i) {
2613  if (ind1) (*i)--;
2614  if (*i<0) *i += ncol;
2615  }
2616  return _enlargeColumns(ncol, cc_mod, false); // Call recursively
2617  }
2618 
2619  // Sparsify the columns
2620  std::vector<casadi_int> new_colind = get_colind();
2621  new_colind.resize(ncol+1, nnz());
2622 
2623  casadi_int ik=cc.back(); // need only to update from the last new index
2624  casadi_int nz=nnz(); // number of nonzeros up till this column
2625  for (casadi_int i=cc.size()-1; i>=0; --i) {
2626  // Update colindex for new columns
2627  for (; ik>cc[i]; --ik) {
2628  new_colind[ik] = nz;
2629  }
2630 
2631  // Update non-zero counter
2632  nz = new_colind[i];
2633 
2634  // Update colindex for old colums
2635  new_colind[cc[i]] = nz;
2636  }
2637 
2638  // Append zeros to the beginning
2639  for (; ik>=0; --ik) {
2640  new_colind[ik] = 0;
2641  }
2642  return Sparsity(size1(), ncol, new_colind, get_row());
2643  }
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:2605
casadi::has_negative
bool has_negative(const std::vector< T > &v)
Check if the vector has negative entries.
Definition: casadi_misc.hpp:801

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 2645 of file sparsity_internal.cpp.

2646  {
2647  casadi_assert_in_range(rr, -nrow+ind1, nrow+ind1);
2648 
2649  // Handle index-1, negative indices
2650  if (ind1 || has_negative(rr)) {
2651  std::vector<casadi_int> rr_mod = rr;
2652  for (std::vector<casadi_int>::iterator i=rr_mod.begin(); i!=rr_mod.end(); ++i) {
2653  if (ind1) (*i)--;
2654  if (*i<0) *i += nrow;
2655  }
2656  return _enlargeRows(nrow, rr_mod, false); // Call recursively
2657  }
2658 
2659  // Assert dimensions
2660  casadi_assert_dev(rr.size() == size1());
2661 
2662  // Begin by sparsify the rows
2663  std::vector<casadi_int> new_row = get_row();
2664  for (casadi_int k=0; k<nnz(); ++k) {
2665  new_row[k] = rr[new_row[k]];
2666  }
2667  return Sparsity(nrow, size2(), get_colind(), new_row);
2668  }
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:2645

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 1931 of file sparsity_internal.cpp.

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

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 2018 of file sparsity_internal.cpp.

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

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 1730 of file sparsity_internal.cpp.

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

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 2797 of file sparsity_internal.cpp.

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

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

◆ _reshape()

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

Definition at line 2714 of file sparsity_internal.cpp.

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

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 2747 of file sparsity_internal.cpp.

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

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

◆ _tril()

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

Definition at line 3861 of file sparsity_internal.cpp.

3861  {
3862  const casadi_int* colind = this->colind();
3863  const casadi_int* row = this->row();
3864  std::vector<casadi_int> ret_colind, ret_row;
3865  ret_colind.reserve(size2()+1);
3866  ret_colind.push_back(0);
3867  for (casadi_int cc=0; cc<size2(); ++cc) {
3868  for (casadi_int el=colind[cc]; el<colind[cc+1]; ++el) {
3869  casadi_int rr=row[el];
3870  if (rr>cc || (includeDiagonal && rr==cc)) {
3871  ret_row.push_back(rr);
3872  }
3873  }
3874  ret_colind.push_back(ret_row.size());
3875  }
3876  return Sparsity(size1(), size2(), ret_colind, ret_row);
3877  }

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

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

◆ _triu()

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

Definition at line 3879 of file sparsity_internal.cpp.

3879  {
3880  const casadi_int* colind = this->colind();
3881  const casadi_int* row = this->row();
3882  std::vector<casadi_int> ret_colind, ret_row;
3883  ret_colind.reserve(size2()+1);
3884  ret_colind.push_back(0);
3885  for (casadi_int cc=0; cc<size2(); ++cc) {
3886  for (casadi_int el=colind[cc]; el<colind[cc+1]; ++el) {
3887  casadi_int rr=row[el];
3888  if (rr<cc || (includeDiagonal && rr==cc)) {
3889  ret_row.push_back(rr);
3890  }
3891  }
3892  ret_colind.push_back(ret_row.size());
3893  }
3894  return Sparsity(size1(), size2(), ret_colind, ret_row);
3895  }

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 661 of file sparsity_internal.cpp.

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

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 1088 of file sparsity_internal.cpp.

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

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 975 of file sparsity_internal.cpp.

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

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 464 of file sparsity_internal.cpp.

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

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 367 of file sparsity_internal.hpp.

370  {
371  T()->dmperm(colperm, rowperm, colblock, rowblock,
372  coarse_colblock, coarse_rowblock);
373  return rowblock.size()-1;
374  }
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:1252

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 3936 of file sparsity_internal.cpp.

3936  {
3937  casadi_int bw = 0;
3938  const casadi_int* colind = this->colind();
3939  const casadi_int* row = this->row();
3940  for (casadi_int cc=0; cc<size2(); ++cc) {
3941  if (colind[cc] != colind[cc+1]) { // if there are any elements of the column
3942  casadi_int rr = row[colind[cc+1]-1];
3943  bw = std::max(bw, rr-cc);
3944  }
3945  }
3946  return bw;
3947  }

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 3923 of file sparsity_internal.cpp.

3923  {
3924  casadi_int bw = 0;
3925  const casadi_int* colind = this->colind();
3926  const casadi_int* row = this->row();
3927  for (casadi_int cc=0; cc<size2(); ++cc) {
3928  if (colind[cc] != colind[cc+1]) { // if there are any elements of the column
3929  casadi_int rr = row[colind[cc]];
3930  bw = std::max(bw, cc-rr);
3931  }
3932  }
3933  return bw;
3934  }

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

◆ class_name()

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

Implements casadi::SharedObjectInternal.

Definition at line 714 of file sparsity_internal.hpp.

714 {return "SparsityInternal";}

◆ colind()

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

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

Definition at line 107 of file sparsity_internal.hpp.

107 { 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(), star_coloring_new(), 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 2328 of file sparsity_internal.cpp.

2329  {
2330  static std::vector<unsigned char> mapping;
2331  return combineGen1<false>(y, f0x_is_zero, function0_is_zero, mapping);
2332  }

◆ 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 2334 of file sparsity_internal.cpp.

2336  {
2337  return combineGen1<true>(y, f0x_is_zero, function0_is_zero, mapping);
2338  }

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 2378 of file sparsity_internal.cpp.

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

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 2351 of file sparsity_internal.cpp.

2353  {
2354 
2355  // Quick return if identical
2356  if (is_equal(y)) {
2357  if (with_mapping) {
2358  mapping.resize(y.nnz());
2359  std::fill(mapping.begin(), mapping.end(), 1 | 2);
2360  }
2361  return y;
2362  }
2363 
2364  if (f0x_is_zero) {
2365  if (function0_is_zero) {
2366  return combineGen<with_mapping, true, true>(y, mapping);
2367  } else {
2368  return combineGen<with_mapping, true, false>(y, mapping);
2369  }
2370  } else if (function0_is_zero) {
2371  return combineGen<with_mapping, false, true>(y, mapping);
2372  } else {
2373  return combineGen<with_mapping, false, false>(y, mapping);
2374  }
2375  }
casadi::SparsityInternal::is_equal
bool is_equal(const Sparsity &y) const
Check if two sparsity patterns are the same.
Definition: sparsity_internal.cpp:2482

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.

175  {
176  // Note: i != this because of something something multiple inheritance
177  return str( (casadi_int)(i)) + "/" + static_cast<const Internal*>(this)->class_name();
178  }

◆ 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 519 of file sparsity_internal.cpp.

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

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 1548 of file sparsity_internal.cpp.

1548  {
1549  /*
1550  Modified version of cs_diag in CSparse
1551  Copyright(c) Timothy A. Davis, 2006-2009
1552  Licensed as a derivative work under the GNU LGPL
1553  */
1554  return (i != j) ;
1555  }

◆ dim()

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

Definition at line 1711 of file sparsity_internal.cpp.

1711  {
1712  std::string ret = str(size1()) + "x" + str(size2());
1713  if (with_nz) ret += "," + str(nnz()) + "nz";
1714  return ret;
1715  }

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 482 of file sparsity_internal.cpp.

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

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 1252 of file sparsity_internal.cpp.

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

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 1504 of file sparsity_internal.cpp.

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

References colind(), and row().

Referenced by dmperm().

◆ drop_diag()

Sparsity casadi::SparsityInternal::drop_diag ( ) const

Definition at line 1691 of file sparsity_internal.cpp.

1691  {
1692  casadi_int nrow = this->size1();
1693  casadi_int ncol = this->size2();
1694  const casadi_int* colind = this->colind();
1695  const casadi_int* row = this->row();
1696  // Return sparsity
1697  std::vector<casadi_int> ret_colind(ncol+1), ret_row;
1698  ret_colind[0] = 0;
1699  ret_row.reserve(nnz());
1700  for (casadi_int c=0; c<ncol; ++c) {
1701  for (casadi_int k=colind[c]; k<colind[c+1]; ++k) {
1702  if (row[k]!=c) {
1703  ret_row.push_back(row[k]);
1704  }
1705  }
1706  ret_colind[c+1] = ret_row.size();
1707  }
1708  return Sparsity(nrow, ncol, ret_colind, ret_row);
1709  }

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 3723 of file sparsity_internal.cpp.

3724  {
3725  casadi_assert(lang=="matlab", "Only matlab language supported for now.");
3726 
3727  // Default values for options
3728  bool opt_inline = false;
3729  std::string name = "sp";
3730  bool as_matrix = true;
3731  casadi_int indent_level = 0;
3732  std::vector<std::string> nonzeros;
3733 
3734  // Read options
3735  for (auto&& op : options) {
3736  if (op.first=="inline") {
3737  opt_inline = op.second;
3738  } else if (op.first=="name") {
3739  name = op.second.to_string();
3740  } else if (op.first=="as_matrix") {
3741  as_matrix = op.second;
3742  } else if (op.first=="indent_level") {
3743  indent_level = op.second;
3744  } else if (op.first=="nonzeros") {
3745  nonzeros = op.second;
3746  } else {
3747  casadi_error("Unknown option '" + op.first + "'.");
3748  }
3749  }
3750 
3751  // Construct indent string
3752  std::string indent;
3753  for (casadi_int i=0; i < indent_level; ++i) indent += " ";
3754  casadi_assert(!opt_inline, "Inline not supported for now.");
3755 
3756  // Export dimensions
3757  stream << indent << name << "_m = " << size1() << ";\n";
3758  stream << indent << name << "_n = " << size2() << ";\n";
3759 
3760  // Matlab indices are one-based
3761  const casadi_int index_offset = 1;
3762 
3763  // Print columns
3764  const casadi_int* colind = this->colind();
3765  const casadi_int* row = this->row();
3766  stream << indent << name<< "_j = [";
3767  bool first = true;
3768  for (casadi_int i=0; i<size2(); ++i) {
3769  for (casadi_int el=colind[i]; el<colind[i+1]; ++el) {
3770  if (!first) stream << ", ";
3771  stream << (i+index_offset);
3772  first = false;
3773  }
3774  }
3775  stream << "];\n";
3776 
3777  // Print rows
3778  stream << indent << name << "_i = [";
3779  first = true;
3780  casadi_int nz = nnz();
3781  for (casadi_int i=0; i<nz; ++i) {
3782  if (!first) stream << ", ";
3783  stream << (row[i]+index_offset);
3784  first = false;
3785  }
3786  stream << "];\n";
3787 
3788  // Print nonzeros
3789  stream << indent << name << "_v = ";
3790  if (nonzeros.empty()) {
3791  stream << "ones(size(" << name << "_i));\n";
3792  } else {
3793  stream << "[";
3794  for (casadi_int i = 0; i < nonzeros.size(); ++i) {
3795  if (i > 0) stream << ", ";
3796  stream << nonzeros.at(i);
3797  }
3798  stream << "];\n";
3799  }
3800 
3801  if (as_matrix) {
3802  // Generate matrix
3803  stream << indent << name << " = sparse(" << name << "_i, " << name << "_j, ";
3804  stream << name << "_v, " << name << "_m, " << name << "_n);\n";
3805  }
3806 
3807  }

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 2848 of file sparsity_internal.cpp.

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

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

◆ get_col()

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

Definition at line 491 of file sparsity_internal.cpp.

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

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 3949 of file sparsity_internal.cpp.

3949  {
3950  const casadi_int* colind = this->colind();
3951  return std::vector<casadi_int>(colind, colind+size2()+1);
3952  }

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 1627 of file sparsity_internal.cpp.

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

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 3897 of file sparsity_internal.cpp.

3897  {
3898  const casadi_int* colind = this->colind();
3899  const casadi_int* row = this->row();
3900  std::vector<casadi_int> ret;
3901  for (casadi_int cc=0; cc<size2(); ++cc) {
3902  for (casadi_int el = colind[cc]; el<colind[cc+1]; ++el) {
3903  if (row[el]>=cc) {
3904  ret.push_back(el);
3905  }
3906  }
3907  }
3908  return ret;
3909  }

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 2684 of file sparsity_internal.cpp.

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

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 2126 of file sparsity_internal.cpp.

2127  {
2128  casadi_assert_bounded(rr, size1());
2129  casadi_assert_bounded(cc, size2());
2130 
2131  std::vector<casadi_int> rr_sorted;
2132  std::vector<casadi_int> rr_sorted_index;
2133 
2134  sort(rr, rr_sorted, rr_sorted_index);
2135 
2136  std::vector<casadi_int> ret(cc.size()*rr.size());
2137 
2138  casadi_int stride = rr.size();
2139  const casadi_int* colind = this->colind();
2140  const casadi_int* row = this->row();
2141 
2142  for (casadi_int i=0;i<cc.size();++i) {
2143  casadi_int it = cc[i];
2144  casadi_int el=colind[it];
2145  for (casadi_int j=0;j<rr_sorted.size();++j) {
2146  casadi_int jt=rr_sorted[j];
2147  // Continue to the non-zero element
2148  for (; el<colind[it+1] && row[el]<jt; ++el) {}
2149  // Add the non-zero element, if there was an element in the location exists
2150  if (el<colind[it+1] && row[el]== jt) {
2151  ret[i*stride+rr_sorted_index[j]] = el;
2152  } else {
2153  ret[i*stride+rr_sorted_index[j]] = -1;
2154  }
2155  }
2156  }
2157  return ret;
2158  }
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:906

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 2875 of file sparsity_internal.cpp.

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

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 3954 of file sparsity_internal.cpp.

3954  {
3955  const casadi_int* row = this->row();
3956  return std::vector<casadi_int>(row, row+nnz());
3957  }

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 3911 of file sparsity_internal.cpp.

3911  {
3912  const casadi_int* colind = this->colind();
3913  const casadi_int* row = this->row();
3914  std::vector<casadi_int> ret;
3915  for (casadi_int cc=0; cc<size2(); ++cc) {
3916  for (casadi_int el = colind[cc]; el<colind[cc+1] && row[el]<=cc; ++el) {
3917  ret.push_back(el);
3918  }
3919  }
3920  return ret;
3921  }

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

◆ getCount()

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

Definition at line 60 of file generic_shared_internal.hpp.

205  {
206  return static_cast<const Internal*>(this)->count;
207  }

◆ has_diag()

bool casadi::SparsityInternal::has_diag ( ) const

Definition at line 1678 of file sparsity_internal.cpp.

1678  {
1679  casadi_int nrow = this->size1();
1680  casadi_int ncol = this->size2();
1681  const casadi_int* colind = this->colind();
1682  const casadi_int* row = this->row();
1683  for (casadi_int c=0; c<ncol && c<nrow; ++c) {
1684  for (casadi_int k=colind[c]; k<colind[c+1]; ++k) {
1685  if (row[k]==c) return true;
1686  }
1687  }
1688  return false;
1689  }

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

◆ hash()

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

Definition at line 3825 of file sparsity_internal.cpp.

3825  {
3826  return hash_sparsity(size1(), size2(), colind(), row());
3827  }
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:1012

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 1458 of file sparsity_internal.cpp.

1458  {
1459  std::vector<casadi_int> pinv(p.size());
1460  for (casadi_int k=0; k<p.size(); ++k) pinv[p[k]] = k;
1461  return pinv;
1462  }

Referenced by dmperm().

◆ is_column()

bool casadi::SparsityInternal::is_column ( ) const

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

Definition at line 1801 of file sparsity_internal.cpp.

1801  {
1802  return size2()==1;
1803  }

References size2().

Referenced by is_vector().

◆ is_compactible()

bool casadi::SparsityInternal::is_compactible ( std::vector< casadi_int > &  row,
std::vector< casadi_int > &  col 
) const

Definition at line 3686 of file sparsity_internal.cpp.

3687  {
3688  // Nonempty rows: rows of sum2(*this) — an m x 1 sparse vector
3689  // Nonempty cols: cols of sum1(*this) — a 1 x n sparse vector
3690  Sparsity self = shared_from_this<Sparsity>();
3691  row = Sparsity::sum2(self).get_row();
3692  col = Sparsity::sum1(self).get_col();
3693  // Compactible iff every (nonempty row) x (nonempty col) is a nonzero
3694  return nnz() == static_cast<casadi_int>(row.size() * col.size());
3695  }
casadi::Sparsity::sum2
static Sparsity sum2(const Sparsity &x)
Enlarge matrix.
Definition: sparsity.cpp:1720
casadi::Sparsity::get_col
std::vector< casadi_int > get_col() const
Get the column for each non-zero entry.
Definition: sparsity.cpp:368
casadi::Sparsity::sum1
static Sparsity sum1(const Sparsity &x)
Enlarge matrix.
Definition: sparsity.cpp:1724
casadi::Sparsity::get_row
std::vector< casadi_int > get_row() const
Get the row for each non-zero entry.
Definition: sparsity.cpp:372

References casadi::Sparsity::get_col(), casadi::Sparsity::get_row(), nnz(), row(), casadi::Sparsity::sum1(), and casadi::Sparsity::sum2().

◆ is_dense()

bool casadi::SparsityInternal::is_dense ( ) const

Definition at line 1793 of file sparsity_internal.cpp.

1793  {
1794  return nnz() == numel();
1795  }

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 1813 of file sparsity_internal.cpp.

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

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 1809 of file sparsity_internal.cpp.

1809  {
1810  return both ? size2()==0 && size1()==0 : size2()==0 || size1()==0;
1811  }

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 2540 of file sparsity_internal.cpp.

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

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 2532 of file sparsity_internal.cpp.

2534  {
2535  casadi_assert_dev(y_colind.size()==y_ncol+1);
2536  casadi_assert_dev(y_row.size()==y_colind.back());
2537  return is_equal(y_nrow, y_ncol, get_ptr(y_colind), get_ptr(y_row));
2538  }

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

◆ is_equal() [3/3]

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

Definition at line 2482 of file sparsity_internal.cpp.

2482  {
2483  // Quick true if the objects are the same
2484  if (this == y.get()) return true;
2485 
2486  // Otherwise, compare the patterns
2487  return is_equal(y.size1(), y.size2(), y.colind(), y.row());
2488  }

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 1851 of file sparsity_internal.cpp.

1851  {
1852  if (!allow_empty) {
1853  if (!is_square()) return false;
1854  if (nnz()!=size1()) return false;
1855  }
1856 
1857  Sparsity sp = shared_from_this<Sparsity>();
1858  if (sum2(sp).nnz()!=nnz()) return false;
1859  if (sum1(sp).nnz()!=nnz()) return false;
1860  return true;
1861  }
casadi::SparsityInternal::is_square
bool is_square() const
Is square?
Definition: sparsity_internal.cpp:1839

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 1875 of file sparsity_internal.cpp.

1875  {
1876  if (!allow_empty) {
1877  if (size2()>size1()) return false;
1878  if (nnz()!=size2()) return false;
1879  }
1880 
1881  Sparsity sp = shared_from_this<Sparsity>();
1882  if (sum2(sp).nnz()!=nnz()) return false;
1883  if (sum1(sp).nnz()!=nnz()) return false;
1884  return true;
1885  }

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 1863 of file sparsity_internal.cpp.

1863  {
1864  if (!allow_empty) {
1865  if (size1()>size2()) return false;
1866  if (nnz()!=size1()) return false;
1867  }
1868 
1869  Sparsity sp = shared_from_this<Sparsity>();
1870  if (sum2(sp).nnz()!=nnz()) return false;
1871  if (sum1(sp).nnz()!=nnz()) return false;
1872  return true;
1873  }

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 1847 of file sparsity_internal.cpp.

1847  {
1848  return is_orthonormal();
1849  }
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:1851

References is_orthonormal().

◆ is_reshape()

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

Definition at line 3647 of file sparsity_internal.cpp.

3647  {
3648  // Quick return if the objects are the same
3649  if (this==&y) return true;
3650 
3651  // Check if same number of entries and nonzeros
3652  if (numel()!=y.numel() || nnz()!=y.nnz()) return false;
3653 
3654  // Quick return if empty interior or dense
3655  if (nnz()==0 || is_dense()) return true;
3656 
3657  // Get Pattern
3658  const casadi_int* colind = this->colind();
3659  const casadi_int* row = this->row();
3660  const casadi_int* y_colind = y.colind();
3661  const casadi_int* y_row = y.row();
3662 
3663  // If same number of rows, check if patterns are identical
3664  if (size1()==y.size1()) return is_equal(y.size1(), y.size2(), y_colind, y_row);
3665 
3666  // Loop over the elements
3667  for (casadi_int cc=0; cc<size2(); ++cc) {
3668  for (casadi_int el=colind[cc]; el<colind[cc+1]; ++el) {
3669  casadi_int rr=row[el];
3670 
3671  // Get row and column of y
3672  casadi_int loc = rr+size1()*cc;
3673  casadi_int rr_y = loc % y.size1();
3674  casadi_int cc_y = loc / y.size1();
3675 
3676  // Make sure matching
3677  if (rr_y != y_row[el] || el<y_colind[cc_y] || el>=y_colind[cc_y+1])
3678  return false;
3679  }
3680  }
3681 
3682  // Reshape if reached this point
3683  return true;
3684  }

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 1797 of file sparsity_internal.cpp.

1797  {
1798  return size1()==1;
1799  }

References size1().

Referenced by is_vector().

◆ is_scalar()

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

Definition at line 1789 of file sparsity_internal.cpp.

1789  {
1790  return size2()==1 && size1()==1 && (!scalar_and_dense || nnz()==1);
1791  }

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 1887 of file sparsity_internal.cpp.

1887  {
1888  return is_orthonormal_rows(allow_empty);
1889  }
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:1863

References is_orthonormal_rows().

◆ is_square()

bool casadi::SparsityInternal::is_square ( ) const

Definition at line 1839 of file sparsity_internal.cpp.

1839  {
1840  return size2() == size1();
1841  }

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 2446 of file sparsity_internal.cpp.

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

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 2340 of file sparsity_internal.cpp.

2340  {
2341  if (is_equal(rhs)) return true;
2342  std::vector<unsigned char> mapping;
2343  shared_from_this<Sparsity>().unite(rhs, mapping);
2344  for (auto e : mapping) {
2345  if (e==1) return false;
2346  }
2347  return true;
2348  }

References is_equal().

◆ is_symmetric()

bool casadi::SparsityInternal::is_symmetric ( ) const

Definition at line 1843 of file sparsity_internal.cpp.

1843  {
1844  return is_transpose(*this);
1845  }
casadi::SparsityInternal::is_transpose
bool is_transpose(const SparsityInternal &y) const
Check if the sparsity is the transpose of another.
Definition: sparsity_internal.cpp:3601

References is_transpose().

Referenced by amd(), and star_coloring_new().

◆ is_transpose()

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

Definition at line 3601 of file sparsity_internal.cpp.

3601  {
3602  // Assert dimensions and number of nonzeros
3603  if (size2()!=y.size1() || size1()!=y.size2() || nnz()!=y.nnz())
3604  return false;
3605 
3606  // Quick return if empty interior or dense
3607  if (nnz()==0 || is_dense())
3608  return true;
3609 
3610  // Run algorithm on the pattern with the least number of rows
3611  if (size1()>size2()) return y.is_transpose(*this);
3612 
3613  // Index counter for columns of the possible transpose
3614  std::vector<casadi_int> y_col_count(y.size2(), 0);
3615  const casadi_int* colind = this->colind();
3616  const casadi_int* row = this->row();
3617  const casadi_int* y_colind = y.colind();
3618  const casadi_int* y_row = y.row();
3619 
3620  // Loop over the columns
3621  for (casadi_int i=0; i<size2(); ++i) {
3622 
3623  // Loop over the nonzeros
3624  for (casadi_int el=colind[i]; el<colind[i+1]; ++el) {
3625 
3626  // Get the row
3627  casadi_int j=row[el];
3628 
3629  // Get the element of the possible transpose
3630  casadi_int el_y = y_colind[j] + y_col_count[j]++;
3631 
3632  // Quick return if element doesn't exist
3633  if (el_y>=y_colind[j+1]) return false;
3634 
3635  // Get the row of the possible transpose
3636  casadi_int j_y = y_row[el_y];
3637 
3638  // Quick return if mismatch
3639  if (j_y != i) return false;
3640  }
3641  }
3642 
3643  // Transpose if reached this point
3644  return true;
3645  }

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 3829 of file sparsity_internal.cpp.

3829  {
3830  const casadi_int* colind = this->colind();
3831  const casadi_int* row = this->row();
3832  // Loop over columns
3833  for (casadi_int i=0; i<size2(); ++i) {
3834  if (colind[i] != colind[i+1]) { // if there are any elements of the column
3835  // check row of the top-most element of the column
3836  casadi_int rr = row[colind[i]];
3837  // Not lower triangular if row > i
3838  if (strictly ? rr <= i : rr < i) return false;
3839  }
3840  }
3841  // all columns ok
3842  return true;
3843  }

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

◆ is_triu()

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

Definition at line 3845 of file sparsity_internal.cpp.

3845  {
3846  const casadi_int* colind = this->colind();
3847  const casadi_int* row = this->row();
3848  // Loop over columns
3849  for (casadi_int i=0; i<size2(); ++i) {
3850  if (colind[i] != colind[i+1]) { // if there are any elements of the column
3851  // check row of the bottom-most element of the column
3852  casadi_int rr = row[colind[i+1]-1];
3853  // Not upper triangular if row>i
3854  if (strictly ? rr >= i : rr > i) return false;
3855  }
3856  }
3857  // all columns ok
3858  return true;
3859  }

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 1805 of file sparsity_internal.cpp.

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

References is_column(), and is_row().

Referenced by get_diag().

◆ largest_first()

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

Definition at line 3511 of file sparsity_internal.cpp.

3511  {
3512  std::vector<casadi_int> degree = get_colind();
3513  casadi_int max_degree = 0;
3514  for (casadi_int k=0; k<size2(); ++k) {
3515  degree[k] = degree[k+1]-degree[k];
3516  max_degree = std::max(max_degree, 1+degree[k]);
3517  }
3518  degree.resize(size2());
3519 
3520  // Vector for binary sort
3521  std::vector<casadi_int> degree_count(max_degree+1, 0);
3522  for (std::vector<casadi_int>::const_iterator it=degree.begin(); it!=degree.end(); ++it) {
3523  degree_count.at(*it+1)++;
3524  }
3525 
3526  // Cumsum to get the offset for each degree
3527  for (casadi_int d=0; d<max_degree; ++d) {
3528  degree_count[d+1] += degree_count[d];
3529  }
3530 
3531  // Now a bucket sort
3532  std::vector<casadi_int> ordering(size2());
3533  for (casadi_int k=size2()-1; k>=0; --k) {
3534  ordering[degree_count[degree[k]]++] = k;
3535  }
3536 
3537  // Invert the ordering
3538  std::vector<casadi_int>& reverse_ordering = degree_count; // reuse memory
3539  reverse_ordering.resize(ordering.size());
3540  std::copy(ordering.begin(), ordering.end(), reverse_ordering.rbegin());
3541 
3542  // Return the ordering
3543  return reverse_ordering;
3544  }

References get_colind(), and size2().

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

◆ 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 378 of file sparsity_internal.cpp.

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

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 413 of file sparsity_internal.cpp.

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

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 2670 of file sparsity_internal.cpp.

2670  {
2671  const casadi_int* colind = this->colind();
2672  const casadi_int* row = this->row();
2673  mapping.resize(nnz());
2674  for (casadi_int i=0; i<size2(); ++i) {
2675  for (casadi_int el=colind[i]; el<colind[i+1]; ++el) {
2676  casadi_int j = row[el];
2677  mapping[el] = j + i*size1();
2678  }
2679  }
2680 
2681  return Sparsity::dense(size1(), size2());
2682  }

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 1040 of file sparsity_internal.cpp.

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

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 1162 of file sparsity_internal.cpp.

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

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 1583 of file sparsity_internal.cpp.

1583  {
1584  /*
1585  Modified version of cs_multiply in CSparse
1586  Copyright(c) Timothy A. Davis, 2006-2009
1587  Licensed as a derivative work under the GNU LGPL
1588  */
1589  casadi_int nz = 0;
1590  casadi_assert(size2() == B.size1(), "Dimension mismatch.");
1591  casadi_int m = size1();
1592  casadi_int anz = nnz();
1593  casadi_int n = B.size2();
1594  const casadi_int* Bp = B.colind();
1595  const casadi_int* Bi = B.row();
1596  casadi_int bnz = Bp[n];
1597 
1598  // get workspace
1599  std::vector<casadi_int> w(m);
1600 
1601  // allocate result
1602  std::vector<casadi_int> C_colind(n+1, 0), C_row;
1603 
1604  C_colind.resize(anz + bnz);
1605 
1606  casadi_int* Cp = &C_colind.front();
1607  for (casadi_int j=0; j<n; ++j) {
1608  if (nz+m > C_row.size()) {
1609  C_row.resize(2*(C_row.size())+m);
1610  }
1611 
1612  // row j of C starts here
1613  Cp[j] = nz;
1614  for (casadi_int p = Bp[j] ; p<Bp[j+1] ; ++p) {
1615  nz = scatter(Bi[p], w, j+1, &C_row.front(), nz);
1616  }
1617  }
1618 
1619  // finalize the last row of C
1620  Cp[n] = nz;
1621  C_row.resize(nz);
1622 
1623  // Success
1624  return Sparsity(m, n, C_colind, C_row);
1625  }
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:1557

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 115 of file sparsity_internal.hpp.

115 { 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_compactible(), 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(), star_coloring_new(), and sub().

◆ nnz_diag()

casadi_int casadi::SparsityInternal::nnz_diag ( ) const

Definition at line 1903 of file sparsity_internal.cpp.

1903  {
1904  const casadi_int* colind = this->colind();
1905  const casadi_int* row = this->row();
1906  casadi_int nnz = 0;
1907  for (casadi_int cc=0; cc<size2(); ++cc) {
1908  for (casadi_int el = colind[cc]; el < colind[cc+1]; ++el) {
1909  nnz += row[el]==cc;
1910  }
1911  }
1912  return nnz;
1913  }

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

◆ nnz_lower()

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

Definition at line 1891 of file sparsity_internal.cpp.

1891  {
1892  const casadi_int* colind = this->colind();
1893  const casadi_int* row = this->row();
1894  casadi_int nnz = 0;
1895  for (casadi_int cc=0; cc<size2(); ++cc) {
1896  for (casadi_int el = colind[cc]; el<colind[cc+1]; ++el) {
1897  if (cc<row[el] || (!strictly && cc==row[el])) nnz++;
1898  }
1899  }
1900  return nnz;
1901  }

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

◆ nnz_upper()

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

Definition at line 1915 of file sparsity_internal.cpp.

1915  {
1916  const casadi_int* colind = this->colind();
1917  const casadi_int* row = this->row();
1918  casadi_int nnz = 0;
1919  for (casadi_int cc=0; cc<size2(); ++cc) {
1920  for (casadi_int el = colind[cc]; el<colind[cc+1]; ++el) {
1921  if (cc>row[el] || (!strictly && cc==row[el])) nnz++;
1922  }
1923  }
1924  return nnz;
1925  }

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

Referenced by star_coloring2().

◆ numel()

casadi_int casadi::SparsityInternal::numel ( ) const

Definition at line 478 of file sparsity_internal.cpp.

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

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 2490 of file sparsity_internal.cpp.

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

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 1464 of file sparsity_internal.cpp.

1465  {
1466  std::vector<casadi_int> colind_C, row_C;
1467  permute(pinv, q, values, colind_C, row_C);
1468  return Sparsity(size1(), size2(), colind_C, row_C);
1469  }

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 1471 of file sparsity_internal.cpp.

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

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 3546 of file sparsity_internal.cpp.

3547  {
3548  // Invert p, possibly
3549  std::vector<casadi_int> p_inv;
3550  if (invert_permutation) {
3551  p_inv.resize(p.size());
3552  for (casadi_int k=0; k<p.size(); ++k) {
3553  p_inv[p[k]] = k;
3554  }
3555  }
3556  const std::vector<casadi_int>& pp = invert_permutation ? p_inv : p;
3557 
3558  // Get columns
3559  std::vector<casadi_int> col = get_col();
3560 
3561  // Get rows
3562  const casadi_int* row = this->row();
3563 
3564  // Sparsity of the return matrix
3565  std::vector<casadi_int> new_row(col.size()), new_col(col.size());
3566 
3567  // Possibly permute columns
3568  if (permute_columns) {
3569  // Assert dimensions
3570  casadi_assert_dev(p.size()==size2());
3571 
3572  // Permute
3573  for (casadi_int k=0; k<col.size(); ++k) {
3574  new_col[k] = pp[col[k]];
3575  }
3576 
3577  } else {
3578  // No permutation of columns
3579  std::copy(col.begin(), col.end(), new_col.begin());
3580  }
3581 
3582  // Possibly permute rows
3583  if (permute_rows) {
3584  // Assert dimensions
3585  casadi_assert_dev(p.size()==size1());
3586 
3587  // Permute
3588  for (casadi_int k=0; k<nnz(); ++k) {
3589  new_row[k] = pp[row[k]];
3590  }
3591 
3592  } else {
3593  // No permutation of rows
3594  std::copy(row, row+nnz(), new_row.begin());
3595  }
3596 
3597  // Return permuted matrix
3598  return Sparsity::triplet(size1(), size2(), new_row, new_col);
3599  }
casadi::SparsityInternal::get_col
std::vector< casadi_int > get_col() const
Get the column for each nonzero.
Definition: sparsity_internal.cpp:491
casadi::invert_permutation
std::vector< casadi_int > invert_permutation(const std::vector< casadi_int > &a)
inverse a permutation vector
Definition: casadi_misc.cpp:182

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 289 of file sparsity_internal.cpp.

292  {
293  // Extract sparsity
294  casadi_int ncol = sp[1];
295  // Calculate elimination tree for A'A
296  etree(sp, parent, w, 1); // len[w] >= nrow+ncol
297  // Calculate postorder
298  casadi_int* post = w; w += ncol;
299  postorder(parent, ncol, post, w); // len[w] >= 3*ncol
300  // Calculate nnz in R
301  *r_nnz = qr_counts(sp_tr, parent, post, w, w+ncol);
302  // Calculate nnz in V
303  *v_nnz = qr_nnz(sp, pinv, leftmost, parent, nrow_ext, w);
304  }
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  pa = parent[c];
276  if (pa != -1) {
277  // Move all rows to parent of c
278  if (nque[pa]==0) tail[pa] = tail[c];
279  next[tail[c]] = head[pa];
280  head[pa] = next[r];
281  nque[pa] += nque[c];
282  }
283  }
284  for (r=0; r<nrow; ++r) if (pinv[r]<0) pinv[r] = c++;
285  if (nrow_ext) *nrow_ext = nrow_new;
286  return v_nnz;
287  }

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 306 of file sparsity_internal.cpp.

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

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 1416 of file sparsity_internal.cpp.

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

Referenced by maxtrans().

◆ repr_el()

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

Definition at line 1717 of file sparsity_internal.cpp.

1717  {
1718  casadi_int start_index = GlobalOptions::start_index;
1719  std::stringstream ss;
1720  if (numel()!=nnz()) {
1721  ss << "nonzero index " << k+start_index << " ";
1722  }
1723  casadi_int r = row()[k];
1724  casadi_int c = get_col()[k];
1725  ss << "(row " << r+start_index << ", col " << c+start_index << ")";
1726 
1727  return ss.str();
1728  }
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 112 of file sparsity_internal.hpp.

112 { 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_compactible(), 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(), star_coloring_new(), sub(), and uni_coloring().

◆ rowsSequential()

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

Definition at line 2773 of file sparsity_internal.cpp.

2773  {
2774  const casadi_int* colind = this->colind();
2775  const casadi_int* row = this->row();
2776  for (casadi_int i=0; i<size2(); ++i) {
2777  casadi_int lastrow = -1;
2778  for (casadi_int k=colind[i]; k<colind[i+1]; ++k) {
2779 
2780  // check if not in sequence
2781  if (row[k] < lastrow)
2782  return false;
2783 
2784  // Check if duplicate
2785  if (strictly && row[k] == lastrow)
2786  return false;
2787 
2788  // update last row of the col
2789  lastrow = row[k];
2790  }
2791  }
2792 
2793  // sequential if reached this point
2794  return true;
2795  }

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 1078 of file sparsity_internal.cpp.

1078  {
1079  /*
1080  Modified version of cs_rprune in CSparse
1081  Copyright(c) Timothy A. Davis, 2006-2009
1082  Licensed as a derivative work under the GNU LGPL
1083  */
1084  std::vector<casadi_int> &rr = *static_cast<std::vector<casadi_int> *>(other);
1085  return (i >= rr[1] && i < rr[2]) ;
1086  }

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 1557 of file sparsity_internal.cpp.

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

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 584 of file sparsity_internal.cpp.

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

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 1927 of file sparsity_internal.cpp.

1927  {
1928  return std::pair<casadi_int, casadi_int>(size1(), size2());
1929  }

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 97 of file sparsity_internal.hpp.

97 { 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(), star_coloring_new(), sub(), and transpose().

◆ size2()

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

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

Definition at line 102 of file sparsity_internal.hpp.

102 { 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 92 of file sparsity_internal.hpp.

92 { 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 3959 of file sparsity_internal.cpp.

3960  {
3961  const Btf& btf = this->btf();
3962  const casadi_int* colind = this->colind();
3963  const casadi_int* row = this->row();
3964 
3965  if (!tr) {
3966  for (casadi_int b = 0; b < btf.nb; ++b) { // loop over the blocks forward
3967 
3968  // Get dependencies from all right-hand-sides in the block ...
3969  bvec_t block_dep = 0;
3970  for (casadi_int el=btf.rowblock[b]; el<btf.rowblock[b+1]; ++el) {
3971  casadi_int rr = btf.rowperm[el];
3972  block_dep |= B[rr];
3973  }
3974 
3975  // ... as well as all other variables in the block
3976  for (casadi_int el=btf.colblock[b]; el<btf.colblock[b+1]; ++el) {
3977  casadi_int cc = btf.colperm[el];
3978  block_dep |= X[cc];
3979  }
3980 
3981  // Propagate ...
3982  for (casadi_int el=btf.colblock[b]; el<btf.colblock[b+1]; ++el) {
3983  casadi_int cc = btf.colperm[el];
3984 
3985  // ... to all variables in the block ...
3986  X[cc] |= block_dep;
3987 
3988  // ... as well as to all right-hand-sides that depend on variables in the block
3989  for (casadi_int k=colind[cc]; k<colind[cc+1]; ++k) {
3990  casadi_int rr=row[k];
3991  B[rr] |= block_dep;
3992  }
3993  }
3994  }
3995 
3996  } else { // transpose
3997  for (casadi_int b = btf.nb; b-- > 0; ) { // loop over the blocks backward
3998 
3999  // Get dependencies ...
4000  bvec_t block_dep = 0;
4001  for (casadi_int el=btf.colblock[b]; el<btf.colblock[b+1]; ++el) {
4002  casadi_int cc = btf.colperm[el];
4003 
4004  // .. from all right-hand-sides in the block ...
4005  block_dep |= B[cc];
4006 
4007  // ... as well as from all depending variables ...
4008  for (casadi_int k=colind[cc]; k<colind[cc+1]; ++k) {
4009  casadi_int rr=row[k];
4010  block_dep |= X[rr];
4011  }
4012  }
4013 
4014  // Propagate to all right-hand-sides in the block ...
4015  for (casadi_int el=btf.colblock[b]; el<btf.colblock[b+1]; ++el) {
4016  casadi_int cc = btf.colperm[el];
4017  B[cc] |= block_dep;
4018  }
4019 
4020  // ... as well as all other variables in the block
4021  for (casadi_int el=btf.rowblock[b]; el<btf.rowblock[b+1]; ++el) {
4022  casadi_int rr = btf.rowperm[el];
4023  X[rr] |= block_dep;
4024  }
4025  }
4026  }
4027  }
casadi::bvec_t
unsigned long long bvec_t
Definition: casadi_common.hpp:77

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 3697 of file sparsity_internal.cpp.

3697  {
3698 
3699  // Index counter for each column
3700  std::vector<casadi_int> cind = get_colind();
3701  const casadi_int* colind = this->colind();
3702  const casadi_int* row = this->row();
3703 
3704  // Loop over rows
3705  for (casadi_int rr=0; rr<size1(); ++rr) {
3706 
3707  // Loop over columns
3708  for (casadi_int cc=0; cc<size2(); ++cc) {
3709  // Check if nonzero
3710  if (cind[cc]<colind[cc+1] && row[cind[cc]]==rr) {
3711  stream << "*";
3712  cind[cc]++;
3713  } else {
3714  stream << ".";
3715  }
3716  }
3717 
3718  // End of row
3719  stream << std::endl;
3720  }
3721  }

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

◆ spy_matlab()

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

Definition at line 3809 of file sparsity_internal.cpp.

3809  {
3810  // Create the .m file
3811  auto mfile_ptr = Filesystem::ofstream_ptr(mfile_name);
3812  std::ostream& mfile = *mfile_ptr;
3813 
3814  // Header
3815  mfile << "% This function was automatically generated by CasADi" << std::endl;
3816 
3817  Dict opts;
3818  opts["name"] = "A";
3819  export_code("matlab", mfile, opts);
3820 
3821  // Issue spy command
3822  mfile << "spy(A);" << std::endl;
3823  }
casadi::Filesystem::ofstream_ptr
static std::unique_ptr< std::ostream > ofstream_ptr(const std::string &path, std::ios_base::openmode mode=std::ios_base::out)
Definition: filesystem.cpp:115
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:3723
casadi::Dict
GenericType::Dict Dict
C++ equivalent of Python's dict or MATLAB's struct.
Definition: generic_type.hpp:258

References export_code(), and casadi::Filesystem::ofstream_ptr().

◆ 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 3398 of file sparsity_internal.cpp.

3398  {
3399  if (!is_square()) {
3400  // NOTE(@jaeandersson) Why warning and not error?
3401  casadi_message("StarColoring requires a square matrix, got " + dim() + ".");
3402  }
3403 
3404  // Reorder, if necessary
3405  if (ordering!=0) {
3406  casadi_assert_dev(ordering==1);
3407 
3408  // Ordering
3409  std::vector<casadi_int> ord = largest_first();
3410 
3411  // Create a new sparsity pattern
3412  Sparsity sp_permuted = pmult(ord, true, true, true);
3413 
3414  // Star coloring for the permuted matrix
3415  Sparsity ret_permuted = sp_permuted.star_coloring(0);
3416 
3417  // Permute result back
3418  return ret_permuted.pmult(ord, true, false, false);
3419  }
3420 
3421  // Allocate temporary vectors
3422  const casadi_int* colind = this->colind();
3423  const casadi_int* row = this->row();
3424  std::vector<casadi_int> forbiddenColors;
3425  forbiddenColors.reserve(size2());
3426  std::vector<casadi_int> color(size2(), -1);
3427 
3428  // 4: for i <- 1 to |V | do
3429  for (casadi_int i=0; i<size2(); ++i) {
3430 
3431  // 5: for each w \in N1 (vi) do
3432  for (casadi_int w_el=colind[i]; w_el<colind[i+1]; ++w_el) {
3433  casadi_int w = row[w_el];
3434 
3435  // 6: if w is colored then
3436  if (color[w]!=-1) {
3437 
3438  // 7: forbiddenColors[color[w]] <- v
3439  forbiddenColors[color[w]] = i;
3440 
3441  } // 8: end if
3442 
3443  // 9: for each colored vertex x \in N1 (w) do
3444  for (casadi_int x_el=colind[w]; x_el<colind[w+1]; ++x_el) {
3445  casadi_int x = row[x_el];
3446  if (color[x]==-1) continue;
3447 
3448  // 10: if w is not colored then
3449  if (color[w]==-1) {
3450 
3451  //11: forbiddenColors[color[x]] <- vi
3452  forbiddenColors[color[x]] = i;
3453 
3454  } else { // 12: else
3455 
3456  // 13: for each colored vertex y \in N1 (x), y != w do
3457  for (casadi_int y_el=colind[x]; y_el<colind[x+1]; ++y_el) {
3458  casadi_int y = row[y_el];
3459  if (color[y]==-1 || y==w) continue;
3460 
3461  // 14: if color[y] = color[w] then
3462  if (color[y]==color[w]) {
3463 
3464  // 15: forbiddenColors[color[x]] <- vi
3465  forbiddenColors[color[x]] = i;
3466 
3467  // 16: break
3468  break;
3469 
3470  } // 17: end if
3471 
3472  } // 18: end for
3473 
3474  } // 19: end if
3475 
3476  } // 20 end for
3477 
3478  } // 21 end for
3479 
3480  // 22: color[v] <- min {c > 0 : forbiddenColors[c] = v}
3481  bool new_color = true;
3482  for (casadi_int color_i=0; color_i<forbiddenColors.size(); ++color_i) {
3483  // Break if color is ok
3484  if (forbiddenColors[color_i]!=i) {
3485  color[i] = color_i;
3486  new_color = false;
3487  break;
3488  }
3489  }
3490 
3491  // New color if reached end
3492  if (new_color) {
3493  color[i] = forbiddenColors.size();
3494  forbiddenColors.push_back(-1);
3495 
3496  // Cutoff if too many colors
3497  if (forbiddenColors.size()>cutoff) {
3498  return Sparsity();
3499  }
3500  }
3501 
3502  } // 23 end for
3503 
3504  // Number of colors used
3505  casadi_int num_colors = forbiddenColors.size();
3506 
3507  // Return sparsity in sparse triplet format
3508  return Sparsity::triplet(size2(), num_colors, range(color.size()), color);
3509  }
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:3546
casadi::SparsityInternal::largest_first
std::vector< casadi_int > largest_first() const
Order the columns by decreasing degree.
Definition: sparsity_internal.cpp:3511

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(), and star_coloring_new().

◆ 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 3048 of file sparsity_internal.cpp.

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

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(), and star_coloring_new().

◆ star_coloring_new()

Sparsity casadi::SparsityInternal::star_coloring_new ( std::vector< casadi_int > &  which_color,
const Dict &  opts 
) const

See description in public class.

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

Definition at line 3301 of file sparsity_internal.cpp.

3302  {
3303  // Read options
3304  bool new_algo = false;
3305  casadi_int cutoff = std::numeric_limits<casadi_int>::max();
3306  bool largest_first = true;
3307  for (auto&& e : opts) {
3308  if (e.first=="new_algo") {
3309  new_algo = e.second;
3310  } else if (e.first=="cutoff") {
3311  cutoff = e.second;
3312  } else if (e.first=="largest_first") {
3313  largest_first = e.second;
3314  } else {
3315  casadi_warning("Unknown option '" + e.first + "' for star_coloring_new.");
3316  }
3317  }
3318 
3319  // Consistency check
3320  casadi_assert(is_symmetric(), "Expected symmetric matrix");
3321 
3322  // Call either of the algorithms
3323  Sparsity coloring;
3324  if (new_algo) {
3325  coloring = star_coloring2(largest_first ? 1 : 0, cutoff);
3326  } else {
3327  coloring = star_coloring(largest_first ? 1 : 0, cutoff);
3328  }
3329 
3330  // Get the coloring
3331  casadi_int ncolor = coloring.size2();
3332  const casadi_int *col_colind = coloring.colind(), *col_row = coloring.row();
3333 
3334  // Get sparsity pattern of the original matrix
3335  casadi_int nrow = size1();
3336  const casadi_int* colind = this->colind();
3337  const casadi_int* row = this->row();
3338 
3339  // How many of the unit vectors influence each output
3340  std::vector<casadi_int> output_ctr(nrow);
3341 
3342  // Loop over colors to find which color to use to calculate each nonzero
3343  which_color.resize(nnz());
3344  std::fill(which_color.begin(), which_color.end(), -1);
3345  for (casadi_int color=0; color<ncolor; ++color) {
3346  // How many times does each output get perturbed for the current color
3347  std::fill(output_ctr.begin(), output_ctr.end(), 0);
3348  for (casadi_int el=col_colind[color]; el<col_colind[color+1]; ++el) {
3349  casadi_int c = col_row[el];
3350  for (casadi_int k=colind[c]; k<colind[c+1]; ++k) output_ctr[row[k]]++;
3351  }
3352 
3353  // Loop over corresponding columns in the original matrix
3354  for (casadi_int el=col_colind[color]; el<col_colind[color+1]; ++el) {
3355  casadi_int c = col_row[el];
3356  for (casadi_int k=colind[c]; k<colind[c+1]; ++k) {
3357  // Set which_color, if output impacted by exactly one unit vector
3358  if (output_ctr[row[k]] == 1 && which_color[k] < 0) which_color[k] = color;
3359  }
3360  }
3361  }
3362 
3363  // Select a color for each nonzero, taking into account symmetry
3364  std::vector<casadi_int> rowind(colind, colind + nrow);
3365  for (casadi_int c=0; c<nrow; ++c) {
3366  for (casadi_int k=colind[c]; k<colind[c+1]; ++k) {
3367  casadi_int r = row[k];
3368  if (r > c) {
3369  // Lower triangular element, get corresponding upper triangular element
3370  casadi_int k_tr = rowind[r]++;
3371  // Pick the color with the smallest index
3372  if (which_color[k] < 0) {
3373  // Pick upper triangular element
3374  casadi_assert(which_color[k_tr] >= 0,
3375  "Expected color for element " + str(c) + ", " + str(r));
3376  } else if (which_color[k_tr] < 0) {
3377  // Pick lower triangular element
3378  casadi_assert(which_color[k] >= 0,
3379  "Expected color for element " + str(r) + ", " + str(c));
3380  } else if (which_color[k] <= which_color[k_tr]) {
3381  // Both can be calculated, but lower triangular has smaller index
3382  which_color[k_tr] = -1;
3383  } else {
3384  // Both can be calculated, but upper triangular has smaller index
3385  which_color[k] = -1;
3386  }
3387  } else if (r == c) {
3388  // Ensure diagonal elements are colored
3389  casadi_assert(which_color[k] >= 0, "No color for diagonal element " + str(r));
3390  }
3391  }
3392  }
3393 
3394  // Return coloring (in addition to which_color mapping)
3395  return coloring;
3396  }
casadi::SparsityInternal::star_coloring2
Sparsity star_coloring2(casadi_int ordering, casadi_int cutoff) const
An improved distance-2 coloring algorithm.
Definition: sparsity_internal.cpp:3048
casadi::SparsityInternal::star_coloring
Sparsity star_coloring(casadi_int ordering, casadi_int cutoff) const
A greedy distance-2 coloring algorithm.
Definition: sparsity_internal.cpp:3398

References casadi::Sparsity::colind(), colind(), is_symmetric(), largest_first(), nnz(), casadi::Sparsity::row(), row(), size1(), casadi::Sparsity::size2(), star_coloring(), star_coloring2(), and casadi::str().

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

◆ 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 2160 of file sparsity_internal.cpp.

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

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 2204 of file sparsity_internal.cpp.

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

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 502 of file sparsity_internal.cpp.

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

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 509 of file sparsity_internal.cpp.

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

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 2959 of file sparsity_internal.cpp.

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

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 1063 of file sparsity_internal.cpp.

1064  {
1065  /*
1066  Modified version of cs_unmatched in CSparse
1067  Copyright(c) Timothy A. Davis, 2006-2009
1068  Licensed as a derivative work under the GNU LGPL
1069  */
1070  casadi_int i, kr = rr[set] ;
1071  for (i=0; i<m; i++)
1072  if (wi[i] == 0)
1073  p[kr++] = i;
1074 
1075  rr[set+1] = kr;
1076  }

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 1533 of file sparsity_internal.cpp.

1534  {
1535  /*
1536  Modified version of cs_wclear in CSparse
1537  Copyright(c) Timothy A. Davis, 2006-2009
1538  Licensed as a derivative work under the GNU LGPL
1539  */
1540  if (mark < 2 || (mark + lemax < 0)) {
1541  for (casadi_int k = 0; k<n; ++k) if (w[k] != 0) w[k] = 1;
1542  mark = 2 ;
1543  }
1544  // at this point, w [0..n-1] < mark holds
1545  return mark;
1546  }

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.

210  {
211 #ifdef CASADI_WITH_THREADSAFE_SYMBOLICS
212  auto* w = weak_ref_.load(std::memory_order_acquire);
213  if (!w) {
214  auto* nw = new GenericWeakRef<Shared, Internal>(static_cast<Internal*>(this));
215  GenericWeakRef<Shared, Internal>* expected = nullptr;
216  if (weak_ref_.compare_exchange_strong(
217  expected, nw, std::memory_order_release, std::memory_order_acquire)) {
218  w = nw;
219  } else {
220  delete nw; // lost the race; another thread published first
221  w = expected;
222  }
223  }
224  return w;
225 #else
226  if (weak_ref_==nullptr) {
227  weak_ref_ = new GenericWeakRef<Shared, Internal>(static_cast<Internal*>(this));
228  }
229  return weak_ref_;
230 #endif // CASADI_WITH_THREADSAFE_SYMBOLICS
231  }
The documentation for this class was generated from the following files: