sparsity_internal.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_INTERNAL_HPP
28 #define CASADI_SPARSITY_INTERNAL_HPP
29 
30 #include "sparsity.hpp"
31 #include "shared_object.hpp"
32 #ifdef CASADI_WITH_THREAD
33 #ifdef CASADI_WITH_THREAD_MINGW
34 #include <mingw.mutex.h>
35 #else // CASADI_WITH_THREAD_MINGW
36 #include <mutex>
37 #endif // CASADI_WITH_THREAD_MINGW
38 #endif //CASADI_WITH_THREAD
39 
40 #include "filesystem_impl.hpp"
41 
43 
44 namespace casadi {
45 
46  class CASADI_EXPORT SparsityInternal : public SharedObjectInternal {
47  private:
57  std::vector<casadi_int> sp_;
58 
62  struct Btf {
63  casadi_int nb;
64  std::vector<casadi_int> rowperm, colperm;
65  std::vector<casadi_int> rowblock, colblock;
66  std::vector<casadi_int> coarse_rowblock, coarse_colblock;
67  };
68 
74  mutable Btf* btf_;
75 
76 #ifdef CASADI_WITH_THREADSAFE_SYMBOLICS
78  mutable std::mutex btf_mtx_;
79 #endif // CASADI_WITH_THREADSAFE_SYMBOLICS
80 
81  public:
83  SparsityInternal(casadi_int nrow, casadi_int ncol,
84  const casadi_int* colind, const casadi_int* row);
85 
87  ~SparsityInternal() override;
88 
92  inline const std::vector<casadi_int>& sp() const { return sp_;}
93 
97  inline casadi_int size1() const { return sp_[0];}
98 
102  inline casadi_int size2() const { return sp_[1];}
103 
107  inline const casadi_int* colind() const { return &sp_.front()+2;}
108 
112  inline const casadi_int* row() const { return colind()+size2()+1;}
113 
115  inline casadi_int nnz() const { return colind()[size2()];}
116 
122  Sparsity get_diag(std::vector<casadi_int>& mapping) const;
123 
125  bool has_diag() const;
126 
128  Sparsity drop_diag() const;
129 
137  casadi_int dfs(casadi_int j, casadi_int top, std::vector<casadi_int>& xi,
138  std::vector<casadi_int>& pstack,
139  const std::vector<casadi_int>& pinv, std::vector<bool>& marked) const;
140 
148  casadi_int scc(std::vector<casadi_int>& p, std::vector<casadi_int>& r) const;
149 
157  std::vector<casadi_int> amd() const;
158 
169  static void etree(const casadi_int* sp, casadi_int* parent, casadi_int *w, casadi_int ata);
170 
179  static casadi_int postorder_dfs(casadi_int j, casadi_int k, casadi_int* head,
180  const casadi_int* next, casadi_int* post, casadi_int* stack);
181 
192  static void postorder(const casadi_int* parent, casadi_int n, casadi_int* post, casadi_int* w);
193 
202  static casadi_int leaf(casadi_int i, casadi_int j, const casadi_int* first,
203  casadi_int* maxfirst,
204  casadi_int* prevleaf, casadi_int* ancestor, casadi_int* jleaf);
205 
216  static casadi_int qr_counts(const casadi_int* tr_sp, const casadi_int* parent,
217  const casadi_int* post, casadi_int* counts, casadi_int* w);
218 
230  static casadi_int qr_nnz(const casadi_int* sp, casadi_int* pinv, casadi_int* leftmost,
231  const casadi_int* parent, casadi_int* nrow_ext, casadi_int* w);
232 
241  static void qr_init(const casadi_int* sp, const casadi_int* sp_tr,
242  casadi_int* leftmost, casadi_int* parent, casadi_int* pinv,
243  casadi_int* nrow_ext, casadi_int* v_nnz, casadi_int* r_nnz, casadi_int* w);
244 
261  static void qr_sparsities(
262  const casadi_int* sp_a, casadi_int nrow_ext, casadi_int* sp_v, casadi_int* sp_r,
263  const casadi_int* leftmost, const casadi_int* parent, const casadi_int* pinv,
264  casadi_int* iw);
265 
278  static void ldl_colind(const casadi_int* sp, casadi_int* parent,
279  casadi_int* l_colind, casadi_int* w);
280 
291  static void ldl_row(const casadi_int* sp, const casadi_int* parent,
292  casadi_int* l_colind, casadi_int* l_row, casadi_int *w);
293 
295  Sparsity T() const;
296 
302  Sparsity transpose(std::vector<casadi_int>& mapping, bool invert_mapping=false) const;
303 
305  bool is_transpose(const SparsityInternal& y) const;
306 
308  bool is_reshape(const SparsityInternal& y) const;
309 
311  bool is_compactible(std::vector<casadi_int>& row, std::vector<casadi_int>& col) const;
312 
320  void bfs(casadi_int n, std::vector<casadi_int>& wi, std::vector<casadi_int>& wj,
321  std::vector<casadi_int>& queue, const std::vector<casadi_int>& imatch,
322  const std::vector<casadi_int>& jmatch, casadi_int mark) const;
323 
331  static void matched(
332  casadi_int n, const std::vector<casadi_int>& wj, const std::vector<casadi_int>& imatch,
333  std::vector<casadi_int>& p, std::vector<casadi_int>& q, std::vector<casadi_int>& cc,
334  std::vector<casadi_int>& rr, casadi_int set, casadi_int mark);
335 
343  static void unmatched(casadi_int m, const std::vector<casadi_int>& wi,
344  std::vector<casadi_int>& p, std::vector<casadi_int>& rr, casadi_int set);
345 
353  static casadi_int rprune(casadi_int i, casadi_int j, double aij, void *other);
354 
362  static casadi_int drop(casadi_int (*fkeep)(casadi_int, casadi_int, double, void *), void *other,
363  casadi_int nrow, casadi_int ncol,
364  std::vector<casadi_int>& colind, std::vector<casadi_int>& row);
365 
367  casadi_int btf(std::vector<casadi_int>& rowperm, std::vector<casadi_int>& colperm,
368  std::vector<casadi_int>& rowblock, std::vector<casadi_int>& colblock,
369  std::vector<casadi_int>& coarse_rowblock,
370  std::vector<casadi_int>& coarse_colblock) const {
371  T()->dmperm(colperm, rowperm, colblock, rowblock,
372  coarse_colblock, coarse_rowblock);
373  return rowblock.size()-1;
374  }
375 
377  const Btf& btf() const;
378 
386  void dmperm(std::vector<casadi_int>& rowperm, std::vector<casadi_int>& colperm,
387  std::vector<casadi_int>& rowblock, std::vector<casadi_int>& colblock,
388  std::vector<casadi_int>& coarse_rowblock,
389  std::vector<casadi_int>& coarse_colblock) const;
390 
398  void maxtrans(std::vector<casadi_int>& imatch,
399  std::vector<casadi_int>& jmatch, Sparsity& trans, casadi_int seed) const;
400 
408  void augment(casadi_int k, std::vector<casadi_int>& jmatch,
409  casadi_int *cheap, std::vector<casadi_int>& w, casadi_int *js,
410  casadi_int *is, casadi_int *ps) const;
411 
421  static std::vector<casadi_int> randperm(casadi_int n, casadi_int seed);
422 
426  static std::vector<casadi_int> invertPermutation(const std::vector<casadi_int>& p);
427 
429  Sparsity permute(const std::vector<casadi_int>& pinv,
430  const std::vector<casadi_int>& q, casadi_int values) const;
431 
439  void permute(const std::vector<casadi_int>& pinv,
440  const std::vector<casadi_int>& q, casadi_int values,
441  std::vector<casadi_int>& colind_C,
442  std::vector<casadi_int>& row_C) const;
443 
451  static casadi_int wclear(casadi_int mark, casadi_int lemax, casadi_int *w, casadi_int n);
452 
460  static casadi_int diag(casadi_int i, casadi_int j, double aij, void *other);
461 
469  Sparsity multiply(const Sparsity& B) const;
470 
478  casadi_int scatter(casadi_int j, std::vector<casadi_int>& w, casadi_int mark,
479  casadi_int* Ci, casadi_int nz) const;
480 
482  std::vector<casadi_int> get_row() const;
483 
485  std::vector<casadi_int> get_colind() const;
486 
488  std::vector<casadi_int> get_col() const;
489 
491  Sparsity _resize(casadi_int nrow, casadi_int ncol) const;
492 
494  Sparsity _reshape(casadi_int nrow, casadi_int ncol) const;
495 
497  casadi_int numel() const;
498 
500  casadi_int nnz_lower(bool strictly=false) const;
501 
503  casadi_int nnz_upper(bool strictly=false) const;
504 
506  casadi_int nnz_diag() const;
507 
511  casadi_int bw_upper() const;
512 
516  casadi_int bw_lower() const;
517 
519  std::pair<casadi_int, casadi_int> size() const;
520 
522  bool is_scalar(bool scalar_and_dense) const;
523 
530  bool is_empty(bool both=false) const;
531 
533  bool is_dense() const;
534 
538  bool is_row() const;
539 
543  bool is_column() const;
544 
548  bool is_vector() const;
549 
551  bool is_diag() const;
552 
554  bool is_square() const;
555 
559  bool is_permutation() const;
560 
564  bool is_selection(bool allow_empty=false) const;
565 
569  bool is_orthonormal(bool allow_empty=false) const;
570 
574  bool is_orthonormal_rows(bool allow_empty=false) const;
575 
579  bool is_orthonormal_columns(bool allow_empty=false) const;
580 
582  bool is_symmetric() const;
583 
585  bool is_tril(bool strictly) const;
586 
588  bool is_triu(bool strictly) const;
589 
591  Sparsity _triu(bool includeDiagonal) const;
592 
594  Sparsity _tril(bool includeDiagonal) const;
595 
597  std::vector<casadi_int> get_lower() const;
598 
600  std::vector<casadi_int> get_upper() const;
601 
603  std::string dim(bool with_nz=false) const;
604 
606  std::string repr_el(casadi_int k) const;
607 
609  Sparsity _mtimes(const Sparsity& y) const;
610 
613  Sparsity combine(const Sparsity& y, bool f0x_is_zero, bool function0_is_zero,
614  std::vector<unsigned char>& mapping) const;
615  Sparsity combine(const Sparsity& y, bool f0x_is_zero, bool function0_is_zero) const;
616 
617  template<bool with_mapping>
618  Sparsity combineGen1(const Sparsity& y, bool f0x_is_zero, bool function0_is_zero,
619  std::vector<unsigned char>& mapping) const;
620 
621  template<bool with_mapping, bool f0x_is_zero, bool function0_is_zero>
622  Sparsity combineGen(const Sparsity& y, std::vector<unsigned char>& mapping) const;
624 
626  bool is_subset(const Sparsity& rhs) const;
627 
629  Sparsity pattern_inverse() const;
630 
632  bool is_equal(const Sparsity& y) const;
633 
635  bool is_equal(casadi_int y_nrow, casadi_int y_ncol, const std::vector<casadi_int>& y_colind,
636  const std::vector<casadi_int>& y_row) const;
637 
639  bool is_equal(casadi_int y_nrow, casadi_int y_ncol,
640  const casadi_int* y_colind, const casadi_int* y_row) const;
641 
643  bool is_stacked(const Sparsity& y, casadi_int n) const;
644 
646  Sparsity _enlargeRows(casadi_int nrow, const std::vector<casadi_int>& rr, bool ind1) const;
647 
649  Sparsity _enlargeColumns(casadi_int ncol, const std::vector<casadi_int>& cc, bool ind1) const;
650 
652  Sparsity makeDense(std::vector<casadi_int>& mapping) const;
653 
655  Sparsity _erase(const std::vector<casadi_int>& rr, const std::vector<casadi_int>& cc,
656  bool ind1, std::vector<casadi_int>& mapping) const;
657 
659  Sparsity _erase(const std::vector<casadi_int>& rr, bool ind1,
660  std::vector<casadi_int>& mapping) const;
661 
663  Sparsity _appendVector(const SparsityInternal& sp) const;
664 
666  Sparsity _appendColumns(const SparsityInternal& sp) const;
667 
674  Sparsity sub(const std::vector<casadi_int>& rr, const std::vector<casadi_int>& cc,
675  std::vector<casadi_int>& mapping, bool ind1) const;
676 
683  Sparsity sub(const std::vector<casadi_int>& rr, const SparsityInternal& sp,
684  std::vector<casadi_int>& mapping, bool ind1) const;
685 
687  casadi_int get_nz(casadi_int rr, casadi_int cc) const;
688 
690  std::vector<casadi_int> get_nz(const std::vector<casadi_int>& rr,
691  const std::vector<casadi_int>& cc) const;
692 
694  void get_nz(std::vector<casadi_int>& indices) const;
695 
697  bool rowsSequential(bool strictly) const;
698 
705  Sparsity _removeDuplicates(std::vector<casadi_int>& mapping) const;
706 
708  void find(std::vector<casadi_int>& loc, bool ind1) const;
709 
711  std::size_t hash() const;
712 
714  std::string class_name() const override {return "SparsityInternal";}
715 
717  void disp(std::ostream& stream, bool more) const override;
718 
725  Sparsity uni_coloring(const Sparsity& AT, casadi_int cutoff) const;
726 
732  Sparsity star_coloring_new(std::vector<casadi_int>& which_color, const Dict& opts) const;
733 
739  Sparsity star_coloring(casadi_int ordering, casadi_int cutoff) const;
740 
746  Sparsity star_coloring2(casadi_int ordering, casadi_int cutoff) const;
747 
749  std::vector<casadi_int> largest_first() const;
750 
752  Sparsity pmult(const std::vector<casadi_int>& p, bool permute_rows=true, bool permute_cols=true,
753  bool invert_permutation=false) const;
754 
758  void spy(std::ostream &stream) const;
759 
761  void spy_matlab(const std::string& mfile) const;
762 
766  void export_code(const std::string& lang, std::ostream &stream,
767  const Dict& options) const;
768 
770  void spsolve(bvec_t* X, bvec_t* B, bool tr) const;
771 };
772 
773 } // namespace casadi
775 
776 #endif // CASADI_SPARSITY_INTERNAL_HPP
casadi
The casadi namespace.
Definition: archiver.hpp:32
casadi::Dict
GenericType::Dict Dict
C++ equivalent of Python's dict or MATLAB's struct.
Definition: generic_type.hpp:258