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 
598 
608 #ifndef SWIG
609  Sparsity combine(const Sparsity& y, bool f0x_is_zero, bool function0_is_zero,
610  std::vector<unsigned char>& mapping) const;
611 #endif // SWIG
612  Sparsity combine(const Sparsity& y, bool f0x_is_zero, bool function0_is_zero) const;
614 
616 
619 #ifndef SWIG
620  Sparsity unite(const Sparsity& y, std::vector<unsigned char>& mapping) const;
621 #endif // SWIG
622  Sparsity unite(const Sparsity& y) const;
623  Sparsity operator+(const Sparsity& b) const;
625 
627 
635 #ifndef SWIG
636  Sparsity intersect(const Sparsity& y,
637  std::vector<unsigned char>& mapping) const;
638 #endif // SWIG
639  Sparsity intersect(const Sparsity& y) const;
640  Sparsity operator*(const Sparsity& b) const;
642 
644  bool is_subset(const Sparsity& rhs) const;
645 
659  Sparsity sparsity_cast_mod(const Sparsity& X, const Sparsity& Y) const;
660 
662  Sparsity pattern_inverse() const;
663 
664 #ifndef SWIG
671  static void mul_sparsityF(const bvec_t* x, const Sparsity& x_sp,
672  const bvec_t* y, const Sparsity& y_sp,
673  bvec_t* z, const Sparsity& z_sp,
674  bvec_t* w);
675 
682  static void mul_sparsityR(bvec_t* x, const Sparsity& x_sp,
683  bvec_t* y, const Sparsity& y_sp,
684  bvec_t* z, const Sparsity& z_sp,
685  bvec_t* w);
686 
689 
692  static Sparsity horzcat(const std::vector<Sparsity> & sp);
693  static Sparsity vertcat(const std::vector<Sparsity> & sp);
694  static Sparsity blockcat(const std::vector< std::vector< Sparsity > > &v);
695  static Sparsity diagcat(const std::vector< Sparsity > &v);
696  static std::vector<Sparsity>
697  horzsplit(const Sparsity& x, const std::vector<casadi_int>& offset);
698  static std::vector<Sparsity>
699  vertsplit(const Sparsity& x, const std::vector<casadi_int>& offset);
700  static std::vector<Sparsity>
701  diagsplit(const Sparsity& x,
702  const std::vector<casadi_int>& offset1,
703  const std::vector<casadi_int>& offset2);
704  static Sparsity mtimes(const Sparsity& x, const Sparsity& y);
705  static Sparsity mac(const Sparsity& x, const Sparsity& y, const Sparsity& z) { return z;}
706  static Sparsity reshape(const Sparsity& x, casadi_int nrow, casadi_int ncol);
707  static Sparsity reshape(const Sparsity& x, const Sparsity& sp);
708  static Sparsity sparsity_cast(const Sparsity& x, const Sparsity& sp);
709  static casadi_int sprank(const Sparsity& x);
710  static casadi_int norm_0_mul(const Sparsity& x, const Sparsity& A);
711  static Sparsity kron(const Sparsity& a, const Sparsity& b);
712  static Sparsity triu(const Sparsity& x, bool includeDiagonal=true);
713  static Sparsity tril(const Sparsity& x, bool includeDiagonal=true);
714 
715  static Sparsity sum2(const Sparsity &x);
716  static Sparsity sum1(const Sparsity &x);
717 #endif //SWIG
718 
733  void enlarge(casadi_int nrow, casadi_int ncol,
734  const std::vector<casadi_int>& rr,
735  const std::vector<casadi_int>& cc, bool ind1=false);
736 
740  void enlargeRows(casadi_int nrow, const std::vector<casadi_int>& rr, bool ind1=false);
741 
745  void enlargeColumns(casadi_int ncol, const std::vector<casadi_int>& cc, bool ind1=false);
746 
750  Sparsity makeDense(std::vector<casadi_int>& SWIG_OUTPUT(mapping)) const;
751 
755  std::vector<casadi_int> erase(const std::vector<casadi_int>& rr,
756  const std::vector<casadi_int>& cc, bool ind1=false);
757 
761  std::vector<casadi_int> erase(const std::vector<casadi_int>& rr, bool ind1=false);
762 
764  void append(const Sparsity& sp);
765 
767  void appendColumns(const Sparsity& sp);
768 
770  bool is_scalar(bool scalar_and_dense=false) const;
771 
773  bool is_dense() const;
774 
778  bool is_row() const;
779 
783  bool is_column() const;
784 
788  bool is_vector() const;
789 
793  bool is_diag() const;
794 
798  bool is_square() const;
799 
803  bool is_symmetric() const;
804 
808  bool is_triu(bool strictly = false) const;
809 
813  bool is_tril(bool strictly = false) const;
814 
818  bool is_singular() const;
819 
830  bool is_permutation() const;
831 
842  bool is_selection(bool allow_empty=false) const;
843 
849  bool is_orthonormal(bool allow_empty=false) const;
850 
856  bool is_orthonormal_rows(bool allow_empty=false) const;
857 
863  bool is_orthonormal_columns(bool allow_empty=false) const;
864 
870  bool rowsSequential(bool strictly=true) const;
871 
878  void removeDuplicates(std::vector<casadi_int>& SWIG_INOUT(mapping));
879 
880 #ifndef SWIG
881  typedef std::unordered_multimap<std::size_t, WeakRef> CachingMap;
882 
884  static CachingMap& getCache();
885 
886 #ifdef CASADI_WITH_THREADSAFE_SYMBOLICS
887  // Safe access to CachingMap
888  static std::mutex cachingmap_mtx;
889 #endif //CASADI_WITH_THREADSAFE_SYMBOLICS
890 
892  static const Sparsity& getScalar();
893 
895  static const Sparsity& getScalarSparse();
896 
898  static const Sparsity& getEmpty();
899 
900 #endif //SWIG
901 
914  std::vector<casadi_int> etree(bool ata=false) const;
915 
925  Sparsity ldl(std::vector<casadi_int>& SWIG_OUTPUT(p), bool amd=true) const;
926 
936  void qr_sparse(Sparsity& SWIG_OUTPUT(V), Sparsity& SWIG_OUTPUT(R),
937  std::vector<casadi_int>& SWIG_OUTPUT(prinv),
938  std::vector<casadi_int>& SWIG_OUTPUT(pc), bool amd=true) const;
939 
945  casadi_int dfs(casadi_int j, casadi_int top, std::vector<casadi_int>& SWIG_INOUT(xi),
946  std::vector<casadi_int>& SWIG_INOUT(pstack),
947  const std::vector<casadi_int>& pinv, std::vector<bool>& SWIG_INOUT(marked)) const;
948 
972  casadi_int scc(std::vector<casadi_int>& SWIG_OUTPUT(index),
973  std::vector<casadi_int>& SWIG_OUTPUT(offset)) const;
974 
996  casadi_int btf(std::vector<casadi_int>& SWIG_OUTPUT(rowperm),
997  std::vector<casadi_int>& SWIG_OUTPUT(colperm),
998  std::vector<casadi_int>& SWIG_OUTPUT(rowblock),
999  std::vector<casadi_int>& SWIG_OUTPUT(colblock),
1000  std::vector<casadi_int>& SWIG_OUTPUT(coarse_rowblock),
1001  std::vector<casadi_int>& SWIG_OUTPUT(coarse_colblock)) const;
1002 
1015  std::vector<casadi_int> amd() const;
1016 
1017 #ifndef SWIG
1021  void spsolve(bvec_t* X, bvec_t* B, bool tr) const;
1022 #endif // SWIG
1023 
1035  std::vector<casadi_int> find(bool ind1=SWIG_IND1) const;
1036 
1037 #ifndef SWIG
1039  void find(std::vector<casadi_int>& loc, bool ind1=false) const;
1040 #endif // SWIG
1041 
1047  Sparsity uni_coloring(const Sparsity& AT=Sparsity(),
1048  casadi_int cutoff = std::numeric_limits<casadi_int>::max()) const;
1049 
1061  Sparsity star_coloring(casadi_int ordering = 1,
1062  casadi_int cutoff = std::numeric_limits<casadi_int>::max()) const;
1063 
1075  Sparsity star_coloring2(casadi_int ordering = 1,
1076  casadi_int cutoff = std::numeric_limits<casadi_int>::max()) const;
1077 
1081  std::vector<casadi_int> largest_first() const;
1082 
1090  Sparsity pmult(const std::vector<casadi_int>& p,
1091  bool permute_rows=true, bool permute_columns=true,
1092  bool invert_permutation=false) const;
1093 
1095  std::string dim(bool with_nz=false) const;
1096 
1107  std::string postfix_dim() const;
1108 
1110  std::string repr_el(casadi_int k) const;
1111 
1115  void spy(std::ostream &stream=casadi::uout()) const;
1116 
1122  void spy_matlab(const std::string& mfile) const;
1123 
1136  void export_code(const std::string& lang, std::ostream &stream=casadi::uout(),
1137  const Dict& options=Dict()) const;
1138 
1140  static std::string type_name() {return "Sparsity";}
1141 
1142  // Hash the sparsity pattern
1143  std::size_t hash() const;
1144 
1146  static bool test_cast(const SharedObjectInternal* ptr);
1147 
1153  static Sparsity kkt(const Sparsity& H, const Sparsity& J,
1154  bool with_x_diag=true, bool with_lam_g_diag=true);
1155 
1156 #ifndef SWIG
1162  template<typename T>
1163  void set(T* data, const T* val_data, const Sparsity& val_sp) const;
1164 
1170  template<typename T>
1171  void add(T* data, const T* val_data, const Sparsity& val_sp) const;
1172 
1178  template<typename T>
1179  void bor(T* data, const T* val_data, const Sparsity& val_sp) const;
1180 
1181  static std::string file_format(const std::string& filename,
1182  const std::string& format_hint, const std::set<std::string>& file_formats);
1183  static std::set<std::string> file_formats;
1184  private:
1185 
1187  void assign_cached(casadi_int nrow, casadi_int ncol, const std::vector<casadi_int>& colind,
1188  const std::vector<casadi_int>& row, bool order_rows=false);
1189 
1191  void assign_cached(casadi_int nrow, casadi_int ncol,
1192  const casadi_int* colind, const casadi_int* row, bool order_rows=false);
1193 
1194 #endif //SWIG
1195  };
1196 
1200  template<typename T>
1201  inline size_t hash_value(T v) { return size_t(v);}
1202 
1206  template<typename T>
1207  inline void hash_combine(std::size_t& seed, T v) {
1208  seed ^= hash_value(v) + 0x9e3779b9 + (seed << 6) + (seed >> 2);
1209  }
1210 
1214  template<typename T>
1215  inline void hash_combine(std::size_t& seed, const T* v, std::size_t sz) {
1216  for (casadi_int i=0; i<sz; ++i) hash_combine(seed, v[i]);
1217  }
1218 
1222  template<typename T>
1223  inline void hash_combine(std::size_t& seed, const std::vector<T>& v) {
1224  hash_combine(seed, get_ptr(v), v.size());
1225  }
1226 
1227  template<>
1228  inline size_t hash_value(std::string v) {
1229  size_t seed = 0;
1230  hash_combine(seed, v.c_str(), v.size());
1231  return seed;
1232  }
1233 
1237  CASADI_EXPORT std::size_t hash_sparsity(casadi_int nrow, casadi_int ncol,
1238  const std::vector<casadi_int>& colind,
1239  const std::vector<casadi_int>& row);
1240 
1241  CASADI_EXPORT std::size_t hash_sparsity(casadi_int nrow, casadi_int ncol,
1242  const casadi_int* colind,
1243  const casadi_int* row);
1244 
1245 #ifndef SWIG
1246  // Template instantiations
1247  template<typename DataType>
1248  void Sparsity::set(DataType* data, const DataType* val_data, const Sparsity& val_sp) const {
1249  // Get dimensions of this
1250  const casadi_int sz = nnz();
1251  const casadi_int sz1 = size1();
1252  const casadi_int sz2 = size2();
1253 
1254  // Get dimensions of assigning matrix
1255  const casadi_int val_sz = val_sp.nnz();
1256  const casadi_int val_sz1 = val_sp.size1();
1257  const casadi_int val_sz2 = val_sp.size2();
1258  const casadi_int val_nel = val_sz1*val_sz2;
1259 
1260  // Check if sparsity matches
1261  if (val_sp==*this) {
1262  std::copy(val_data, val_data+sz, data);
1263  } else if (this->is_empty()) {
1264  // Quick return
1265  return;
1266  } else if (val_sp.is_empty()) {
1267  // Quick return
1268  return;
1269  } else if (val_nel==1) { // if scalar
1270  std::fill(data, data+sz, val_sz==0 ? DataType(0) : val_data[0]);
1271  } else if (sz2==val_sz2 && sz1==val_sz1) {
1272  // Matching dimensions
1273  // Sparsity
1274  const casadi_int* c = row();
1275  const casadi_int* rind = colind();
1276  const casadi_int* v_c = val_sp.row();
1277  const casadi_int* v_rind = val_sp.colind();
1278 
1279  // For all columns
1280  for (casadi_int i=0; i<sz2; ++i) {
1281 
1282  // Nonzero of the assigning matrix
1283  casadi_int v_el = v_rind[i];
1284 
1285  // First nonzero of the following column
1286  casadi_int v_el_end = v_rind[i+1];
1287 
1288  // Next row of the assigning matrix
1289  casadi_int v_j = v_el<v_el_end ? v_c[v_el] : sz1;
1290 
1291  // Assign all nonzeros
1292  for (casadi_int el=rind[i]; el!=rind[i+1]; ++el) {
1293 
1294  // Get row
1295  casadi_int j=c[el];
1296 
1297  // Forward the assigning nonzero
1298  while (v_j<j) {
1299  v_el++;
1300  v_j = v_el<v_el_end ? v_c[v_el] : sz1;
1301  }
1302 
1303  // Assign nonzero
1304  if (v_j==j) {
1305  data[el] = val_data[v_el++];
1306  v_j = v_el<v_el_end ? v_c[v_el] : sz1;
1307  } else {
1308  data[el] = 0;
1309  }
1310  }
1311  }
1312  } else if (sz1==val_sz2 && sz2==val_sz1 && sz2 == 1) {
1313  // Assign transposed (this is column)
1314  const casadi_int* v_cind = val_sp.colind();
1315  const casadi_int* r = row();
1316  for (casadi_int el=0; el<sz; ++el) {
1317  casadi_int rr=r[el];
1318  data[el] = v_cind[rr]==v_cind[rr+1] ? 0 : val_data[v_cind[rr]];
1319  }
1320  } else if (sz1==val_sz2 && sz2==val_sz1 && sz1 == 1) {
1321  // Assign transposed (this is row)
1322  for (casadi_int el=0; el<sz; ++el) data[el] = 0;
1323  const casadi_int* cind = colind();
1324  const casadi_int* v_r = val_sp.row();
1325  for (casadi_int el=0; el<val_sz; ++el) {
1326  casadi_int rr=v_r[el];
1327  if (cind[rr]!=cind[rr+1]) {
1328  data[cind[rr]] = val_data[el];
1329  }
1330  }
1331  } else {
1332  // Make sure that dimension matches
1333  casadi_error("Sparsity::set<DataType>: shape mismatch. lhs is "
1334  + dim() + ", while rhs is " + val_sp.dim() + ".");
1335  }
1336  }
1337 
1338  template<typename DataType>
1339  void Sparsity::add(DataType* data, const DataType* val_data, const Sparsity& val_sp) const {
1340  // Get dimensions of this
1341  const casadi_int sz = nnz();
1342  const casadi_int sz1 = size1();
1343  const casadi_int sz2 = size2();
1344  const casadi_int nel = sz1*sz2;
1345 
1346  // Get dimensions of assigning matrix
1347  const casadi_int val_sz = val_sp.nnz();
1348  const casadi_int val_sz1 = val_sp.size1();
1349  const casadi_int val_sz2 = val_sp.size2();
1350  const casadi_int val_nel = val_sz1*val_sz2;
1351 
1352  // Check if sparsity matches
1353  if (val_sp==*this) {
1354  for (casadi_int k=0; k<sz; ++k) {
1355  data[k] += val_data[k];
1356  }
1357  } else if (this->is_empty()) {
1358  // Quick return
1359  return;
1360  } else if (val_sp.is_empty()) {
1361  // Quick return
1362  return;
1363  } else if (val_nel==1) { // if scalar
1364  if (val_sz!=0) {
1365  for (casadi_int k=0; k<sz; ++k) {
1366  data[k] += val_data[0];
1367  }
1368  }
1369  } else {
1370  // Quick return if empty
1371  if (nel==0 && val_nel==0) return;
1372 
1373  // Make sure that dimension matches
1374  casadi_assert(sz2==val_sz2 && sz1==val_sz1,
1375  "Sparsity::add<DataType>: shape mismatch. lhs is "
1376  + dim() + ", while rhs is " + val_sp.dim() + ".");
1377 
1378  // Sparsity
1379  const casadi_int* c = row();
1380  const casadi_int* rind = colind();
1381  const casadi_int* v_c = val_sp.row();
1382  const casadi_int* v_rind = val_sp.colind();
1383 
1384  // For all columns
1385  for (casadi_int i=0; i<sz2; ++i) {
1386 
1387  // Nonzero of the assigning matrix
1388  casadi_int v_el = v_rind[i];
1389 
1390  // First nonzero of the following column
1391  casadi_int v_el_end = v_rind[i+1];
1392 
1393  // Next row of the assigning matrix
1394  casadi_int v_j = v_el<v_el_end ? v_c[v_el] : sz1;
1395 
1396  // Assign all nonzeros
1397  for (casadi_int el=rind[i]; el!=rind[i+1]; ++el) {
1398 
1399  // Get row
1400  casadi_int j=c[el];
1401 
1402  // Forward the assigning nonzero
1403  while (v_j<j) {
1404  v_el++;
1405  v_j = v_el<v_el_end ? v_c[v_el] : sz1;
1406  }
1407 
1408  // Assign nonzero
1409  if (v_j==j) {
1410  data[el] += val_data[v_el++];
1411  v_j = v_el<v_el_end ? v_c[v_el] : sz1;
1412  }
1413  }
1414  }
1415  }
1416  }
1417 
1418  template<typename DataType>
1419  void Sparsity::bor(DataType* data, const DataType* val_data, const Sparsity& val_sp) const {
1420  // Get dimensions of this
1421  const casadi_int sz = nnz();
1422  const casadi_int sz1 = size1();
1423  const casadi_int sz2 = size2();
1424  const casadi_int nel = sz1*sz2;
1425 
1426  // Get dimensions of assigning matrix
1427  const casadi_int val_sz = val_sp.nnz();
1428  const casadi_int val_sz1 = val_sp.size1();
1429  const casadi_int val_sz2 = val_sp.size2();
1430  const casadi_int val_nel = val_sz1*val_sz2;
1431 
1432  // Check if sparsity matches
1433  if (val_sp==*this) {
1434  for (casadi_int k=0; k<sz; ++k) {
1435  data[k] |= val_data[k];
1436  }
1437  } else if (this->is_empty()) {
1438  // Quick return
1439  return;
1440  } else if (val_sp.is_empty()) {
1441  // Quick return
1442  return;
1443  } else if (val_nel==1) { // if scalar
1444  if (val_sz!=0) {
1445  for (casadi_int k=0; k<sz; ++k) {
1446  data[k] |= val_data[0];
1447  }
1448  }
1449  } else {
1450  // Quick return if empty
1451  if (nel==0 && val_nel==0) return;
1452 
1453  // Make sure that dimension matches
1454  casadi_assert(sz2==val_sz2 && sz1==val_sz1,
1455  "Sparsity::add<DataType>: shape mismatch. lhs is "
1456  + dim() + ", while rhs is " + val_sp.dim() + ".");
1457 
1458  // Sparsity
1459  const casadi_int* c = row();
1460  const casadi_int* rind = colind();
1461  const casadi_int* v_c = val_sp.row();
1462  const casadi_int* v_rind = val_sp.colind();
1463 
1464  // For all columns
1465  for (casadi_int i=0; i<sz2; ++i) {
1466 
1467  // Nonzero of the assigning matrix
1468  casadi_int v_el = v_rind[i];
1469 
1470  // First nonzero of the following column
1471  casadi_int v_el_end = v_rind[i+1];
1472 
1473  // Next row of the assigning matrix
1474  casadi_int v_j = v_el<v_el_end ? v_c[v_el] : sz1;
1475 
1476  // Assign all nonzeros
1477  for (casadi_int el=rind[i]; el!=rind[i+1]; ++el) {
1478 
1479  // Get row
1480  casadi_int j=c[el];
1481 
1482  // Forward the assigning nonzero
1483  while (v_j<j) {
1484  v_el++;
1485  v_j = v_el<v_el_end ? v_c[v_el] : sz1;
1486  }
1487 
1488  // Assign nonzero
1489  if (v_j==j) {
1490  data[el] |= val_data[v_el++];
1491  v_j = v_el<v_el_end ? v_c[v_el] : sz1;
1492  }
1493  }
1494  }
1495  }
1496  }
1497 
1498 #endif //SWIG
1499 
1502  typedef std::map<std::string, Sparsity> SpDict;
1504 
1505 } // namespace casadi
1506 
1507 #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:215
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:1140
casadi::Sparsity::CachingMap
std::unordered_multimap< std::size_t, WeakRef > CachingMap
Enlarge matrix.
Definition: sparsity.hpp:881
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:1183
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:587
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::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::mac
static Sparsity mac(const Sparsity &x, const Sparsity &y, const Sparsity &z)
Enlarge matrix.
Definition: sparsity.hpp:705
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:281
casadi::invert_permutation
std::vector< casadi_int > invert_permutation(const std::vector< casadi_int > &a)
inverse a permutation vector
Definition: casadi_misc.cpp:178
casadi::bvec_t
unsigned long long bvec_t
Definition: casadi_common.hpp:76
casadi::combine
Dict combine(const Dict &first, const Dict &second, bool recurse)
Combine two dicts. First has priority.
Definition: generic_type.cpp:796
casadi::find
std::vector< casadi_int > find(const std::vector< T > &v)
find nonzeros
Definition: casadi_misc.hpp:618
casadi::Dict
GenericType::Dict Dict
C++ equivalent of Python's dict or MATLAB's struct.
Definition: generic_type.hpp:258
casadi::hash_combine
void hash_combine(std::size_t &seed, T v)
Generate a hash value incrementally (function taken from boost)
Definition: sparsity.hpp:1207
casadi::hash_sparsity
std::size_t hash_sparsity(casadi_int nrow, casadi_int ncol, const std::vector< casadi_int > &colind, const std::vector< casadi_int > &row)
Hash a sparsity pattern.
Definition: sparsity.cpp:996
casadi::hash_value
size_t hash_value(T v)
Hash value of an integer.
Definition: sparsity.hpp:1201
casadi::SpDict
std::map< std::string, Sparsity > SpDict
Definition: sparsity.hpp:1502
casadi::get_ptr
T * get_ptr(std::vector< T > &v)
Get a pointer to the data contained in the vector.
Definition: casadi_misc.hpp:857
casadi::is_permutation
bool is_permutation(const std::vector< casadi_int > &order)
Does the list represent a permutation?
Definition: casadi_misc.cpp:171
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