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