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 
103  class CASADI_EXPORT Sparsity
104  : public SharedObject,
105  public SWIG_IF_ELSE(SparsityInterfaceCommon, SparsityInterface<Sparsity>),
106  public SWIG_IF_ELSE(PrintableCommon, Printable<Sparsity>) {
107  public:
108 
110  explicit Sparsity(casadi_int dummy=0);
111 
115  Sparsity(casadi_int nrow, casadi_int ncol);
116 
118  Sparsity(casadi_int nrow, casadi_int ncol,
119  const std::vector<casadi_int>& colind, const std::vector<casadi_int>& row,
120  bool order_rows=false);
121 
125  explicit Sparsity(const std::pair<casadi_int, casadi_int>& rc);
126 
127 #ifndef SWIG
129  Sparsity(casadi_int nrow, casadi_int ncol, const casadi_int* colind, const casadi_int* row,
130  bool order_rows=false);
131 
135  static Sparsity create(SparsityInternal *node);
136 
138  typedef SparsityInterface<Sparsity> B;
139 
141  using B::horzsplit;
142  using B::diagsplit;
143  using B::vertsplit;
144  using B::mtimes;
145 
146  SparsityInternal* get() const;
147 #endif
148 
153  static Sparsity scalar(bool dense_scalar=true)
154  { return dense_scalar ? dense(1, 1) : Sparsity(1, 1); }
156 
161  static Sparsity dense(casadi_int nrow, casadi_int ncol=1);
162  static Sparsity dense(const std::pair<casadi_int, casadi_int> &rc) {
163  return dense(rc.first, rc.second);
164  }
166 
173  static Sparsity unit(casadi_int n, casadi_int el);
175 
179  static Sparsity upper(casadi_int n);
180 
184  static Sparsity lower(casadi_int n);
185 
190  static Sparsity diag(casadi_int nrow) { return diag(nrow, nrow);}
191  static Sparsity diag(casadi_int nrow, casadi_int ncol);
192  static Sparsity diag(const std::pair<casadi_int, casadi_int> &rc) {
193  return diag(rc.first, rc.second);
194  }
196 
204  static Sparsity band(casadi_int n, casadi_int p);
205 
212  static Sparsity banded(casadi_int n, casadi_int p);
213 
217  static Sparsity rowcol(const std::vector<casadi_int>& row,
218  const std::vector<casadi_int>& col,
219  casadi_int nrow, casadi_int ncol);
220 
222 
225  static Sparsity triplet(casadi_int nrow, casadi_int ncol,
226  const std::vector<casadi_int>& row, const std::vector<casadi_int>& col,
227  std::vector<casadi_int>& SWIG_OUTPUT(mapping), bool invert_mapping);
228  static Sparsity triplet(casadi_int nrow, casadi_int ncol, const std::vector<casadi_int>& row,
229  const std::vector<casadi_int>& col);
231 
237  static Sparsity nonzeros(casadi_int nrow, casadi_int ncol, const std::vector<casadi_int>& nz,
238  bool ind1=SWIG_IND1);
239 
248  static Sparsity compressed(const std::vector<casadi_int>& v, bool order_rows=false);
249 #ifndef SWIG
250  static Sparsity compressed(const casadi_int* v, bool order_rows=false);
251 #endif // SWIG
253 
264  static Sparsity permutation(const std::vector<casadi_int>& p, bool invert=false);
265 
269  const std::vector<casadi_int> permutation_vector(bool invert=false) const;
270 
276  Sparsity get_diag(std::vector<casadi_int>& SWIG_OUTPUT(mapping)) const;
277 
279  std::vector<casadi_int> compress(bool canonical=true) const;
280 
281 #ifndef SWIG
283  const SparsityInternal* operator->() const;
284 
286  const SparsityInternal& operator*() const;
287 #endif // SWIG
290  bool is_equal(const Sparsity& y) const;
291  bool is_equal(casadi_int nrow, casadi_int ncol, const std::vector<casadi_int>& colind,
292  const std::vector<casadi_int>& row) const;
293 #ifndef SWIG
294  bool is_equal(casadi_int nrow, casadi_int ncol,
295  const casadi_int* colind, const casadi_int* row) const;
296 #endif // SWIG
297 
298  bool operator==(const Sparsity& y) const { return is_equal(y);}
300 
302  bool operator!=(const Sparsity& y) const {return !is_equal(y);}
303 
305  bool is_stacked(const Sparsity& y, casadi_int n) const;
306 
307 #ifndef SWIG
314  operator const casadi_int*() const;
315 
319  operator const std::vector<casadi_int>&() const;
320 
324  operator SparsityStruct() const;
325 #endif // SWIG
326 
329 
331  casadi_int size1() const;
332 
334  casadi_int rows() const {return size1();}
335 
337  casadi_int size2() const;
338 
340  casadi_int columns() const {return size2();}
341 
348  casadi_int numel() const;
349 
355  double density() const;
356 
363  bool is_empty(bool both=false) const;
364 
370  casadi_int nnz() const;
371 
377  casadi_int nnz_upper(bool strictly=false) const;
378 
384  casadi_int nnz_lower(bool strictly=false) const;
385 
389  casadi_int nnz_diag() const;
390 
394  casadi_int bw_upper() const;
395 
399  casadi_int bw_lower() const;
400 
404  std::pair<casadi_int, casadi_int> size() const;
405 
409  casadi_int size(casadi_int axis) const;
411 
413  Dict info() const;
414 
420  void to_file(const std::string& filename, const std::string& format_hint="") const;
421 
422  static Sparsity from_file(const std::string& filename, const std::string& format_hint="");
423 
424 #ifndef SWIG
428  void serialize(std::ostream &stream) const;
429 #endif
430 
434  std::string serialize() const;
435 
439  static Sparsity deserialize(std::istream& stream);
440 
444  static Sparsity deserialize(const std::string& s);
445 
449  void serialize(SerializingStream& s) const;
450 
454  static Sparsity deserialize(DeserializingStream& s);
455 
456 #ifndef SWIG
462  const casadi_int* row() const;
463 
467  const casadi_int* colind() const;
468 #endif
469 
477  std::vector<casadi_int> get_row() const;
478 
485  std::vector<casadi_int> get_colind() const;
486 
490  casadi_int colind(casadi_int cc) const;
491 
495  casadi_int row(casadi_int el) const;
496 
503  std::vector<casadi_int> get_col() const;
504 
506  void resize(casadi_int nrow, casadi_int ncol);
507 
513  casadi_int add_nz(casadi_int rr, casadi_int cc);
514 
520  casadi_int get_nz(casadi_int rr, casadi_int cc) const;
521 
523  bool has_nz(casadi_int rr, casadi_int cc) const;
524 
530  std::vector<casadi_int> get_nz(const std::vector<casadi_int>& rr,
531  const std::vector<casadi_int>& cc) const;
532 
540  void get_nz(std::vector<casadi_int>& SWIG_INOUT(indices)) const;
541 
543  std::vector<casadi_int> get_lower() const;
544 
546  std::vector<casadi_int> get_upper() const;
547 
549  void get_ccs(std::vector<casadi_int>& SWIG_OUTPUT(colind),
550  std::vector<casadi_int>& SWIG_OUTPUT(row)) const;
551 
553  void get_crs(std::vector<casadi_int>& SWIG_OUTPUT(rowind),
554  std::vector<casadi_int>& SWIG_OUTPUT(col)) const;
555 
557  void get_triplet(std::vector<casadi_int>& SWIG_OUTPUT(row),
558  std::vector<casadi_int>& SWIG_OUTPUT(col)) const;
559 
566  Sparsity sub(const std::vector<casadi_int>& rr,
567  const std::vector<casadi_int>& cc,
568  std::vector<casadi_int>& SWIG_OUTPUT(mapping), bool ind1=false) const;
569 
576  Sparsity sub(const std::vector<casadi_int>& rr, const Sparsity& sp,
577  std::vector<casadi_int>& SWIG_OUTPUT(mapping), bool ind1=false) const;
578 
580  Sparsity T() const;
581 
588  Sparsity transpose(std::vector<casadi_int>& SWIG_OUTPUT(mapping),
589  bool invert_mapping=false) const;
590 
592  bool is_transpose(const Sparsity& y) const;
593 
595  bool is_reshape(const Sparsity& y) const;
596 
606  bool is_compactible(std::vector<casadi_int>& SWIG_OUTPUT(row),
607  std::vector<casadi_int>& SWIG_OUTPUT(col)) const;
608 
610 
620 #ifndef SWIG
621  Sparsity combine(const Sparsity& y, bool f0x_is_zero, bool function0_is_zero,
622  std::vector<unsigned char>& mapping) const;
623 #endif // SWIG
624  Sparsity combine(const Sparsity& y, bool f0x_is_zero, bool function0_is_zero) const;
626 
628 
631 #ifndef SWIG
632  Sparsity unite(const Sparsity& y, std::vector<unsigned char>& mapping) const;
633 #endif // SWIG
634  Sparsity unite(const Sparsity& y) const;
635  Sparsity operator+(const Sparsity& b) const;
637 
639 
647 #ifndef SWIG
648  Sparsity intersect(const Sparsity& y,
649  std::vector<unsigned char>& mapping) const;
650 #endif // SWIG
651  Sparsity intersect(const Sparsity& y) const;
652  Sparsity operator*(const Sparsity& b) const;
654 
656  bool is_subset(const Sparsity& rhs) const;
657 
671  Sparsity sparsity_cast_mod(const Sparsity& X, const Sparsity& Y) const;
672 
674  Sparsity pattern_inverse() const;
675 
676 #ifndef SWIG
683  static void mul_sparsityF(const bvec_t* x, const Sparsity& x_sp,
684  const bvec_t* y, const Sparsity& y_sp,
685  bvec_t* z, const Sparsity& z_sp,
686  bvec_t* w);
687 
695  static void mul_activityF(const bvec_t* x, const Sparsity& x_sp,
696  const bvec_t* y, const Sparsity& y_sp,
697  bvec_t* z, const Sparsity& z_sp,
698  bvec_t* w);
699 
706  static void mul_sparsityR(bvec_t* x, const Sparsity& x_sp,
707  bvec_t* y, const Sparsity& y_sp,
708  bvec_t* z, const Sparsity& z_sp,
709  bvec_t* w);
710 
713 
716  static Sparsity horzcat(const std::vector<Sparsity> & sp);
717  static Sparsity vertcat(const std::vector<Sparsity> & sp);
718  static Sparsity blockcat(const std::vector< std::vector< Sparsity > > &v);
719  static Sparsity diagcat(const std::vector< Sparsity > &v);
720  static std::vector<Sparsity>
721  horzsplit(const Sparsity& x, const std::vector<casadi_int>& offset);
722  static std::vector<Sparsity>
723  vertsplit(const Sparsity& x, const std::vector<casadi_int>& offset);
724  static std::vector<Sparsity>
725  diagsplit(const Sparsity& x,
726  const std::vector<casadi_int>& offset1,
727  const std::vector<casadi_int>& offset2);
728  static Sparsity mtimes(const Sparsity& x, const Sparsity& y,
729  const std::string& blas = "reference");
730  static Sparsity mac(const Sparsity& x, const Sparsity& y, const Sparsity& z,
731  const std::string& /*blas*/ = "reference") { return z;}
732  static Sparsity reshape(const Sparsity& x, casadi_int nrow, casadi_int ncol);
733  static Sparsity reshape(const Sparsity& x, const Sparsity& sp);
734  static Sparsity sparsity_cast(const Sparsity& x, const Sparsity& sp);
735  static casadi_int sprank(const Sparsity& x);
736  static casadi_int norm_0_mul(const Sparsity& x, const Sparsity& A);
737  static Sparsity kron(const Sparsity& a, const Sparsity& b);
738 
754  static Sparsity kron_contract(const Sparsity& sp_m, const Sparsity& sp_x, bool inner);
755  static Sparsity triu(const Sparsity& x, bool includeDiagonal=true);
756  static Sparsity tril(const Sparsity& x, bool includeDiagonal=true);
757 
758  static Sparsity sum2(const Sparsity &x);
759  static Sparsity sum1(const Sparsity &x);
760 #endif //SWIG
761 
776  void enlarge(casadi_int nrow, casadi_int ncol,
777  const std::vector<casadi_int>& rr,
778  const std::vector<casadi_int>& cc, bool ind1=false);
779 
783  void enlargeRows(casadi_int nrow, const std::vector<casadi_int>& rr, bool ind1=false);
784 
788  void enlargeColumns(casadi_int ncol, const std::vector<casadi_int>& cc, bool ind1=false);
789 
793  Sparsity makeDense(std::vector<casadi_int>& SWIG_OUTPUT(mapping)) const;
794 
798  std::vector<casadi_int> erase(const std::vector<casadi_int>& rr,
799  const std::vector<casadi_int>& cc, bool ind1=false);
800 
804  std::vector<casadi_int> erase(const std::vector<casadi_int>& rr, bool ind1=false);
805 
807  void append(const Sparsity& sp);
808 
810  void appendColumns(const Sparsity& sp);
811 
813  bool is_scalar(bool scalar_and_dense=false) const;
814 
816  bool is_dense() const;
817 
821  bool is_row() const;
822 
826  bool is_column() const;
827 
831  bool is_vector() const;
832 
836  bool is_diag() const;
837 
841  bool is_square() const;
842 
846  bool is_symmetric() const;
847 
851  bool is_triu(bool strictly = false) const;
852 
856  bool is_tril(bool strictly = false) const;
857 
861  bool is_singular() const;
862 
873  bool is_permutation() const;
874 
885  bool is_selection(bool allow_empty=false) const;
886 
892  bool is_orthonormal(bool allow_empty=false) const;
893 
899  bool is_orthonormal_rows(bool allow_empty=false) const;
900 
906  bool is_orthonormal_columns(bool allow_empty=false) const;
907 
913  bool rowsSequential(bool strictly=true) const;
914 
921  void removeDuplicates(std::vector<casadi_int>& SWIG_INOUT(mapping));
922 
923 #ifndef SWIG
924  typedef std::unordered_multimap<std::size_t, WeakRef> CachingMap;
925 
927  static CachingMap& getCache();
928 
929 #ifdef CASADI_WITH_THREADSAFE_SYMBOLICS
930  // Safe access to CachingMap
931  static std::mutex cachingmap_mtx;
932 #endif //CASADI_WITH_THREADSAFE_SYMBOLICS
933 
935  static const Sparsity& getScalar();
936 
938  static const Sparsity& getScalarSparse();
939 
941  static const Sparsity& getEmpty();
942 
943 #endif //SWIG
944 
957  std::vector<casadi_int> etree(bool ata=false) const;
958 
968  Sparsity ldl(std::vector<casadi_int>& SWIG_OUTPUT(p), bool amd=true) const;
969 
979  void qr_sparse(Sparsity& SWIG_OUTPUT(V), Sparsity& SWIG_OUTPUT(R),
980  std::vector<casadi_int>& SWIG_OUTPUT(prinv),
981  std::vector<casadi_int>& SWIG_OUTPUT(pc), bool amd=true) const;
982 
988  casadi_int dfs(casadi_int j, casadi_int top, std::vector<casadi_int>& SWIG_INOUT(xi),
989  std::vector<casadi_int>& SWIG_INOUT(pstack),
990  const std::vector<casadi_int>& pinv, std::vector<bool>& SWIG_INOUT(marked)) const;
991 
1015  casadi_int scc(std::vector<casadi_int>& SWIG_OUTPUT(index),
1016  std::vector<casadi_int>& SWIG_OUTPUT(offset)) const;
1017 
1039  casadi_int btf(std::vector<casadi_int>& SWIG_OUTPUT(rowperm),
1040  std::vector<casadi_int>& SWIG_OUTPUT(colperm),
1041  std::vector<casadi_int>& SWIG_OUTPUT(rowblock),
1042  std::vector<casadi_int>& SWIG_OUTPUT(colblock),
1043  std::vector<casadi_int>& SWIG_OUTPUT(coarse_rowblock),
1044  std::vector<casadi_int>& SWIG_OUTPUT(coarse_colblock)) const;
1045 
1058  std::vector<casadi_int> amd() const;
1059 
1060 #ifndef SWIG
1064  void spsolve(bvec_t* X, bvec_t* B, bool tr) const;
1065 #endif // SWIG
1066 
1078  std::vector<casadi_int> find(bool ind1=SWIG_IND1) const;
1079 
1080 #ifndef SWIG
1082  void find(std::vector<casadi_int>& loc, bool ind1=false) const;
1083 #endif // SWIG
1084 
1090  Sparsity uni_coloring(const Sparsity& AT=Sparsity(),
1091  casadi_int cutoff = std::numeric_limits<casadi_int>::max()) const;
1092 
1106  Sparsity star_coloring_new(std::vector<casadi_int>& SWIG_OUTPUT(which_color),
1107  const Dict& opts=Dict()) const;
1108 
1120  Sparsity star_coloring(casadi_int ordering = 1,
1121  casadi_int cutoff = std::numeric_limits<casadi_int>::max()) const;
1122 
1134  Sparsity star_coloring2(casadi_int ordering = 1,
1135  casadi_int cutoff = std::numeric_limits<casadi_int>::max()) const;
1136 
1140  std::vector<casadi_int> largest_first() const;
1141 
1149  Sparsity pmult(const std::vector<casadi_int>& p,
1150  bool permute_rows=true, bool permute_columns=true,
1151  bool invert_permutation=false) const;
1152 
1154  std::string dim(bool with_nz=false) const;
1155 
1166  std::string postfix_dim() const;
1167 
1169  std::string repr_el(casadi_int k) const;
1170 
1174  void spy(std::ostream &stream=casadi::uout()) const;
1175 
1181  void spy_matlab(const std::string& mfile) const;
1182 
1195  void export_code(const std::string& lang, std::ostream &stream=casadi::uout(),
1196  const Dict& options=Dict()) const;
1197 
1199  static std::string type_name() {return "Sparsity";}
1200 
1201  // Hash the sparsity pattern
1202  std::size_t hash() const;
1203 
1205  static bool test_cast(const SharedObjectInternal* ptr);
1206 
1212  static Sparsity kkt(const Sparsity& H, const Sparsity& J,
1213  bool with_x_diag=true, bool with_lam_g_diag=true);
1214 
1215 #ifndef SWIG
1221  template<typename T>
1222  void set(T* data, const T* val_data, const Sparsity& val_sp) const;
1223 
1229  template<typename T>
1230  void add(T* data, const T* val_data, const Sparsity& val_sp) const;
1231 
1237  template<typename T>
1238  void bor(T* data, const T* val_data, const Sparsity& val_sp) const;
1239 
1240  static std::string file_format(const std::string& filename,
1241  const std::string& format_hint, const std::set<std::string>& file_formats);
1242  static std::set<std::string> file_formats;
1243  private:
1244 
1246  void assign_cached(casadi_int nrow, casadi_int ncol, const std::vector<casadi_int>& colind,
1247  const std::vector<casadi_int>& row, bool order_rows=false);
1248 
1250  void assign_cached(casadi_int nrow, casadi_int ncol,
1251  const casadi_int* colind, const casadi_int* row, bool order_rows=false);
1252 
1253 #endif //SWIG
1254  };
1255 
1259  CASADI_EXPORT std::size_t hash_sparsity(casadi_int nrow, casadi_int ncol,
1260  const std::vector<casadi_int>& colind,
1261  const std::vector<casadi_int>& row);
1262 
1263  CASADI_EXPORT std::size_t hash_sparsity(casadi_int nrow, casadi_int ncol,
1264  const casadi_int* colind,
1265  const casadi_int* row);
1266 
1267 #ifndef SWIG
1268  // Template instantiations
1269  template<typename DataType>
1270  void Sparsity::set(DataType* data, const DataType* val_data, const Sparsity& val_sp) const {
1271  // Get dimensions of this
1272  const casadi_int sz = nnz();
1273  const casadi_int sz1 = size1();
1274  const casadi_int sz2 = size2();
1275 
1276  // Get dimensions of assigning matrix
1277  const casadi_int val_sz = val_sp.nnz();
1278  const casadi_int val_sz1 = val_sp.size1();
1279  const casadi_int val_sz2 = val_sp.size2();
1280  const casadi_int val_nel = val_sz1*val_sz2;
1281 
1282  // Check if sparsity matches
1283  if (val_sp==*this) {
1284  std::copy(val_data, val_data+sz, data);
1285  } else if (this->is_empty()) {
1286  // Quick return
1287  return;
1288  } else if (val_sp.is_empty()) {
1289  // Quick return
1290  return;
1291  } else if (val_nel==1) { // if scalar
1292  std::fill(data, data+sz, val_sz==0 ? DataType(0) : val_data[0]);
1293  } else if (sz2==val_sz2 && sz1==val_sz1) {
1294  // Matching dimensions
1295  // Sparsity
1296  const casadi_int* c = row();
1297  const casadi_int* rind = colind();
1298  const casadi_int* v_c = val_sp.row();
1299  const casadi_int* v_rind = val_sp.colind();
1300 
1301  // For all columns
1302  for (casadi_int i=0; i<sz2; ++i) {
1303 
1304  // Nonzero of the assigning matrix
1305  casadi_int v_el = v_rind[i];
1306 
1307  // First nonzero of the following column
1308  casadi_int v_el_end = v_rind[i+1];
1309 
1310  // Next row of the assigning matrix
1311  casadi_int v_j = v_el<v_el_end ? v_c[v_el] : sz1;
1312 
1313  // Assign all nonzeros
1314  for (casadi_int el=rind[i]; el!=rind[i+1]; ++el) {
1315 
1316  // Get row
1317  casadi_int j=c[el];
1318 
1319  // Forward the assigning nonzero
1320  while (v_j<j) {
1321  v_el++;
1322  v_j = v_el<v_el_end ? v_c[v_el] : sz1;
1323  }
1324 
1325  // Assign nonzero
1326  if (v_j==j) {
1327  data[el] = val_data[v_el++];
1328  v_j = v_el<v_el_end ? v_c[v_el] : sz1;
1329  } else {
1330  data[el] = 0;
1331  }
1332  }
1333  }
1334  } else if (sz1==val_sz2 && sz2==val_sz1 && sz2 == 1) {
1335  // Assign transposed (this is column)
1336  const casadi_int* v_cind = val_sp.colind();
1337  const casadi_int* r = row();
1338  for (casadi_int el=0; el<sz; ++el) {
1339  casadi_int rr=r[el];
1340  data[el] = v_cind[rr]==v_cind[rr+1] ? 0 : val_data[v_cind[rr]];
1341  }
1342  } else if (sz1==val_sz2 && sz2==val_sz1 && sz1 == 1) {
1343  // Assign transposed (this is row)
1344  for (casadi_int el=0; el<sz; ++el) data[el] = 0;
1345  const casadi_int* cind = colind();
1346  const casadi_int* v_r = val_sp.row();
1347  for (casadi_int el=0; el<val_sz; ++el) {
1348  casadi_int rr=v_r[el];
1349  if (cind[rr]!=cind[rr+1]) {
1350  data[cind[rr]] = val_data[el];
1351  }
1352  }
1353  } else {
1354  // Make sure that dimension matches
1355  casadi_error("Sparsity::set<DataType>: shape mismatch. lhs is "
1356  + dim() + ", while rhs is " + val_sp.dim() + ".");
1357  }
1358  }
1359 
1360  template<typename DataType>
1361  void Sparsity::add(DataType* data, const DataType* val_data, const Sparsity& val_sp) const {
1362  // Get dimensions of this
1363  const casadi_int sz = nnz();
1364  const casadi_int sz1 = size1();
1365  const casadi_int sz2 = size2();
1366  const casadi_int nel = sz1*sz2;
1367 
1368  // Get dimensions of assigning matrix
1369  const casadi_int val_sz = val_sp.nnz();
1370  const casadi_int val_sz1 = val_sp.size1();
1371  const casadi_int val_sz2 = val_sp.size2();
1372  const casadi_int val_nel = val_sz1*val_sz2;
1373 
1374  // Check if sparsity matches
1375  if (val_sp==*this) {
1376  for (casadi_int k=0; k<sz; ++k) {
1377  data[k] += val_data[k];
1378  }
1379  } else if (this->is_empty()) {
1380  // Quick return
1381  return;
1382  } else if (val_sp.is_empty()) {
1383  // Quick return
1384  return;
1385  } else if (val_nel==1) { // if scalar
1386  if (val_sz!=0) {
1387  for (casadi_int k=0; k<sz; ++k) {
1388  data[k] += val_data[0];
1389  }
1390  }
1391  } else {
1392  // Quick return if empty
1393  if (nel==0 && val_nel==0) return;
1394 
1395  // Make sure that dimension matches
1396  casadi_assert(sz2==val_sz2 && sz1==val_sz1,
1397  "Sparsity::add<DataType>: shape mismatch. lhs is "
1398  + dim() + ", while rhs is " + val_sp.dim() + ".");
1399 
1400  // Sparsity
1401  const casadi_int* c = row();
1402  const casadi_int* rind = colind();
1403  const casadi_int* v_c = val_sp.row();
1404  const casadi_int* v_rind = val_sp.colind();
1405 
1406  // For all columns
1407  for (casadi_int i=0; i<sz2; ++i) {
1408 
1409  // Nonzero of the assigning matrix
1410  casadi_int v_el = v_rind[i];
1411 
1412  // First nonzero of the following column
1413  casadi_int v_el_end = v_rind[i+1];
1414 
1415  // Next row of the assigning matrix
1416  casadi_int v_j = v_el<v_el_end ? v_c[v_el] : sz1;
1417 
1418  // Assign all nonzeros
1419  for (casadi_int el=rind[i]; el!=rind[i+1]; ++el) {
1420 
1421  // Get row
1422  casadi_int j=c[el];
1423 
1424  // Forward the assigning nonzero
1425  while (v_j<j) {
1426  v_el++;
1427  v_j = v_el<v_el_end ? v_c[v_el] : sz1;
1428  }
1429 
1430  // Assign nonzero
1431  if (v_j==j) {
1432  data[el] += val_data[v_el++];
1433  v_j = v_el<v_el_end ? v_c[v_el] : sz1;
1434  }
1435  }
1436  }
1437  }
1438  }
1439 
1440  template<typename DataType>
1441  void Sparsity::bor(DataType* data, const DataType* val_data, const Sparsity& val_sp) const {
1442  // Get dimensions of this
1443  const casadi_int sz = nnz();
1444  const casadi_int sz1 = size1();
1445  const casadi_int sz2 = size2();
1446  const casadi_int nel = sz1*sz2;
1447 
1448  // Get dimensions of assigning matrix
1449  const casadi_int val_sz = val_sp.nnz();
1450  const casadi_int val_sz1 = val_sp.size1();
1451  const casadi_int val_sz2 = val_sp.size2();
1452  const casadi_int val_nel = val_sz1*val_sz2;
1453 
1454  // Check if sparsity matches
1455  if (val_sp==*this) {
1456  for (casadi_int k=0; k<sz; ++k) {
1457  data[k] |= val_data[k];
1458  }
1459  } else if (this->is_empty()) {
1460  // Quick return
1461  return;
1462  } else if (val_sp.is_empty()) {
1463  // Quick return
1464  return;
1465  } else if (val_nel==1) { // if scalar
1466  if (val_sz!=0) {
1467  for (casadi_int k=0; k<sz; ++k) {
1468  data[k] |= val_data[0];
1469  }
1470  }
1471  } else {
1472  // Quick return if empty
1473  if (nel==0 && val_nel==0) return;
1474 
1475  // Make sure that dimension matches
1476  casadi_assert(sz2==val_sz2 && sz1==val_sz1,
1477  "Sparsity::add<DataType>: shape mismatch. lhs is "
1478  + dim() + ", while rhs is " + val_sp.dim() + ".");
1479 
1480  // Sparsity
1481  const casadi_int* c = row();
1482  const casadi_int* rind = colind();
1483  const casadi_int* v_c = val_sp.row();
1484  const casadi_int* v_rind = val_sp.colind();
1485 
1486  // For all columns
1487  for (casadi_int i=0; i<sz2; ++i) {
1488 
1489  // Nonzero of the assigning matrix
1490  casadi_int v_el = v_rind[i];
1491 
1492  // First nonzero of the following column
1493  casadi_int v_el_end = v_rind[i+1];
1494 
1495  // Next row of the assigning matrix
1496  casadi_int v_j = v_el<v_el_end ? v_c[v_el] : sz1;
1497 
1498  // Assign all nonzeros
1499  for (casadi_int el=rind[i]; el!=rind[i+1]; ++el) {
1500 
1501  // Get row
1502  casadi_int j=c[el];
1503 
1504  // Forward the assigning nonzero
1505  while (v_j<j) {
1506  v_el++;
1507  v_j = v_el<v_el_end ? v_c[v_el] : sz1;
1508  }
1509 
1510  // Assign nonzero
1511  if (v_j==j) {
1512  data[el] |= val_data[v_el++];
1513  v_j = v_el<v_el_end ? v_c[v_el] : sz1;
1514  }
1515  }
1516  }
1517  }
1518  }
1519 
1520 #endif //SWIG
1521 
1524  typedef std::map<std::string, Sparsity> SpDict;
1526 
1527 } // namespace casadi
1528 
1529 #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:106
casadi::Sparsity::columns
casadi_int columns() const
Get the number of columns, Octave-style syntax.
Definition: sparsity.hpp:340
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:1199
casadi::Sparsity::CachingMap
std::unordered_multimap< std::size_t, WeakRef > CachingMap
Enlarge matrix.
Definition: sparsity.hpp:924
casadi::Sparsity::size1
casadi_int size1() const
Get the number of rows.
Definition: sparsity.cpp:124
casadi::Sparsity::file_formats
static std::set< std::string > file_formats
Enlarge matrix.
Definition: sparsity.hpp:1242
casadi::Sparsity::dense
static Sparsity dense(const std::pair< casadi_int, casadi_int > &rc)
Create a dense rectangular sparsity pattern *.
Definition: sparsity.hpp:162
casadi::Sparsity::diag
static Sparsity diag(casadi_int nrow)
Create diagonal sparsity pattern *.
Definition: sparsity.hpp:190
casadi::Sparsity::diag
static Sparsity diag(const std::pair< casadi_int, casadi_int > &rc)
Create diagonal sparsity pattern *.
Definition: sparsity.hpp:192
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
Get the dimension as a string.
Definition: sparsity.cpp:588
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")
Enlarge matrix.
Definition: sparsity.hpp:730
casadi::Sparsity::rows
casadi_int rows() const
Get the number of rows, Octave-style syntax.
Definition: sparsity.hpp:334
casadi::Sparsity::nnz
casadi_int nnz() const
Get the number of (structural) non-zeros.
Definition: sparsity.cpp:148
casadi::Sparsity::size2
casadi_int size2() const
Get the number of columns.
Definition: sparsity.cpp:128
casadi::Sparsity::row
const casadi_int * row() const
Get a reference to row-vector,.
Definition: sparsity.cpp:164
casadi::Sparsity::scalar
static Sparsity scalar(bool dense_scalar=true)
Create a scalar sparsity pattern *.
Definition: sparsity.hpp:153
casadi::Sparsity::is_empty
bool is_empty(bool both=false) const
Check if the sparsity is empty.
Definition: sparsity.cpp:144
casadi::Sparsity::operator!=
bool operator!=(const Sparsity &y) const
Check if two sparsity patterns are difference.
Definition: sparsity.hpp:302
casadi::Sparsity::colind
const casadi_int * colind() const
Get a reference to the colindex of all column element (see class description)
Definition: sparsity.cpp:168
casadi::Sparsity::operator==
bool operator==(const Sparsity &y) const
Definition: sparsity.hpp:298
casadi::Sparsity::B
SparsityInterface< Sparsity > B
Base class.
Definition: sparsity.hpp:138
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:182
casadi::bvec_t
unsigned long long bvec_t
Definition: casadi_common.hpp:77
casadi::combine
Dict combine(const Dict &first, const Dict &second, bool recurse)
Combine two dicts. First has priority.
Definition: generic_type.cpp:827
casadi::find
std::vector< casadi_int > find(const std::vector< T > &v)
find nonzeros
Definition: casadi_misc.hpp:643
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:1012
casadi::SpDict
std::map< std::string, Sparsity > SpDict
Definition: sparsity.hpp:1524
casadi::kron_contract
MX kron_contract(const MX &m, const MX &x, bool inner)
Kronecker contraction.
Definition: mx.hpp:1121
casadi::is_permutation
bool is_permutation(const std::vector< casadi_int > &order)
Does the list represent a permutation?
Definition: casadi_misc.cpp:175
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