sparsity.hpp
1 /*
2  * This file is part of CasADi.
3  *
4  * CasADi -- A symbolic framework for dynamic optimization.
5  * Copyright (C) 2010-2023 Joel Andersson, Joris Gillis, Moritz Diehl,
6  * KU Leuven. All rights reserved.
7  * Copyright (C) 2011-2014 Greg Horn
8  * Copyright (C) 2005-2013 Timothy A. Davis
9  *
10  * CasADi is free software; you can redistribute it and/or
11  * modify it under the terms of the GNU Lesser General Public
12  * License as published by the Free Software Foundation; either
13  * version 3 of the License, or (at your option) any later version.
14  *
15  * CasADi is distributed in the hope that it will be useful,
16  * but WITHOUT ANY WARRANTY; without even the implied warranty of
17  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18  * Lesser General Public License for more details.
19  *
20  * You should have received a copy of the GNU Lesser General Public
21  * License along with CasADi; if not, write to the Free Software
22  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
23  *
24  */
25 
26 
27 #ifndef CASADI_SPARSITY_HPP
28 #define CASADI_SPARSITY_HPP
29 
30 #include "shared_object.hpp"
31 #include "printable.hpp"
32 #include "casadi_common.hpp"
33 #include "sparsity_interface.hpp"
34 #include "generic_type.hpp"
35 #include <vector>
36 #include <list>
37 #include <limits>
38 #include <unordered_map>
39 #ifdef CASADI_WITH_THREAD
40 #ifdef CASADI_WITH_THREAD_MINGW
41 #include <mingw.mutex.h>
42 #else // CASADI_WITH_THREAD_MINGW
43 #include <mutex>
44 #endif // CASADI_WITH_THREAD_MINGW
45 #endif //CASADI_WITH_THREAD
46 
47 namespace casadi {
48  // Forward declaration
49  class SparsityInternal;
50  class SerializingStream;
51  class DeserializingStream;
52 
53  #ifndef SWIG
57  struct CASADI_EXPORT SparsityStruct {
58  casadi_int nrow;
59  casadi_int ncol;
60  const casadi_int* colind;
61  const casadi_int* row;
62  };
63  #endif // SWIG
64 
197  class CASADI_EXPORT Sparsity
198  : public SharedObject,
199  public SWIG_IF_ELSE(SparsityInterfaceCommon, SparsityInterface<Sparsity>),
200  public SWIG_IF_ELSE(PrintableCommon, Printable<Sparsity>) {
201  public:
202 
248  explicit Sparsity(casadi_int dummy=0);
249 
303  Sparsity(casadi_int nrow, casadi_int ncol);
304 
379  Sparsity(casadi_int nrow, casadi_int ncol,
380  const std::vector<casadi_int>& colind, const std::vector<casadi_int>& row,
381  bool order_rows=false);
382 
386  explicit Sparsity(const std::pair<casadi_int, casadi_int>& rc);
387 
388 #ifndef SWIG
390  Sparsity(casadi_int nrow, casadi_int ncol, const casadi_int* colind, const casadi_int* row,
391  bool order_rows=false);
392 
396  static Sparsity create(SparsityInternal *node);
397 
399  typedef SparsityInterface<Sparsity> B;
400 
402  using B::horzsplit;
403  using B::diagsplit;
404  using B::vertsplit;
405  using B::mtimes;
406 
407  SparsityInternal* get() const;
408 #endif
409 
454  static Sparsity scalar(bool dense_scalar=true)
455  { return dense_scalar ? dense(1, 1) : Sparsity(1, 1); }
457 
520  static Sparsity dense(casadi_int nrow, casadi_int ncol=1);
521  static Sparsity dense(const std::pair<casadi_int, casadi_int> &rc) {
522  return dense(rc.first, rc.second);
523  }
525 
563  static Sparsity unit(casadi_int n, casadi_int el);
565 
612  static Sparsity upper(casadi_int n);
613 
664  static Sparsity lower(casadi_int n);
665 
723  static Sparsity diag(casadi_int nrow) { return diag(nrow, nrow);}
724  static Sparsity diag(casadi_int nrow, casadi_int ncol);
725  static Sparsity diag(const std::pair<casadi_int, casadi_int> &rc) {
726  return diag(rc.first, rc.second);
727  }
729 
797  static Sparsity band(casadi_int n, casadi_int p);
798 
849  static Sparsity banded(casadi_int n, casadi_int p);
850 
899  static Sparsity rowcol(const std::vector<casadi_int>& row,
900  const std::vector<casadi_int>& col,
901  casadi_int nrow, casadi_int ncol);
902 
904 
992  static Sparsity triplet(casadi_int nrow, casadi_int ncol,
993  const std::vector<casadi_int>& row, const std::vector<casadi_int>& col,
994  std::vector<casadi_int>& SWIG_OUTPUT(mapping), bool invert_mapping);
995  static Sparsity triplet(casadi_int nrow, casadi_int ncol, const std::vector<casadi_int>& row,
996  const std::vector<casadi_int>& col);
998 
1057  static Sparsity nonzeros(casadi_int nrow, casadi_int ncol, const std::vector<casadi_int>& nz,
1058  bool ind1=SWIG_IND1);
1059 
1119  static Sparsity compressed(const std::vector<casadi_int>& v, bool order_rows=false);
1120 #ifndef SWIG
1121  static Sparsity compressed(const casadi_int* v, bool order_rows=false);
1122 #endif // SWIG
1124 
1200  static Sparsity permutation(const std::vector<casadi_int>& p, bool invert=false);
1201 
1256  const std::vector<casadi_int> permutation_vector(bool invert=false) const;
1257 
1315  Sparsity get_diag(std::vector<casadi_int>& SWIG_OUTPUT(mapping)) const;
1316 
1366  std::vector<casadi_int> compress(bool canonical=true) const;
1367 
1368 #ifndef SWIG
1370  const SparsityInternal* operator->() const;
1371 
1373  const SparsityInternal& operator*() const;
1374 #endif // SWIG
1430  bool is_equal(const Sparsity& y) const;
1432  bool is_equal(casadi_int nrow, casadi_int ncol, const std::vector<casadi_int>& colind,
1433  const std::vector<casadi_int>& row) const;
1434 #ifndef SWIG
1435  bool is_equal(casadi_int nrow, casadi_int ncol,
1436  const casadi_int* colind, const casadi_int* row) const;
1437 #endif // SWIG
1438 
1439  bool operator==(const Sparsity& y) const { return is_equal(y);}
1441 
1443  bool operator!=(const Sparsity& y) const {return !is_equal(y);}
1444 
1486  bool is_stacked(const Sparsity& y, casadi_int n) const;
1487 
1488 #ifndef SWIG
1495  operator const casadi_int*() const;
1496 
1500  operator const std::vector<casadi_int>&() const;
1501 
1505  operator SparsityStruct() const;
1506 #endif // SWIG
1507 
1510 
1540  casadi_int size1() const;
1541 
1543  casadi_int rows() const {return size1();}
1544 
1546  casadi_int size2() const;
1547 
1549  casadi_int columns() const {return size2();}
1550 
1582  casadi_int numel() const;
1583 
1619  double density() const;
1620 
1666  bool is_empty(bool both=false) const;
1667 
1676  casadi_int nnz() const;
1677 
1732  casadi_int nnz_upper(bool strictly=false) const;
1733 
1741  casadi_int nnz_lower(bool strictly=false) const;
1742 
1748  casadi_int nnz_diag() const;
1749 
1790  casadi_int bw_upper() const;
1791 
1832  casadi_int bw_lower() const;
1833 
1867  std::pair<casadi_int, casadi_int> size() const;
1868 
1875  casadi_int size(casadi_int axis) const;
1877 
1919  Dict info() const;
1920 
1998  void to_file(const std::string& filename, const std::string& format_hint="") const;
1999 
2008  static Sparsity from_file(const std::string& filename, const std::string& format_hint="");
2009 
2010 #ifndef SWIG
2014  void serialize(std::ostream &stream) const;
2015 #endif
2016 
2060  std::string serialize() const;
2061 
2065  static Sparsity deserialize(std::istream& stream);
2066 
2072  static Sparsity deserialize(const std::string& s);
2073 
2077  void serialize(SerializingStream& s) const;
2078 
2082  static Sparsity deserialize(DeserializingStream& s);
2083 
2084 #ifndef SWIG
2090  const casadi_int* row() const;
2091 
2095  const casadi_int* colind() const;
2096 #endif
2097 
2156  std::vector<casadi_int> get_row() const;
2157 
2167  std::vector<casadi_int> get_colind() const;
2168 
2212  casadi_int colind(casadi_int cc) const;
2213 
2221  casadi_int row(casadi_int el) const;
2222 
2231  std::vector<casadi_int> get_col() const;
2232 
2289  void resize(casadi_int nrow, casadi_int ncol);
2290 
2351  casadi_int add_nz(casadi_int rr, casadi_int cc);
2352 
2397  casadi_int get_nz(casadi_int rr, casadi_int cc) const;
2398 
2430  bool has_nz(casadi_int rr, casadi_int cc) const;
2431 
2465  std::vector<casadi_int> get_nz(const std::vector<casadi_int>& rr,
2466  const std::vector<casadi_int>& cc) const;
2467 
2514  void get_nz(std::vector<casadi_int>& SWIG_INOUT(indices)) const;
2515 
2553  std::vector<casadi_int> get_lower() const;
2554 
2560  std::vector<casadi_int> get_upper() const;
2561 
2594  void get_ccs(std::vector<casadi_int>& SWIG_OUTPUT(colind),
2595  std::vector<casadi_int>& SWIG_OUTPUT(row)) const;
2596 
2633  void get_crs(std::vector<casadi_int>& SWIG_OUTPUT(rowind),
2634  std::vector<casadi_int>& SWIG_OUTPUT(col)) const;
2635 
2670  void get_triplet(std::vector<casadi_int>& SWIG_OUTPUT(row),
2671  std::vector<casadi_int>& SWIG_OUTPUT(col)) const;
2672 
2745  Sparsity sub(const std::vector<casadi_int>& rr,
2746  const std::vector<casadi_int>& cc,
2747  std::vector<casadi_int>& SWIG_OUTPUT(mapping), bool ind1=false) const;
2748 
2815  Sparsity sub(const std::vector<casadi_int>& rr, const Sparsity& sp,
2816  std::vector<casadi_int>& SWIG_OUTPUT(mapping), bool ind1=false) const;
2817 
2875  Sparsity T() const;
2876 
2914  Sparsity transpose(std::vector<casadi_int>& SWIG_OUTPUT(mapping),
2915  bool invert_mapping=false) const;
2916 
2955  bool is_transpose(const Sparsity& y) const;
2956 
2998  bool is_reshape(const Sparsity& y) const;
2999 
3046  bool is_compactible(std::vector<casadi_int>& SWIG_OUTPUT(row),
3047  std::vector<casadi_int>& SWIG_OUTPUT(col)) const;
3048 
3050 
3121 #ifndef SWIG
3122  Sparsity combine(const Sparsity& y, bool f0x_is_zero, bool function0_is_zero,
3123  std::vector<unsigned char>& mapping) const;
3124 #endif // SWIG
3125  Sparsity combine(const Sparsity& y, bool f0x_is_zero, bool function0_is_zero) const;
3127 
3129 
3185 #ifndef SWIG
3186  Sparsity unite(const Sparsity& y, std::vector<unsigned char>& mapping) const;
3187 #endif // SWIG
3188  Sparsity unite(const Sparsity& y) const;
3189  Sparsity operator+(const Sparsity& b) const;
3191 
3193 
3234 #ifndef SWIG
3235  Sparsity intersect(const Sparsity& y,
3236  std::vector<unsigned char>& mapping) const;
3237 #endif // SWIG
3238  Sparsity intersect(const Sparsity& y) const;
3239  Sparsity operator*(const Sparsity& b) const;
3241 
3289  bool is_subset(const Sparsity& rhs) const;
3290 
3336  Sparsity sparsity_cast_mod(const Sparsity& X, const Sparsity& Y) const;
3337 
3379  Sparsity pattern_inverse() const;
3380 
3381 #ifndef SWIG
3388  static void mul_sparsityF(const bvec_t* x, const Sparsity& x_sp,
3389  const bvec_t* y, const Sparsity& y_sp,
3390  bvec_t* z, const Sparsity& z_sp,
3391  bvec_t* w);
3392 
3400  static void mul_activityF(const bvec_t* x, const Sparsity& x_sp,
3401  const bvec_t* y, const Sparsity& y_sp,
3402  bvec_t* z, const Sparsity& z_sp,
3403  bvec_t* w);
3404 
3411  static void mul_sparsityR(bvec_t* x, const Sparsity& x_sp,
3412  bvec_t* y, const Sparsity& y_sp,
3413  bvec_t* z, const Sparsity& z_sp,
3414  bvec_t* w);
3415 
3418 
3421  static Sparsity horzcat(const std::vector<Sparsity> & sp);
3422  static Sparsity vertcat(const std::vector<Sparsity> & sp);
3423  static Sparsity blockcat(const std::vector< std::vector< Sparsity > > &v);
3424  static Sparsity diagcat(const std::vector< Sparsity > &v);
3425  static std::vector<Sparsity>
3426  horzsplit(const Sparsity& x, const std::vector<casadi_int>& offset);
3427  static std::vector<Sparsity>
3428  vertsplit(const Sparsity& x, const std::vector<casadi_int>& offset);
3429  static std::vector<Sparsity>
3430  diagsplit(const Sparsity& x,
3431  const std::vector<casadi_int>& offset1,
3432  const std::vector<casadi_int>& offset2);
3433  static Sparsity mtimes(const Sparsity& x, const Sparsity& y,
3434  const std::string& blas = "reference");
3435  static Sparsity mac(const Sparsity& x, const Sparsity& y, const Sparsity& z,
3436  const std::string& /*blas*/ = "reference") { return z;}
3437  static Sparsity reshape(const Sparsity& x, casadi_int nrow, casadi_int ncol);
3438  static Sparsity reshape(const Sparsity& x, const Sparsity& sp);
3439  static Sparsity sparsity_cast(const Sparsity& x, const Sparsity& sp);
3440  static casadi_int sprank(const Sparsity& x);
3441  static casadi_int norm_0_mul(const Sparsity& x, const Sparsity& A);
3442  static Sparsity kron(const Sparsity& a, const Sparsity& b);
3443 
3459  static Sparsity kron_contract(const Sparsity& sp_m, const Sparsity& sp_x, bool inner);
3460  static Sparsity triu(const Sparsity& x, bool includeDiagonal=true);
3461  static Sparsity tril(const Sparsity& x, bool includeDiagonal=true);
3462 
3463  static Sparsity sum2(const Sparsity &x);
3464  static Sparsity sum1(const Sparsity &x);
3467 #endif //SWIG
3468 
3538  void enlarge(casadi_int nrow, casadi_int ncol,
3539  const std::vector<casadi_int>& rr,
3540  const std::vector<casadi_int>& cc, bool ind1=false);
3541 
3581  void enlargeRows(casadi_int nrow, const std::vector<casadi_int>& rr, bool ind1=false);
3582 
3619  void enlargeColumns(casadi_int ncol, const std::vector<casadi_int>& cc, bool ind1=false);
3620 
3655  Sparsity makeDense(std::vector<casadi_int>& SWIG_OUTPUT(mapping)) const;
3656 
3711  std::vector<casadi_int> erase(const std::vector<casadi_int>& rr,
3712  const std::vector<casadi_int>& cc, bool ind1=false);
3713 
3777  std::vector<casadi_int> erase(const std::vector<casadi_int>& rr, bool ind1=false);
3778 
3836  void append(const Sparsity& sp);
3837 
3882  void appendColumns(const Sparsity& sp);
3883 
3919  bool is_scalar(bool scalar_and_dense=false) const;
3920 
3962  bool is_dense() const;
3963 
3998  bool is_row() const;
3999 
4032  bool is_column() const;
4033 
4074  bool is_vector() const;
4075 
4117  bool is_diag() const;
4118 
4149  bool is_square() const;
4150 
4186  bool is_symmetric() const;
4187 
4230  bool is_triu(bool strictly = false) const;
4231 
4272  bool is_tril(bool strictly = false) const;
4273 
4323  bool is_singular() const;
4324 
4366  bool is_permutation() const;
4367 
4422  bool is_selection(bool allow_empty=false) const;
4423 
4481  bool is_orthonormal(bool allow_empty=false) const;
4482 
4524  bool is_orthonormal_rows(bool allow_empty=false) const;
4525 
4537  bool is_orthonormal_columns(bool allow_empty=false) const;
4538 
4576  bool rowsSequential(bool strictly=true) const;
4577 
4617  void removeDuplicates(std::vector<casadi_int>& SWIG_INOUT(mapping));
4618 
4619 #ifndef SWIG
4620  typedef std::unordered_multimap<std::size_t, WeakRef> CachingMap;
4621 
4623  static CachingMap& getCache();
4624 
4625 #ifdef CASADI_WITH_THREADSAFE_SYMBOLICS
4626  // Safe access to CachingMap
4627  static std::mutex cachingmap_mtx;
4628 #endif //CASADI_WITH_THREADSAFE_SYMBOLICS
4629 
4631  static const Sparsity& getScalar();
4632 
4634  static const Sparsity& getScalarSparse();
4635 
4637  static const Sparsity& getEmpty();
4638 
4639 #endif //SWIG
4640 
4684  std::vector<casadi_int> etree(bool ata=false) const;
4685 
4757  Sparsity ldl(std::vector<casadi_int>& SWIG_OUTPUT(p), bool amd=true) const;
4758 
4817  void qr_sparse(Sparsity& SWIG_OUTPUT(V), Sparsity& SWIG_OUTPUT(R),
4818  std::vector<casadi_int>& SWIG_OUTPUT(prinv),
4819  std::vector<casadi_int>& SWIG_OUTPUT(pc), bool amd=true) const;
4820 
4866  casadi_int dfs(casadi_int j, casadi_int top, std::vector<casadi_int>& SWIG_INOUT1(xi),
4867  std::vector<casadi_int>& SWIG_INOUT2(pstack),
4868  const std::vector<casadi_int>& pinv, std::vector<bool>& SWIG_INOUT3(marked)) const;
4869 
4963  casadi_int scc(std::vector<casadi_int>& SWIG_OUTPUT(index),
4964  std::vector<casadi_int>& SWIG_OUTPUT(offset)) const;
4965 
5059  casadi_int btf(std::vector<casadi_int>& SWIG_OUTPUT(rowperm),
5060  std::vector<casadi_int>& SWIG_OUTPUT(colperm),
5061  std::vector<casadi_int>& SWIG_OUTPUT(rowblock),
5062  std::vector<casadi_int>& SWIG_OUTPUT(colblock),
5063  std::vector<casadi_int>& SWIG_OUTPUT(coarse_rowblock),
5064  std::vector<casadi_int>& SWIG_OUTPUT(coarse_colblock)) const;
5065 
5117  std::vector<casadi_int> amd() const;
5118 
5119 #ifndef SWIG
5123  void spsolve(bvec_t* X, bvec_t* B, bool tr) const;
5124 #endif // SWIG
5125 
5173  std::vector<casadi_int> find(bool ind1=SWIG_IND1) const;
5174 
5175 #ifndef SWIG
5177  void find(std::vector<casadi_int>& loc, bool ind1=false) const;
5178 #endif // SWIG
5179 
5274  Sparsity uni_coloring(const Sparsity& AT=Sparsity(),
5275  casadi_int cutoff = std::numeric_limits<casadi_int>::max()) const;
5276 
5360  Sparsity star_coloring_new(std::vector<casadi_int>& SWIG_OUTPUT(which_color),
5361  const Dict& opts=Dict()) const;
5362 
5426  Sparsity star_coloring(casadi_int ordering = 1,
5427  casadi_int cutoff = std::numeric_limits<casadi_int>::max()) const;
5428 
5473  Sparsity star_coloring2(casadi_int ordering = 1,
5474  casadi_int cutoff = std::numeric_limits<casadi_int>::max()) const;
5475 
5509  std::vector<casadi_int> largest_first() const;
5510 
5606  Sparsity pmult(const std::vector<casadi_int>& p,
5607  bool permute_rows=true, bool permute_columns=true,
5608  bool invert_permutation=false) const;
5609 
5658  std::string dim(bool with_nz=false) const;
5659 
5714  std::string postfix_dim() const;
5715 
5755  std::string repr_el(casadi_int k) const;
5756 
5766  void spy(std::ostream &stream=casadi::uout()) const;
5767 
5832  void spy_matlab(const std::string& mfile) const;
5833 
5921  void export_code(const std::string& lang, std::ostream &stream=casadi::uout(),
5922  const Dict& options=Dict()) const;
5923 
5925  static std::string type_name() {return "Sparsity";}
5926 
5964  std::size_t hash() const;
5965 
5967  static bool test_cast(const SharedObjectInternal* ptr);
5968 
6038  static Sparsity kkt(const Sparsity& H, const Sparsity& J,
6039  bool with_x_diag=true, bool with_lam_g_diag=true);
6040 
6041 #ifndef SWIG
6047  template<typename T>
6048  void set(T* data, const T* val_data, const Sparsity& val_sp) const;
6049 
6055  template<typename T>
6056  void add(T* data, const T* val_data, const Sparsity& val_sp) const;
6057 
6063  template<typename T>
6064  void bor(T* data, const T* val_data, const Sparsity& val_sp) const;
6065 
6066  static std::string file_format(const std::string& filename,
6067  const std::string& format_hint, const std::set<std::string>& file_formats);
6068  static std::set<std::string> file_formats;
6069  private:
6070 
6072  void assign_cached(casadi_int nrow, casadi_int ncol, const std::vector<casadi_int>& colind,
6073  const std::vector<casadi_int>& row, bool order_rows=false);
6074 
6076  void assign_cached(casadi_int nrow, casadi_int ncol,
6077  const casadi_int* colind, const casadi_int* row, bool order_rows=false);
6078 
6079 #endif //SWIG
6080  };
6081 
6085  CASADI_EXPORT std::size_t hash_sparsity(casadi_int nrow, casadi_int ncol,
6086  const std::vector<casadi_int>& colind,
6087  const std::vector<casadi_int>& row);
6088 
6089  CASADI_EXPORT std::size_t hash_sparsity(casadi_int nrow, casadi_int ncol,
6090  const casadi_int* colind,
6091  const casadi_int* row);
6092 
6093 #ifndef SWIG
6094  // Template instantiations
6095  template<typename DataType>
6096  void Sparsity::set(DataType* data, const DataType* val_data, const Sparsity& val_sp) const {
6097  // Get dimensions of this
6098  const casadi_int sz = nnz();
6099  const casadi_int sz1 = size1();
6100  const casadi_int sz2 = size2();
6101 
6102  // Get dimensions of assigning matrix
6103  const casadi_int val_sz = val_sp.nnz();
6104  const casadi_int val_sz1 = val_sp.size1();
6105  const casadi_int val_sz2 = val_sp.size2();
6106  const casadi_int val_nel = val_sz1*val_sz2;
6107 
6108  // Check if sparsity matches
6109  if (val_sp==*this) {
6110  std::copy(val_data, val_data+sz, data);
6111  } else if (this->is_empty()) {
6112  // Quick return
6113  return;
6114  } else if (val_sp.is_empty()) {
6115  // Quick return
6116  return;
6117  } else if (val_nel==1) { // if scalar
6118  std::fill(data, data+sz, val_sz==0 ? DataType(0) : val_data[0]);
6119  } else if (sz2==val_sz2 && sz1==val_sz1) {
6120  // Matching dimensions
6121  // Sparsity
6122  const casadi_int* c = row();
6123  const casadi_int* rind = colind();
6124  const casadi_int* v_c = val_sp.row();
6125  const casadi_int* v_rind = val_sp.colind();
6126 
6127  // For all columns
6128  for (casadi_int i=0; i<sz2; ++i) {
6129 
6130  // Nonzero of the assigning matrix
6131  casadi_int v_el = v_rind[i];
6132 
6133  // First nonzero of the following column
6134  casadi_int v_el_end = v_rind[i+1];
6135 
6136  // Next row of the assigning matrix
6137  casadi_int v_j = v_el<v_el_end ? v_c[v_el] : sz1;
6138 
6139  // Assign all nonzeros
6140  for (casadi_int el=rind[i]; el!=rind[i+1]; ++el) {
6141 
6142  // Get row
6143  casadi_int j=c[el];
6144 
6145  // Forward the assigning nonzero
6146  while (v_j<j) {
6147  v_el++;
6148  v_j = v_el<v_el_end ? v_c[v_el] : sz1;
6149  }
6150 
6151  // Assign nonzero
6152  if (v_j==j) {
6153  data[el] = val_data[v_el++];
6154  v_j = v_el<v_el_end ? v_c[v_el] : sz1;
6155  } else {
6156  data[el] = 0;
6157  }
6158  }
6159  }
6160  } else if (sz1==val_sz2 && sz2==val_sz1 && sz2 == 1) {
6161  // Assign transposed (this is column)
6162  const casadi_int* v_cind = val_sp.colind();
6163  const casadi_int* r = row();
6164  for (casadi_int el=0; el<sz; ++el) {
6165  casadi_int rr=r[el];
6166  data[el] = v_cind[rr]==v_cind[rr+1] ? 0 : val_data[v_cind[rr]];
6167  }
6168  } else if (sz1==val_sz2 && sz2==val_sz1 && sz1 == 1) {
6169  // Assign transposed (this is row)
6170  for (casadi_int el=0; el<sz; ++el) data[el] = 0;
6171  const casadi_int* cind = colind();
6172  const casadi_int* v_r = val_sp.row();
6173  for (casadi_int el=0; el<val_sz; ++el) {
6174  casadi_int rr=v_r[el];
6175  if (cind[rr]!=cind[rr+1]) {
6176  data[cind[rr]] = val_data[el];
6177  }
6178  }
6179  } else {
6180  // Make sure that dimension matches
6181  casadi_error("Sparsity::set<DataType>: shape mismatch. lhs is "
6182  + dim() + ", while rhs is " + val_sp.dim() + ".");
6183  }
6184  }
6185 
6186  template<typename DataType>
6187  void Sparsity::add(DataType* data, const DataType* val_data, const Sparsity& val_sp) const {
6188  // Get dimensions of this
6189  const casadi_int sz = nnz();
6190  const casadi_int sz1 = size1();
6191  const casadi_int sz2 = size2();
6192  const casadi_int nel = sz1*sz2;
6193 
6194  // Get dimensions of assigning matrix
6195  const casadi_int val_sz = val_sp.nnz();
6196  const casadi_int val_sz1 = val_sp.size1();
6197  const casadi_int val_sz2 = val_sp.size2();
6198  const casadi_int val_nel = val_sz1*val_sz2;
6199 
6200  // Check if sparsity matches
6201  if (val_sp==*this) {
6202  for (casadi_int k=0; k<sz; ++k) {
6203  data[k] += val_data[k];
6204  }
6205  } else if (this->is_empty()) {
6206  // Quick return
6207  return;
6208  } else if (val_sp.is_empty()) {
6209  // Quick return
6210  return;
6211  } else if (val_nel==1) { // if scalar
6212  if (val_sz!=0) {
6213  for (casadi_int k=0; k<sz; ++k) {
6214  data[k] += val_data[0];
6215  }
6216  }
6217  } else {
6218  // Quick return if empty
6219  if (nel==0 && val_nel==0) return;
6220 
6221  // Make sure that dimension matches
6222  casadi_assert(sz2==val_sz2 && sz1==val_sz1,
6223  "Sparsity::add<DataType>: shape mismatch. lhs is "
6224  + dim() + ", while rhs is " + val_sp.dim() + ".");
6225 
6226  // Sparsity
6227  const casadi_int* c = row();
6228  const casadi_int* rind = colind();
6229  const casadi_int* v_c = val_sp.row();
6230  const casadi_int* v_rind = val_sp.colind();
6231 
6232  // For all columns
6233  for (casadi_int i=0; i<sz2; ++i) {
6234 
6235  // Nonzero of the assigning matrix
6236  casadi_int v_el = v_rind[i];
6237 
6238  // First nonzero of the following column
6239  casadi_int v_el_end = v_rind[i+1];
6240 
6241  // Next row of the assigning matrix
6242  casadi_int v_j = v_el<v_el_end ? v_c[v_el] : sz1;
6243 
6244  // Assign all nonzeros
6245  for (casadi_int el=rind[i]; el!=rind[i+1]; ++el) {
6246 
6247  // Get row
6248  casadi_int j=c[el];
6249 
6250  // Forward the assigning nonzero
6251  while (v_j<j) {
6252  v_el++;
6253  v_j = v_el<v_el_end ? v_c[v_el] : sz1;
6254  }
6255 
6256  // Assign nonzero
6257  if (v_j==j) {
6258  data[el] += val_data[v_el++];
6259  v_j = v_el<v_el_end ? v_c[v_el] : sz1;
6260  }
6261  }
6262  }
6263  }
6264  }
6265 
6266  template<typename DataType>
6267  void Sparsity::bor(DataType* data, const DataType* val_data, const Sparsity& val_sp) const {
6268  // Get dimensions of this
6269  const casadi_int sz = nnz();
6270  const casadi_int sz1 = size1();
6271  const casadi_int sz2 = size2();
6272  const casadi_int nel = sz1*sz2;
6273 
6274  // Get dimensions of assigning matrix
6275  const casadi_int val_sz = val_sp.nnz();
6276  const casadi_int val_sz1 = val_sp.size1();
6277  const casadi_int val_sz2 = val_sp.size2();
6278  const casadi_int val_nel = val_sz1*val_sz2;
6279 
6280  // Check if sparsity matches
6281  if (val_sp==*this) {
6282  for (casadi_int k=0; k<sz; ++k) {
6283  data[k] |= val_data[k];
6284  }
6285  } else if (this->is_empty()) {
6286  // Quick return
6287  return;
6288  } else if (val_sp.is_empty()) {
6289  // Quick return
6290  return;
6291  } else if (val_nel==1) { // if scalar
6292  if (val_sz!=0) {
6293  for (casadi_int k=0; k<sz; ++k) {
6294  data[k] |= val_data[0];
6295  }
6296  }
6297  } else {
6298  // Quick return if empty
6299  if (nel==0 && val_nel==0) return;
6300 
6301  // Make sure that dimension matches
6302  casadi_assert(sz2==val_sz2 && sz1==val_sz1,
6303  "Sparsity::add<DataType>: shape mismatch. lhs is "
6304  + dim() + ", while rhs is " + val_sp.dim() + ".");
6305 
6306  // Sparsity
6307  const casadi_int* c = row();
6308  const casadi_int* rind = colind();
6309  const casadi_int* v_c = val_sp.row();
6310  const casadi_int* v_rind = val_sp.colind();
6311 
6312  // For all columns
6313  for (casadi_int i=0; i<sz2; ++i) {
6314 
6315  // Nonzero of the assigning matrix
6316  casadi_int v_el = v_rind[i];
6317 
6318  // First nonzero of the following column
6319  casadi_int v_el_end = v_rind[i+1];
6320 
6321  // Next row of the assigning matrix
6322  casadi_int v_j = v_el<v_el_end ? v_c[v_el] : sz1;
6323 
6324  // Assign all nonzeros
6325  for (casadi_int el=rind[i]; el!=rind[i+1]; ++el) {
6326 
6327  // Get row
6328  casadi_int j=c[el];
6329 
6330  // Forward the assigning nonzero
6331  while (v_j<j) {
6332  v_el++;
6333  v_j = v_el<v_el_end ? v_c[v_el] : sz1;
6334  }
6335 
6336  // Assign nonzero
6337  if (v_j==j) {
6338  data[el] |= val_data[v_el++];
6339  v_j = v_el<v_el_end ? v_c[v_el] : sz1;
6340  }
6341  }
6342  }
6343  }
6344  }
6345 
6346 #endif //SWIG
6347 
6350  typedef std::map<std::string, Sparsity> SpDict;
6352 
6353 } // namespace casadi
6354 
6355 #endif // CASADI_SPARSITY_HPP
casadi::DeserializingStream
Helper class for Serialization.
Definition: serializing_stream.hpp:73
casadi::SerializingStream
Helper class for Serialization.
Definition: serializing_stream.hpp:227
casadi::SharedObject
GenericShared implements a reference counting framework similar for efficient and.
Definition: shared_object.hpp:79
casadi::Sparsity
General sparsity class.
Definition: sparsity.hpp:200
casadi::Sparsity::columns
casadi_int columns() const
Number of columns, Octave-style alias of size2.
Definition: sparsity.hpp:1549
casadi::Sparsity::bor
void bor(T *data, const T *val_data, const Sparsity &val_sp) const
Bitwise or of the nonzero entries of one sparsity pattern and the nonzero.
casadi::Sparsity::type_name
static std::string type_name()
Readable name of the public class.
Definition: sparsity.hpp:5925
casadi::Sparsity::CachingMap
std::unordered_multimap< std::size_t, WeakRef > CachingMap
Definition: sparsity.hpp:4620
casadi::Sparsity::size1
casadi_int size1() const
Number of rows.
Definition: sparsity.cpp:124
casadi::Sparsity::file_formats
static std::set< std::string > file_formats
Definition: sparsity.hpp:6068
casadi::Sparsity::dense
static Sparsity dense(const std::pair< casadi_int, casadi_int > &rc)
Create a dense rectangular sparsity pattern.
Definition: sparsity.hpp:521
casadi::Sparsity::diag
static Sparsity diag(casadi_int nrow)
Create a diagonal sparsity pattern.
Definition: sparsity.hpp:723
casadi::Sparsity::diag
static Sparsity diag(const std::pair< casadi_int, casadi_int > &rc)
Create a diagonal sparsity pattern.
Definition: sparsity.hpp:725
casadi::Sparsity::add
void add(T *data, const T *val_data, const Sparsity &val_sp) const
Add the nonzero entries of one sparsity pattern to the nonzero entries.
casadi::Sparsity::dim
std::string dim(bool with_nz=false) const
Dimensions as a string "nrow x ncol", optionally with the number of nonzeros.
Definition: sparsity.cpp:595
casadi::Sparsity::set
void set(T *data, const T *val_data, const Sparsity &val_sp) const
Assign the nonzero entries of one sparsity pattern to the nonzero.
casadi::Sparsity::mac
static Sparsity mac(const Sparsity &x, const Sparsity &y, const Sparsity &z, const std::string &="reference")
Accessed by SparsityInterface.
Definition: sparsity.hpp:3435
casadi::Sparsity::rows
casadi_int rows() const
Number of rows, Octave-style alias of size1.
Definition: sparsity.hpp:1543
casadi::Sparsity::nnz
casadi_int nnz() const
Number of structural nonzeros.
Definition: sparsity.cpp:150
casadi::Sparsity::size2
casadi_int size2() const
Number of columns.
Definition: sparsity.cpp:128
casadi::Sparsity::row
const casadi_int * row() const
Get a reference to row-vector,.
Definition: sparsity.cpp:166
casadi::Sparsity::scalar
static Sparsity scalar(bool dense_scalar=true)
Create a 1-by-1 sparsity pattern.
Definition: sparsity.hpp:454
casadi::Sparsity::is_empty
bool is_empty(bool both=false) const
Check whether the pattern has a zero dimension.
Definition: sparsity.cpp:146
casadi::Sparsity::operator!=
bool operator!=(const Sparsity &y) const
Check whether two sparsity patterns differ (negation of is_equal)
Definition: sparsity.hpp:1443
casadi::Sparsity::colind
const casadi_int * colind() const
Get a reference to the colindex of all column element (see class description)
Definition: sparsity.cpp:170
casadi::Sparsity::operator==
bool operator==(const Sparsity &y) const
Check whether two sparsity patterns are identical.
Definition: sparsity.hpp:1439
casadi::Sparsity::B
SparsityInterface< Sparsity > B
Base class.
Definition: sparsity.hpp:399
casadi
The casadi namespace.
Definition: archiver.cpp:28
casadi::is_equal
bool is_equal(double x, double y, casadi_int depth=0)
Definition: calculus.hpp:287
casadi::invert_permutation
std::vector< casadi_int > invert_permutation(const std::vector< casadi_int > &a)
inverse a permutation vector
Definition: casadi_misc.cpp:189
casadi::bvec_t
unsigned long long bvec_t
Definition: casadi_common.hpp:93
casadi::combine
Dict combine(const Dict &first, const Dict &second, bool recurse)
Combine two dicts. First has priority.
Definition: generic_type.cpp:829
casadi::find
std::vector< casadi_int > find(const std::vector< T > &v)
find nonzeros
Definition: casadi_misc.hpp:647
casadi::Dict
GenericType::Dict Dict
C++ equivalent of Python's dict or MATLAB's struct.
Definition: generic_type.hpp:258
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:1025
casadi::SpDict
std::map< std::string, Sparsity > SpDict
Definition: sparsity.hpp:6350
casadi::kron_contract
MX kron_contract(const MX &m, const MX &x, bool inner)
Kronecker contraction, the adjoint of kron.
Definition: mx.hpp:5213
casadi::is_permutation
bool is_permutation(const std::vector< casadi_int > &order)
Does the list represent a permutation?
Definition: casadi_misc.cpp:182
casadi::uout
std::ostream & uout()
Definition: casadi_logger.cpp:61
casadi::filename
std::string filename(const std::string &path)
Definition: ghc.cpp:55
casadi::SparsityStruct
Compact representation of a sparsity pattern.
Definition: sparsity.hpp:57
casadi::SparsityStruct::colind
const casadi_int * colind
Definition: sparsity.hpp:60
casadi::SparsityStruct::row
const casadi_int * row
Definition: sparsity.hpp:61