matrix.h
1 #ifndef SYMENGINE_MATRIX_H
2 #define SYMENGINE_MATRIX_H
3 
4 #include <symengine/basic.h>
5 #include <symengine/sets.h>
6 
7 namespace SymEngine
8 {
9 
10 // Base class for matrices
12 {
13 public:
14  virtual ~MatrixBase(){};
15 
16  bool is_square() const
17  {
18  return ncols() == nrows();
19  }
20 
21  // Below methods should be implemented by the derived classes. If not
22  // applicable, raise an exception
23 
24  // Get the # of rows and # of columns
25  virtual unsigned nrows() const = 0;
26  virtual unsigned ncols() const = 0;
27  virtual bool eq(const MatrixBase &other) const;
28 
29  // Get and set elements
30  virtual RCP<const Basic> get(unsigned i, unsigned j) const = 0;
31  virtual void set(unsigned i, unsigned j, const RCP<const Basic> &e) = 0;
32 
33  // Print Matrix, very mundane version, should be overriden derived
34  // class if better printing is available
35  virtual std::string __str__() const;
36 
37  virtual unsigned rank() const = 0;
38  virtual RCP<const Basic> det() const = 0;
39  virtual void inv(MatrixBase &result) const = 0;
40 
41  // Matrix addition
42  virtual void add_matrix(const MatrixBase &other,
43  MatrixBase &result) const = 0;
44 
45  // Matrix Multiplication
46  virtual void mul_matrix(const MatrixBase &other,
47  MatrixBase &result) const = 0;
48 
49  // Matrix elementwise Multiplication
50  virtual void elementwise_mul_matrix(const MatrixBase &other,
51  MatrixBase &result) const = 0;
52 
53  // Add a scalar
54  virtual void add_scalar(const RCP<const Basic> &k,
55  MatrixBase &result) const = 0;
56 
57  // Multiply by a scalar
58  virtual void mul_scalar(const RCP<const Basic> &k,
59  MatrixBase &result) const = 0;
60 
61  // Matrix conjugate
62  virtual void conjugate(MatrixBase &result) const = 0;
63 
64  // Matrix transpose
65  virtual void transpose(MatrixBase &result) const = 0;
66 
67  // Matrix conjugate transpose
68  virtual void conjugate_transpose(MatrixBase &result) const = 0;
69 
70  // Extract out a submatrix
71  virtual void submatrix(MatrixBase &result, unsigned row_start,
72  unsigned col_start, unsigned row_end,
73  unsigned col_end, unsigned row_step = 1,
74  unsigned col_step = 1) const = 0;
75  // LU factorization
76  virtual void LU(MatrixBase &L, MatrixBase &U) const = 0;
77 
78  // LDL factorization
79  virtual void LDL(MatrixBase &L, MatrixBase &D) const = 0;
80 
81  // Fraction free LU factorization
82  virtual void FFLU(MatrixBase &LU) const = 0;
83 
84  // Fraction free LDU factorization
85  virtual void FFLDU(MatrixBase &L, MatrixBase &D, MatrixBase &U) const = 0;
86 
87  // QR factorization
88  virtual void QR(MatrixBase &Q, MatrixBase &R) const = 0;
89 
90  // Cholesky decomposition
91  virtual void cholesky(MatrixBase &L) const = 0;
92 
93  // Solve Ax = b using LU factorization
94  virtual void LU_solve(const MatrixBase &b, MatrixBase &x) const = 0;
95 };
96 
98 
99 class CSRMatrix;
100 
101 // ----------------------------- Dense Matrix --------------------------------//
102 class DenseMatrix : public MatrixBase
103 {
104 public:
105  // Constructors
106  DenseMatrix();
107  DenseMatrix(const DenseMatrix &) = default;
108  DenseMatrix(unsigned row, unsigned col);
109  DenseMatrix(unsigned row, unsigned col, const vec_basic &l);
110  DenseMatrix(const vec_basic &column_elements);
111  DenseMatrix &operator=(const DenseMatrix &other) = default;
112  // Resize
113  void resize(unsigned i, unsigned j);
114 
115  // Should implement all the virtual methods from MatrixBase
116  // and throw an exception if a method is not applicable.
117 
118  // Get and set elements
119  virtual RCP<const Basic> get(unsigned i, unsigned j) const;
120  virtual void set(unsigned i, unsigned j, const RCP<const Basic> &e);
121  virtual vec_basic as_vec_basic() const;
122 
123  virtual unsigned nrows() const
124  {
125  return row_;
126  }
127  virtual unsigned ncols() const
128  {
129  return col_;
130  }
131 
132  virtual bool is_lower() const;
133  virtual bool is_upper() const;
134  virtual tribool is_zero() const;
135  virtual tribool is_diagonal() const;
136  virtual tribool is_real(const Assumptions *assumptions = nullptr) const;
137  virtual tribool is_symmetric() const;
138  virtual tribool is_hermitian() const;
139  virtual tribool is_weakly_diagonally_dominant() const;
140  virtual tribool is_strictly_diagonally_dominant() const;
141  virtual tribool is_positive_definite() const;
142  virtual tribool is_negative_definite() const;
143 
144  RCP<const Basic> trace() const;
145  virtual unsigned rank() const;
146  virtual RCP<const Basic> det() const;
147  virtual void inv(MatrixBase &result) const;
148 
149  // Matrix addition
150  virtual void add_matrix(const MatrixBase &other, MatrixBase &result) const;
151 
152  // Matrix multiplication
153  virtual void mul_matrix(const MatrixBase &other, MatrixBase &result) const;
154 
155  // Matrix elementwise Multiplication
156  virtual void elementwise_mul_matrix(const MatrixBase &other,
157  MatrixBase &result) const;
158 
159  // Add a scalar
160  virtual void add_scalar(const RCP<const Basic> &k,
161  MatrixBase &result) const;
162 
163  // Multiply by a scalar
164  virtual void mul_scalar(const RCP<const Basic> &k,
165  MatrixBase &result) const;
166 
167  // Matrix conjugate
168  virtual void conjugate(MatrixBase &result) const;
169 
170  // Matrix transpose
171  virtual void transpose(MatrixBase &result) const;
172 
173  // Matrix conjugate transpose
174  virtual void conjugate_transpose(MatrixBase &result) const;
175 
176  // Extract out a submatrix
177  virtual void submatrix(MatrixBase &result, unsigned row_start,
178  unsigned col_start, unsigned row_end,
179  unsigned col_end, unsigned row_step = 1,
180  unsigned col_step = 1) const;
181 
182  // LU factorization
183  virtual void LU(MatrixBase &L, MatrixBase &U) const;
184 
185  // LDL factorization
186  virtual void LDL(MatrixBase &L, MatrixBase &D) const;
187 
188  // Solve Ax = b using LU factorization
189  virtual void LU_solve(const MatrixBase &b, MatrixBase &x) const;
190 
191  // Fraction free LU factorization
192  virtual void FFLU(MatrixBase &LU) const;
193 
194  // Fraction free LDU factorization
195  virtual void FFLDU(MatrixBase &L, MatrixBase &D, MatrixBase &U) const;
196 
197  // QR factorization
198  virtual void QR(MatrixBase &Q, MatrixBase &R) const;
199 
200  // Cholesky decomposition
201  virtual void cholesky(MatrixBase &L) const;
202 
203  // Return the Jacobian of the matrix
204  friend void jacobian(const DenseMatrix &A, const DenseMatrix &x,
205  DenseMatrix &result, bool diff_cache);
206  // Return the Jacobian of the matrix using sdiff
207  friend void sjacobian(const DenseMatrix &A, const DenseMatrix &x,
208  DenseMatrix &result, bool diff_cache);
209 
210  // Differentiate the matrix element-wise
211  friend void diff(const DenseMatrix &A, const RCP<const Symbol> &x,
212  DenseMatrix &result, bool diff_cache);
213  // Differentiate the matrix element-wise using SymPy compatible diff
214  friend void sdiff(const DenseMatrix &A, const RCP<const Basic> &x,
215  DenseMatrix &result, bool diff_cache);
216 
217  // Friend functions related to Matrix Operations
218  friend void add_dense_dense(const DenseMatrix &A, const DenseMatrix &B,
219  DenseMatrix &C);
220  friend void add_dense_scalar(const DenseMatrix &A,
221  const RCP<const Basic> &k, DenseMatrix &B);
222  friend void mul_dense_dense(const DenseMatrix &A, const DenseMatrix &B,
223  DenseMatrix &C);
224  friend void elementwise_mul_dense_dense(const DenseMatrix &A,
225  const DenseMatrix &B,
226  DenseMatrix &C);
227  friend void mul_dense_scalar(const DenseMatrix &A,
228  const RCP<const Basic> &k, DenseMatrix &C);
229  friend void conjugate_dense(const DenseMatrix &A, DenseMatrix &B);
230  friend void transpose_dense(const DenseMatrix &A, DenseMatrix &B);
231  friend void conjugate_transpose_dense(const DenseMatrix &A, DenseMatrix &B);
232  friend void submatrix_dense(const DenseMatrix &A, DenseMatrix &B,
233  unsigned row_start, unsigned col_start,
234  unsigned row_end, unsigned col_end,
235  unsigned row_step, unsigned col_step);
236  void row_join(const DenseMatrix &B);
237  void col_join(const DenseMatrix &B);
238  void row_insert(const DenseMatrix &B, unsigned pos);
239  void col_insert(const DenseMatrix &B, unsigned pos);
240  void row_del(unsigned k);
241  void col_del(unsigned k);
242 
243  // Row operations
244  friend void row_exchange_dense(DenseMatrix &A, unsigned i, unsigned j);
245  friend void row_mul_scalar_dense(DenseMatrix &A, unsigned i,
246  RCP<const Basic> &c);
247  friend void row_add_row_dense(DenseMatrix &A, unsigned i, unsigned j,
248  RCP<const Basic> &c);
249  friend void permuteFwd(DenseMatrix &A, permutelist &pl);
250 
251  // Column operations
252  friend void column_exchange_dense(DenseMatrix &A, unsigned i, unsigned j);
253 
254  // Gaussian elimination
255  friend void pivoted_gaussian_elimination(const DenseMatrix &A,
256  DenseMatrix &B,
257  permutelist &pivotlist);
258  friend void fraction_free_gaussian_elimination(const DenseMatrix &A,
259  DenseMatrix &B);
260  friend void pivoted_fraction_free_gaussian_elimination(
261  const DenseMatrix &A, DenseMatrix &B, permutelist &pivotlist);
262  friend void pivoted_gauss_jordan_elimination(const DenseMatrix &A,
263  DenseMatrix &B,
264  permutelist &pivotlist);
265  friend void fraction_free_gauss_jordan_elimination(const DenseMatrix &A,
266  DenseMatrix &B);
267  friend void pivoted_fraction_free_gauss_jordan_elimination(
268  const DenseMatrix &A, DenseMatrix &B, permutelist &pivotlist);
269  friend unsigned pivot(DenseMatrix &B, unsigned r, unsigned c);
270 
271  friend void reduced_row_echelon_form(const DenseMatrix &A, DenseMatrix &B,
272  vec_uint &pivot_cols,
273  bool normalize_last);
274 
275  // Ax = b
276  friend void diagonal_solve(const DenseMatrix &A, const DenseMatrix &b,
277  DenseMatrix &x);
278  friend void back_substitution(const DenseMatrix &U, const DenseMatrix &b,
279  DenseMatrix &x);
280  friend void forward_substitution(const DenseMatrix &A, const DenseMatrix &b,
281  DenseMatrix &x);
282  friend void fraction_free_gaussian_elimination_solve(const DenseMatrix &A,
283  const DenseMatrix &b,
284  DenseMatrix &x);
285  friend void fraction_free_gauss_jordan_solve(const DenseMatrix &A,
286  const DenseMatrix &b,
287  DenseMatrix &x, bool pivot);
288 
289  // Matrix Decomposition
290  friend void fraction_free_LU(const DenseMatrix &A, DenseMatrix &LU);
291  friend void LU(const DenseMatrix &A, DenseMatrix &L, DenseMatrix &U);
292  friend void pivoted_LU(const DenseMatrix &A, DenseMatrix &LU,
293  permutelist &pl);
294  friend void pivoted_LU(const DenseMatrix &A, DenseMatrix &L, DenseMatrix &U,
295  permutelist &pl);
296  friend void fraction_free_LDU(const DenseMatrix &A, DenseMatrix &L,
297  DenseMatrix &D, DenseMatrix &U);
298  friend void QR(const DenseMatrix &A, DenseMatrix &Q, DenseMatrix &R);
299  friend void LDL(const DenseMatrix &A, DenseMatrix &L, DenseMatrix &D);
300  friend void cholesky(const DenseMatrix &A, DenseMatrix &L);
301 
302  // Matrix queries
303  friend bool is_symmetric_dense(const DenseMatrix &A);
304 
305  // Determinant
306  friend RCP<const Basic> det_bareis(const DenseMatrix &A);
307  friend void berkowitz(const DenseMatrix &A,
308  std::vector<DenseMatrix> &polys);
309 
310  // Inverse
311  friend void inverse_fraction_free_LU(const DenseMatrix &A, DenseMatrix &B);
312  friend void inverse_LU(const DenseMatrix &A, DenseMatrix &B);
313  friend void inverse_pivoted_LU(const DenseMatrix &A, DenseMatrix &B);
314  friend void inverse_gauss_jordan(const DenseMatrix &A, DenseMatrix &B);
315 
316  // Vector-specific methods
317  friend void dot(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &C);
318  friend void cross(const DenseMatrix &A, const DenseMatrix &B,
319  DenseMatrix &C);
320 
321  // NumPy-like functions
322  friend void eye(DenseMatrix &A, int k);
323  friend void diag(DenseMatrix &A, vec_basic &v, int k);
324  friend void ones(DenseMatrix &A);
325  friend void zeros(DenseMatrix &A);
326 
327  friend CSRMatrix;
328 
329 private:
330  // Matrix elements are stored in row-major order
331  vec_basic m_;
332  // Stores the dimension of the Matrix
333  unsigned row_;
334  unsigned col_;
335 
336  tribool shortcut_to_posdef() const;
337  tribool is_positive_definite_GE();
338 };
339 
340 // ----------------------------- Sparse Matrices -----------------------------//
341 class CSRMatrix : public MatrixBase
342 {
343 public:
344  CSRMatrix();
345  CSRMatrix(unsigned row, unsigned col);
346  CSRMatrix(unsigned row, unsigned col, const std::vector<unsigned> &p,
347  const std::vector<unsigned> &j, const vec_basic &x);
348  CSRMatrix(unsigned row, unsigned col, std::vector<unsigned> &&p,
350  CSRMatrix &operator=(CSRMatrix &&other);
351  CSRMatrix(const CSRMatrix &) = default;
353  as_vectors() const;
354 
355  bool is_canonical() const;
356 
357  virtual bool eq(const MatrixBase &other) const;
358 
359  // Get and set elements
360  virtual RCP<const Basic> get(unsigned i, unsigned j) const;
361  virtual void set(unsigned i, unsigned j, const RCP<const Basic> &e);
362 
363  virtual unsigned nrows() const
364  {
365  return row_;
366  }
367  virtual unsigned ncols() const
368  {
369  return col_;
370  }
371  virtual unsigned rank() const;
372  virtual RCP<const Basic> det() const;
373  virtual void inv(MatrixBase &result) const;
374 
375  // Matrix addition
376  virtual void add_matrix(const MatrixBase &other, MatrixBase &result) const;
377 
378  // Matrix Multiplication
379  virtual void mul_matrix(const MatrixBase &other, MatrixBase &result) const;
380 
381  // Matrix elementwise Multiplication
382  virtual void elementwise_mul_matrix(const MatrixBase &other,
383  MatrixBase &result) const;
384 
385  // Add a scalar
386  virtual void add_scalar(const RCP<const Basic> &k,
387  MatrixBase &result) const;
388 
389  // Multiply by a scalar
390  virtual void mul_scalar(const RCP<const Basic> &k,
391  MatrixBase &result) const;
392 
393  // Matrix conjugate
394  virtual void conjugate(MatrixBase &result) const;
395 
396  // Matrix transpose
397  virtual void transpose(MatrixBase &result) const;
398  CSRMatrix transpose(bool conjugate = false) const;
399 
400  // Matrix conjugate transpose
401  virtual void conjugate_transpose(MatrixBase &result) const;
402 
403  // Extract out a submatrix
404  virtual void submatrix(MatrixBase &result, unsigned row_start,
405  unsigned col_start, unsigned row_end,
406  unsigned col_end, unsigned row_step = 1,
407  unsigned col_step = 1) const;
408 
409  // LU factorization
410  virtual void LU(MatrixBase &L, MatrixBase &U) const;
411 
412  // LDL factorization
413  virtual void LDL(MatrixBase &L, MatrixBase &D) const;
414 
415  // Solve Ax = b using LU factorization
416  virtual void LU_solve(const MatrixBase &b, MatrixBase &x) const;
417 
418  // Fraction free LU factorization
419  virtual void FFLU(MatrixBase &LU) const;
420 
421  // Fraction free LDU factorization
422  virtual void FFLDU(MatrixBase &L, MatrixBase &D, MatrixBase &U) const;
423 
424  // QR factorization
425  virtual void QR(MatrixBase &Q, MatrixBase &R) const;
426 
427  // Cholesky decomposition
428  virtual void cholesky(MatrixBase &L) const;
429 
430  static void csr_sum_duplicates(std::vector<unsigned> &p_,
432  unsigned row_);
433 
434  static void csr_sort_indices(std::vector<unsigned> &p_,
436  unsigned row_);
437 
438  static bool csr_has_sorted_indices(const std::vector<unsigned> &p_,
439  const std::vector<unsigned> &j_,
440  unsigned row_);
441 
442  static bool csr_has_duplicates(const std::vector<unsigned> &p_,
443  const std::vector<unsigned> &j_,
444  unsigned row_);
445 
446  static bool csr_has_canonical_format(const std::vector<unsigned> &p_,
447  const std::vector<unsigned> &j_,
448  unsigned row_);
449 
450  static CSRMatrix from_coo(unsigned row, unsigned col,
451  const std::vector<unsigned> &i,
452  const std::vector<unsigned> &j,
453  const vec_basic &x);
454  static CSRMatrix jacobian(const vec_basic &exprs, const vec_sym &x,
455  bool diff_cache = true);
456  static CSRMatrix jacobian(const DenseMatrix &A, const DenseMatrix &x,
457  bool diff_cache = true);
458 
459  friend void csr_matmat_pass1(const CSRMatrix &A, const CSRMatrix &B,
460  CSRMatrix &C);
461  friend void csr_matmat_pass2(const CSRMatrix &A, const CSRMatrix &B,
462  CSRMatrix &C);
463  friend void csr_diagonal(const CSRMatrix &A, DenseMatrix &D);
464  friend void csr_scale_rows(CSRMatrix &A, const DenseMatrix &X);
465  friend void csr_scale_columns(CSRMatrix &A, const DenseMatrix &X);
466 
467  friend void csr_binop_csr_canonical(
468  const CSRMatrix &A, const CSRMatrix &B, CSRMatrix &C,
469  RCP<const Basic> (&bin_op)(const RCP<const Basic> &,
470  const RCP<const Basic> &));
471 
472 private:
475  vec_basic x_;
476  // Stores the dimension of the Matrix
477  unsigned row_;
478  unsigned col_;
479 };
480 
481 // Return the Jacobian of the matrix
482 void jacobian(const DenseMatrix &A, const DenseMatrix &x, DenseMatrix &result,
483  bool diff_cache = true);
484 // Return the Jacobian of the matrix using sdiff
485 void sjacobian(const DenseMatrix &A, const DenseMatrix &x, DenseMatrix &result,
486  bool diff_cache = true);
487 
488 // Differentiate all the elements
489 void diff(const DenseMatrix &A, const RCP<const Symbol> &x, DenseMatrix &result,
490  bool diff_cache = true);
491 // Differentiate all the elements using SymPy compatible diff
492 void sdiff(const DenseMatrix &A, const RCP<const Basic> &x, DenseMatrix &result,
493  bool diff_cache = true);
494 
495 // Get submatrix from a DenseMatrix
496 void submatrix_dense(const DenseMatrix &A, DenseMatrix &B, unsigned row_start,
497  unsigned col_start, unsigned row_end, unsigned col_end,
498  unsigned row_step = 1, unsigned col_step = 1);
499 
500 // Row operations
501 void row_exchange_dense(DenseMatrix &A, unsigned i, unsigned j);
502 void row_mul_scalar_dense(DenseMatrix &A, unsigned i, RCP<const Basic> &c);
503 void row_add_row_dense(DenseMatrix &A, unsigned i, unsigned j,
504  RCP<const Basic> &c);
505 
506 // Column operations
507 void column_exchange_dense(DenseMatrix &A, unsigned i, unsigned j);
508 
509 // Vector-specific methods
510 void dot(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &C);
511 void cross(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &C);
512 
513 // Matrix Factorization
514 void LU(const DenseMatrix &A, DenseMatrix &L, DenseMatrix &U);
515 void LDL(const DenseMatrix &A, DenseMatrix &L, DenseMatrix &D);
516 void QR(const DenseMatrix &A, DenseMatrix &Q, DenseMatrix &R);
517 void cholesky(const DenseMatrix &A, DenseMatrix &L);
518 
519 // Inverse
520 void inverse_fraction_free_LU(const DenseMatrix &A, DenseMatrix &B);
521 
522 void inverse_gauss_jordan(const DenseMatrix &A, DenseMatrix &B);
523 
524 // Solving Ax = b
525 void fraction_free_LU_solve(const DenseMatrix &A, const DenseMatrix &b,
526  DenseMatrix &x);
527 
528 void fraction_free_gauss_jordan_solve(const DenseMatrix &A,
529  const DenseMatrix &b, DenseMatrix &x,
530  bool pivot = true);
531 
532 void LU_solve(const DenseMatrix &A, const DenseMatrix &b, DenseMatrix &x);
533 void pivoted_LU_solve(const DenseMatrix &A, const DenseMatrix &b,
534  DenseMatrix &x);
535 
536 void LDL_solve(const DenseMatrix &A, const DenseMatrix &b, DenseMatrix &x);
537 
538 // Determinant
539 RCP<const Basic> det_berkowitz(const DenseMatrix &A);
540 
541 // Characteristic polynomial: Only the coefficients of monomials in decreasing
542 // order of monomial powers is returned, i.e. if `B = transpose([1, -2, 3])`
543 // then the corresponding polynomial is `x**2 - 2x + 3`.
544 void char_poly(const DenseMatrix &A, DenseMatrix &B);
545 
546 // returns a finiteset of eigenvalues of a matrix
547 RCP<const Set> eigen_values(const DenseMatrix &A);
548 
549 // Mimic `eye` function in NumPy
550 void eye(DenseMatrix &A, int k = 0);
551 
552 // Create diagonal matrices directly
553 void diag(DenseMatrix &A, vec_basic &v, int k = 0);
554 
555 // Create a matrix filled with ones
556 void ones(DenseMatrix &A);
557 
558 // Create a matrix filled with zeros
559 void zeros(DenseMatrix &A);
560 
561 // Reduced row echelon form and returns the cols with pivots
562 void reduced_row_echelon_form(const DenseMatrix &A, DenseMatrix &B,
563  vec_uint &pivot_cols,
564  bool normalize_last = false);
565 
566 // Returns true if `b` is exactly the type T.
567 // Here T can be a DenseMatrix, CSRMatrix, etc.
568 template <class T>
569 inline bool is_a(const MatrixBase &b)
570 {
571  return typeid(T) == typeid(b);
572 }
573 
574 // Test two matrices for equality
575 inline bool operator==(const SymEngine::MatrixBase &lhs,
576  const SymEngine::MatrixBase &rhs)
577 {
578  return lhs.eq(rhs);
579 }
580 
581 // Test two matrices for equality
582 inline bool operator!=(const SymEngine::MatrixBase &lhs,
583  const SymEngine::MatrixBase &rhs)
584 {
585  return not lhs.eq(rhs);
586 }
587 
588 } // namespace SymEngine
589 
590 // Print Matrix
592  const SymEngine::MatrixBase &A)
593 {
594  return out << A.__str__();
595 }
596 
597 #endif
The base class for SymEngine.
Main namespace for SymEngine package.
Definition: add.cpp:19
bool is_a(const Basic &b)
Templatised version to check is_a type.
Definition: basic-inl.h:36
std::ostream & operator<<(std::ostream &out, const SymEngine::Basic &p)
<< Operator
Definition: basic-inl.h:95
RCP< const Basic > conjugate(const RCP< const Basic > &arg)
Canonicalize Conjugate.
Definition: functions.cpp:106
T operator!=(T... args)