Loading...
Searching...
No Matches
dense_matrix.cpp
1#include <symengine/matrix.h>
2#include <symengine/number.h>
3#include <symengine/add.h>
5#include <symengine/pow.h>
6#include <symengine/subs.h>
7#include <symengine/symengine_exception.h>
8#include <symengine/polys/uexprpoly.h>
9#include <symengine/solve.h>
10#include <symengine/test_visitors.h>
11
12namespace SymEngine
13{
14
15// Constructors
16DenseMatrix::DenseMatrix() {}
17
18DenseMatrix::DenseMatrix(unsigned row, unsigned col) : row_(row), col_(col)
19{
20 m_ = std::vector<RCP<const Basic>>(row * col);
21}
22
23DenseMatrix::DenseMatrix(unsigned row, unsigned col, const vec_basic &l)
24 : m_{l}, row_(row), col_(col){SYMENGINE_ASSERT(m_.size() == row * col)}
25
26 DenseMatrix::DenseMatrix(const vec_basic &column_elements)
27 : m_(column_elements), row_(static_cast<unsigned>(column_elements.size())),
28 col_(1)
29{
30}
31
32// Resize the Matrix
33void DenseMatrix::resize(unsigned row, unsigned col)
34{
35 row_ = row;
36 col_ = col;
37 m_.resize(row * col);
38}
39
40// Get and set elements
41RCP<const Basic> DenseMatrix::get(unsigned i, unsigned j) const
42{
43 SYMENGINE_ASSERT(i < row_ and j < col_);
44 return m_[i * col_ + j];
45}
46
47void DenseMatrix::set(unsigned i, unsigned j, const RCP<const Basic> &e)
48{
49 SYMENGINE_ASSERT(i < row_ and j < col_);
50 m_[i * col_ + j] = e;
51}
52
53// Flat (row major) view of Matrix
54vec_basic DenseMatrix::as_vec_basic() const
55{
56 return m_;
57}
58
59RCP<const Basic> DenseMatrix::trace() const
60{
61 SYMENGINE_ASSERT(row_ == col_);
62 unsigned offset = 0;
63 vec_basic diag;
64 for (unsigned i = 0; i < row_; i++) {
65 diag.push_back(m_[offset]);
66 offset += row_ + 1;
67 }
68 auto sum = add(diag);
69 return sum;
70}
71
72unsigned DenseMatrix::rank() const
73{
74 throw NotImplementedError("Not Implemented");
75}
76
77bool DenseMatrix::is_lower() const
78{
79 auto A = *this;
80 unsigned n = A.nrows();
81 for (unsigned i = 1; i < n; ++i) {
82 for (unsigned j = 0; j < i; ++j) {
83 if (not is_number_and_zero(*A.get(i, j))) {
84 return false;
85 }
86 }
87 }
88 return true;
89}
90
91bool DenseMatrix::is_upper() const
92{
93 auto A = *this;
94 unsigned n = A.nrows();
95 for (unsigned i = 0; i < n - 1; ++i) {
96 for (unsigned j = i + 1; j < n; ++j) {
97 if (not is_number_and_zero(*A.get(i, j))) {
98 return false;
99 }
100 }
101 }
102 return true;
103}
104
105tribool DenseMatrix::is_zero() const
106{
107 auto A = *this;
108 tribool cur = tribool::tritrue;
109 for (auto &e : m_) {
110 cur = and_tribool(cur, SymEngine::is_zero(*e));
111 if (is_false(cur)) {
112 return cur;
113 }
114 }
115 return cur;
116}
117
118tribool DenseMatrix::is_diagonal() const
119{
120 auto A = *this;
121 if (not A.is_square()) {
122 return tribool::trifalse;
123 }
124 unsigned ncols = A.ncols();
125 unsigned offset;
126 tribool cur = tribool::tritrue;
127 for (unsigned i = 0; i < ncols; i++) {
128 offset = i * ncols;
129 for (unsigned j = 0; j < ncols; j++) {
130 if (j != i) {
131 auto &e = m_[offset];
132 cur = and_tribool(cur, SymEngine::is_zero(*e));
133 if (is_false(cur)) {
134 return cur;
135 }
136 }
137 offset++;
138 }
139 }
140 return cur;
141}
142
143tribool DenseMatrix::is_real(const Assumptions *assumptions) const
144{
145 RealVisitor visitor(assumptions);
146 tribool cur = tribool::tritrue;
147 for (auto &e : m_) {
148 cur = and_tribool(cur, visitor.apply(*e));
149 if (is_false(cur)) {
150 return cur;
151 }
152 }
153 return cur;
154}
155
156tribool DenseMatrix::is_symmetric() const
157{
158 auto A = *this;
159 if (not A.is_square()) {
160 return tribool::trifalse;
161 }
162 unsigned ncols = A.ncols();
163 tribool cur = tribool::tritrue;
164 for (unsigned i = 0; i < ncols; i++) {
165 for (unsigned j = 0; j <= i; j++) {
166 if (j != i) {
167 auto &e = m_[i * ncols + j];
168 auto &e2 = m_[j * ncols + i];
169 cur = and_tribool(cur, SymEngine::is_zero(*sub(e, e2)));
170 }
171 if (is_false(cur)) {
172 return cur;
173 }
174 }
175 }
176 return cur;
177}
178
179tribool DenseMatrix::is_hermitian() const
180{
181 auto A = *this;
182 if (not A.is_square()) {
183 return tribool::trifalse;
184 }
185 unsigned ncols = A.ncols();
186 tribool cur = tribool::tritrue;
187 for (unsigned i = 0; i < ncols; i++) {
188 for (unsigned j = 0; j <= i; j++) {
189 auto &e = m_[i * ncols + j];
190 if (j != i) {
191 auto &e2 = m_[j * ncols + i];
192 cur = and_tribool(
194 } else {
195 cur = and_tribool(cur, SymEngine::is_real(*e));
196 }
197 if (is_false(cur)) {
198 return cur;
199 }
200 }
201 }
202 return cur;
203}
204
205tribool DenseMatrix::is_weakly_diagonally_dominant() const
206{
207 auto A = *this;
208 if (not A.is_square()) {
209 return tribool::trifalse;
210 }
211
212 unsigned ncols = A.ncols();
213 RCP<const Basic> diag;
214 RCP<const Basic> sum;
215 tribool diagdom = tribool::tritrue;
216 for (unsigned i = 0; i < ncols; i++) {
217 sum = zero;
218 for (unsigned j = 0; j < ncols; j++) {
219 auto &e = m_[i * ncols + j];
220 if (i == j) {
221 diag = abs(e);
222 } else {
223 sum = add(sum, abs(e));
224 }
225 }
226 diagdom = and_tribool(diagdom, is_nonnegative(*sub(diag, sum)));
227 if (is_false(diagdom)) {
228 return diagdom;
229 }
230 }
231 return diagdom;
232}
233
234tribool DenseMatrix::is_strictly_diagonally_dominant() const
235{
236 auto A = *this;
237 if (not A.is_square()) {
238 return tribool::trifalse;
239 }
240
241 unsigned ncols = A.ncols();
242 RCP<const Basic> diag;
243 RCP<const Basic> sum;
244 tribool diagdom = tribool::tritrue;
245 for (unsigned i = 0; i < ncols; i++) {
246 sum = zero;
247 for (unsigned j = 0; j < ncols; j++) {
248 auto &e = m_[i * ncols + j];
249 if (i == j) {
250 diag = abs(e);
251 } else {
252 sum = add(sum, abs(e));
253 }
254 }
255 diagdom = and_tribool(diagdom, is_positive(*sub(diag, sum)));
256 if (is_false(diagdom)) {
257 return diagdom;
258 }
259 }
260 return diagdom;
261}
262
263tribool DenseMatrix::shortcut_to_posdef() const
264{
265 tribool is_diagonal_positive = tribool::tritrue;
266 unsigned offset = 0;
267 for (unsigned i = 0; i < row_; i++) {
268 is_diagonal_positive
269 = and_tribool(is_diagonal_positive, is_positive(*m_[offset]));
270 if (is_false(is_diagonal_positive))
271 return is_diagonal_positive;
272 offset += row_ + 1;
273 }
274 if (is_true(and_tribool(is_diagonal_positive,
275 this->is_strictly_diagonally_dominant())))
276 return tribool::tritrue;
277 return tribool::indeterminate;
278}
279
280tribool DenseMatrix::is_positive_definite_GE()
281{
282 auto size = row_;
283 for (unsigned i = 0; i < size; i++) {
284 auto ispos = is_positive(*m_[i * size + i]);
285 if (!is_true(ispos))
286 return ispos;
287 for (unsigned j = i + 1; j < size; j++) {
288 for (unsigned k = i + 1; k < size; k++) {
289 m_[j * size + k] = sub(mul(m_[i * size + i], m_[j * size + k]),
290 mul(m_[j * size + i], m_[i * size + k]));
291 }
292 }
293 }
294 return tribool::tritrue;
295}
296
297tribool DenseMatrix::is_positive_definite() const
298{
299
300 auto A = *this;
302 const DenseMatrix *H;
303 if (!is_true(A.is_hermitian())) {
304 if (!A.is_square())
305 return tribool::trifalse;
306 DenseMatrix tmp1 = DenseMatrix(A.row_, A.col_);
307 B = std::unique_ptr<DenseMatrix>(new DenseMatrix(A.row_, A.col_));
308 A.conjugate_transpose(tmp1);
309 add_dense_dense(A, tmp1, *B.get());
310 H = B.get();
311 } else {
312 H = this;
313 }
314
315 tribool shortcut = H->shortcut_to_posdef();
316 if (!is_indeterminate(shortcut))
317 return shortcut;
318
319 if (!B) {
320 B = std::unique_ptr<DenseMatrix>(new DenseMatrix(A));
321 }
322 return B->is_positive_definite_GE();
323}
324
325tribool DenseMatrix::is_negative_definite() const
326{
327 auto res = DenseMatrix(row_, col_);
328 mul_dense_scalar(*this, integer(-1), res);
329 return res.is_positive_definite();
330}
331
332RCP<const Basic> DenseMatrix::det() const
333{
334 return det_bareis(*this);
335}
336
337void DenseMatrix::inv(MatrixBase &result) const
338{
339 if (is_a<DenseMatrix>(result)) {
340 DenseMatrix &r = down_cast<DenseMatrix &>(result);
341 inverse_pivoted_LU(*this, r);
342 }
343}
344
345void DenseMatrix::add_matrix(const MatrixBase &other, MatrixBase &result) const
346{
347 SYMENGINE_ASSERT(row_ == result.nrows() and col_ == result.ncols());
348
349 if (is_a<DenseMatrix>(other) and is_a<DenseMatrix>(result)) {
350 const DenseMatrix &o = down_cast<const DenseMatrix &>(other);
351 DenseMatrix &r = down_cast<DenseMatrix &>(result);
352 add_dense_dense(*this, o, r);
353 }
354}
355
356void DenseMatrix::mul_matrix(const MatrixBase &other, MatrixBase &result) const
357{
358 SYMENGINE_ASSERT(row_ == result.nrows()
359 and other.ncols() == result.ncols());
360
361 if (is_a<DenseMatrix>(other) and is_a<DenseMatrix>(result)) {
362 const DenseMatrix &o = down_cast<const DenseMatrix &>(other);
363 DenseMatrix &r = down_cast<DenseMatrix &>(result);
364 mul_dense_dense(*this, o, r);
365 }
366}
367
368void DenseMatrix::elementwise_mul_matrix(const MatrixBase &other,
369 MatrixBase &result) const
370{
371 SYMENGINE_ASSERT(row_ == result.nrows() and col_ == result.ncols()
372 and row_ == other.nrows() and col_ == other.ncols());
373
374 if (is_a<DenseMatrix>(other) and is_a<DenseMatrix>(result)) {
375 const DenseMatrix &o = down_cast<const DenseMatrix &>(other);
376 DenseMatrix &r = down_cast<DenseMatrix &>(result);
377 elementwise_mul_dense_dense(*this, o, r);
378 }
379}
380
381// Add a scalar
382void DenseMatrix::add_scalar(const RCP<const Basic> &k,
383 MatrixBase &result) const
384{
385 if (is_a<DenseMatrix>(result)) {
386 DenseMatrix &r = down_cast<DenseMatrix &>(result);
387 add_dense_scalar(*this, k, r);
388 }
389}
390
391// Multiply by a scalar
392void DenseMatrix::mul_scalar(const RCP<const Basic> &k,
393 MatrixBase &result) const
394{
395 if (is_a<DenseMatrix>(result)) {
396 DenseMatrix &r = down_cast<DenseMatrix &>(result);
397 mul_dense_scalar(*this, k, r);
398 }
399}
400
401// Matrix conjugate
402void DenseMatrix::conjugate(MatrixBase &result) const
403{
404 if (is_a<DenseMatrix>(result)) {
405 DenseMatrix &r = down_cast<DenseMatrix &>(result);
406 conjugate_dense(*this, r);
407 }
408}
409
410// Matrix transpose
411void DenseMatrix::transpose(MatrixBase &result) const
412{
413 if (is_a<DenseMatrix>(result)) {
414 DenseMatrix &r = down_cast<DenseMatrix &>(result);
415 transpose_dense(*this, r);
416 }
417}
418
419// Matrix conjugate transpose
420void DenseMatrix::conjugate_transpose(MatrixBase &result) const
421{
422 if (is_a<DenseMatrix>(result)) {
423 DenseMatrix &r = down_cast<DenseMatrix &>(result);
424 conjugate_transpose_dense(*this, r);
425 }
426}
427
428// Extract out a submatrix
429void DenseMatrix::submatrix(MatrixBase &result, unsigned row_start,
430 unsigned col_start, unsigned row_end,
431 unsigned col_end, unsigned row_step,
432 unsigned col_step) const
433{
434 if (is_a<DenseMatrix>(result)) {
435 DenseMatrix &r = down_cast<DenseMatrix &>(result);
436 submatrix_dense(*this, r, row_start, col_start, row_end, col_end,
437 row_step, col_step);
438 }
439}
440
441// LU factorization
442void DenseMatrix::LU(MatrixBase &L, MatrixBase &U) const
443{
444 if (is_a<DenseMatrix>(L) and is_a<DenseMatrix>(U)) {
445 DenseMatrix &L_ = down_cast<DenseMatrix &>(L);
446 DenseMatrix &U_ = down_cast<DenseMatrix &>(U);
447 SymEngine::LU(*this, L_, U_);
448 }
449}
450
451// LDL factorization
452void DenseMatrix::LDL(MatrixBase &L, MatrixBase &D) const
453{
454 if (is_a<DenseMatrix>(L) and is_a<DenseMatrix>(D)) {
455 DenseMatrix &L_ = down_cast<DenseMatrix &>(L);
456 DenseMatrix &D_ = down_cast<DenseMatrix &>(D);
457 SymEngine::LDL(*this, L_, D_);
458 }
459}
460
461// Solve Ax = b using LU factorization
462void DenseMatrix::LU_solve(const MatrixBase &b, MatrixBase &x) const
463{
464 if (is_a<DenseMatrix>(b) and is_a<DenseMatrix>(x)) {
465 const DenseMatrix &b_ = down_cast<const DenseMatrix &>(b);
466 DenseMatrix &x_ = down_cast<DenseMatrix &>(x);
467 SymEngine::LU_solve(*this, b_, x_);
468 }
469}
470
471// Fraction free LU factorization
472void DenseMatrix::FFLU(MatrixBase &LU) const
473{
474 if (is_a<DenseMatrix>(LU)) {
475 DenseMatrix &LU_ = down_cast<DenseMatrix &>(LU);
476 fraction_free_LU(*this, LU_);
477 }
478}
479
480// Fraction free LDU factorization
481void DenseMatrix::FFLDU(MatrixBase &L, MatrixBase &D, MatrixBase &U) const
482{
483 if (is_a<DenseMatrix>(L) and is_a<DenseMatrix>(D)
484 and is_a<DenseMatrix>(U)) {
485 DenseMatrix &L_ = down_cast<DenseMatrix &>(L);
486 DenseMatrix &D_ = down_cast<DenseMatrix &>(D);
487 DenseMatrix &U_ = down_cast<DenseMatrix &>(U);
488 fraction_free_LDU(*this, L_, D_, U_);
489 }
490}
491
492// QR factorization
493void DenseMatrix::QR(MatrixBase &Q, MatrixBase &R) const
494{
495 if (is_a<DenseMatrix>(Q) and is_a<DenseMatrix>(R)) {
496 DenseMatrix &Q_ = down_cast<DenseMatrix &>(Q);
497 DenseMatrix &R_ = down_cast<DenseMatrix &>(R);
498 SymEngine::QR(*this, Q_, R_);
499 }
500}
501
502// Cholesky decomposition
503void DenseMatrix::cholesky(MatrixBase &L) const
504{
505 if (is_a<DenseMatrix>(L)) {
506 DenseMatrix &L_ = down_cast<DenseMatrix &>(L);
507 SymEngine::cholesky(*this, L_);
508 }
509}
510
511// ---------------------------- Jacobian -------------------------------------//
512
513void jacobian(const DenseMatrix &A, const DenseMatrix &x, DenseMatrix &result,
514 bool diff_cache)
515{
516 SYMENGINE_ASSERT(A.col_ == 1);
517 SYMENGINE_ASSERT(x.col_ == 1);
518 SYMENGINE_ASSERT(A.row_ == result.nrows() and x.row_ == result.ncols());
519 bool error = false;
520#pragma omp parallel for
521 for (unsigned i = 0; i < result.row_; i++) {
522 for (unsigned j = 0; j < result.col_; j++) {
523 if (is_a<Symbol>(*(x.m_[j]))) {
524 const RCP<const Symbol> x_
525 = rcp_static_cast<const Symbol>(x.m_[j]);
526 result.m_[i * result.col_ + j] = A.m_[i]->diff(x_, diff_cache);
527 } else {
528 error = true;
529 break;
530 }
531 }
532 }
533 if (error) {
534 throw SymEngineException(
535 "'x' must contain Symbols only. "
536 "Use sjacobian for SymPy style differentiation");
537 }
538}
539
540void sjacobian(const DenseMatrix &A, const DenseMatrix &x, DenseMatrix &result,
541 bool diff_cache)
542{
543 SYMENGINE_ASSERT(A.col_ == 1);
544 SYMENGINE_ASSERT(x.col_ == 1);
545 SYMENGINE_ASSERT(A.row_ == result.nrows() and x.row_ == result.ncols());
546#pragma omp parallel for
547 for (unsigned i = 0; i < result.row_; i++) {
548 for (unsigned j = 0; j < result.col_; j++) {
549 if (is_a<Symbol>(*(x.m_[j]))) {
550 const RCP<const Symbol> x_
551 = rcp_static_cast<const Symbol>(x.m_[j]);
552 result.m_[i * result.col_ + j] = A.m_[i]->diff(x_, diff_cache);
553 } else {
554 // TODO: Use a dummy symbol
555 const RCP<const Symbol> x_ = symbol("x_");
556 result.m_[i * result.col_ + j] = ssubs(
557 ssubs(A.m_[i], {{x.m_[j], x_}})->diff(x_, diff_cache),
558 {{x_, x.m_[j]}});
559 }
560 }
561 }
562}
563
564// ---------------------------- Diff -------------------------------------//
565
566void diff(const DenseMatrix &A, const RCP<const Symbol> &x, DenseMatrix &result,
567 bool diff_cache)
568{
569 SYMENGINE_ASSERT(A.row_ == result.nrows() and A.col_ == result.ncols());
570#pragma omp parallel for
571 for (unsigned i = 0; i < result.row_; i++) {
572 for (unsigned j = 0; j < result.col_; j++) {
573 result.m_[i * result.col_ + j]
574 = A.m_[i * result.col_ + j]->diff(x, diff_cache);
575 }
576 }
577}
578
579void sdiff(const DenseMatrix &A, const RCP<const Basic> &x, DenseMatrix &result,
580 bool diff_cache)
581{
582 SYMENGINE_ASSERT(A.row_ == result.nrows() and A.col_ == result.ncols());
583#pragma omp parallel for
584 for (unsigned i = 0; i < result.row_; i++) {
585 for (unsigned j = 0; j < result.col_; j++) {
586 if (is_a<Symbol>(*x)) {
587 const RCP<const Symbol> x_ = rcp_static_cast<const Symbol>(x);
588 result.m_[i * result.col_ + j]
589 = A.m_[i * result.col_ + j]->diff(x_, diff_cache);
590 } else {
591 // TODO: Use a dummy symbol
592 const RCP<const Symbol> x_ = symbol("_x");
593 result.m_[i * result.col_ + j]
594 = ssubs(ssubs(A.m_[i * result.col_ + j], {{x, x_}})
595 ->diff(x_, diff_cache),
596 {{x_, x}});
597 }
598 }
599 }
600}
601
602// ----------------------------- Matrix Conjugate ----------------------------//
603void conjugate_dense(const DenseMatrix &A, DenseMatrix &B)
604{
605 SYMENGINE_ASSERT(B.col_ == A.col_ and B.row_ == A.row_);
606
607 for (unsigned i = 0; i < A.row_; i++)
608 for (unsigned j = 0; j < A.col_; j++)
609 B.m_[i * B.col_ + j] = conjugate(A.m_[i * A.col_ + j]);
610}
611
612// ----------------------------- Matrix Transpose ----------------------------//
613void transpose_dense(const DenseMatrix &A, DenseMatrix &B)
614{
615 SYMENGINE_ASSERT(B.row_ == A.col_ and B.col_ == A.row_);
616
617 for (unsigned i = 0; i < A.row_; i++)
618 for (unsigned j = 0; j < A.col_; j++)
619 B.m_[j * B.col_ + i] = A.m_[i * A.col_ + j];
620}
621
622// ----------------------------- Matrix Conjugate Transpose -----------------//
623void conjugate_transpose_dense(const DenseMatrix &A, DenseMatrix &B)
624{
625 SYMENGINE_ASSERT(B.row_ == A.col_ and B.col_ == A.row_);
626
627 for (unsigned i = 0; i < A.row_; i++)
628 for (unsigned j = 0; j < A.col_; j++)
629 B.m_[j * B.col_ + i] = conjugate(A.m_[i * A.col_ + j]);
630}
631
632// ------------------------------- Submatrix ---------------------------------//
633void submatrix_dense(const DenseMatrix &A, DenseMatrix &B, unsigned row_start,
634 unsigned col_start, unsigned row_end, unsigned col_end,
635 unsigned row_step, unsigned col_step)
636{
637 SYMENGINE_ASSERT(row_end >= row_start and col_end >= col_start);
638 SYMENGINE_ASSERT(row_end < A.row_);
639 SYMENGINE_ASSERT(col_end < A.col_);
640 SYMENGINE_ASSERT(B.row_ == row_end - row_start + 1
641 and B.col_ == col_end - col_start + 1);
642
643 unsigned row = B.row_, col = B.col_;
644
645 for (unsigned i = 0; i < row; i += row_step)
646 for (unsigned j = 0; j < col; j += col_step)
647 B.m_[i * col + j] = A.m_[(row_start + i) * A.col_ + col_start + j];
648}
649
650// ------------------------------- Matrix Addition ---------------------------//
651void add_dense_dense(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &C)
652{
653 SYMENGINE_ASSERT(A.row_ == B.row_ and A.col_ == B.col_ and A.row_ == C.row_
654 and A.col_ == C.col_);
655
656 unsigned row = A.row_, col = A.col_;
657
658 for (unsigned i = 0; i < row; i++) {
659 for (unsigned j = 0; j < col; j++) {
660 C.m_[i * col + j] = add(A.m_[i * col + j], B.m_[i * col + j]);
661 }
662 }
663}
664
665void add_dense_scalar(const DenseMatrix &A, const RCP<const Basic> &k,
666 DenseMatrix &B)
667{
668 SYMENGINE_ASSERT(A.row_ == B.row_ and A.col_ == B.col_);
669
670 unsigned row = A.row_, col = A.col_;
671
672 for (unsigned i = 0; i < row; i++) {
673 for (unsigned j = 0; j < col; j++) {
674 B.m_[i * col + j] = add(A.m_[i * col + j], k);
675 }
676 }
677}
678
679// ------------------------------- Matrix Multiplication ---------------------//
680void mul_dense_dense(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &C)
681{
682 SYMENGINE_ASSERT(A.col_ == B.row_ and C.row_ == A.row_
683 and C.col_ == B.col_);
684
685 unsigned row = A.row_, col = B.col_;
686
687 if (&A != &C and &B != &C) {
688 for (unsigned r = 0; r < row; r++) {
689 for (unsigned c = 0; c < col; c++) {
690 C.m_[r * col + c] = zero; // Integer Zero
691 for (unsigned k = 0; k < A.col_; k++)
692 C.m_[r * col + c]
693 = add(C.m_[r * col + c],
694 mul(A.m_[r * A.col_ + k], B.m_[k * col + c]));
695 }
696 }
697 } else {
698 DenseMatrix tmp = DenseMatrix(A.row_, B.col_);
699 mul_dense_dense(A, B, tmp);
700 C = tmp;
701 }
702}
703
704void elementwise_mul_dense_dense(const DenseMatrix &A, const DenseMatrix &B,
705 DenseMatrix &C)
706{
707 SYMENGINE_ASSERT(A.row_ == B.row_ and A.col_ == B.col_ and A.row_ == C.row_
708 and A.col_ == C.col_);
709
710 unsigned row = A.row_, col = A.col_;
711
712 for (unsigned i = 0; i < row; i++) {
713 for (unsigned j = 0; j < col; j++) {
714 C.m_[i * col + j] = mul(A.m_[i * col + j], B.m_[i * col + j]);
715 }
716 }
717}
718
719void mul_dense_scalar(const DenseMatrix &A, const RCP<const Basic> &k,
720 DenseMatrix &B)
721{
722 SYMENGINE_ASSERT(A.col_ == B.col_ and A.row_ == B.row_);
723
724 unsigned row = A.row_, col = A.col_;
725
726 for (unsigned i = 0; i < row; i++) {
727 for (unsigned j = 0; j < col; j++) {
728 B.m_[i * col + j] = mul(A.m_[i * col + j], k);
729 }
730 }
731}
732
733// ---------------------------- Joining Operations ---------------------------//
734void DenseMatrix::row_join(const DenseMatrix &B)
735{
736 this->col_insert(B, col_);
737}
738
739void DenseMatrix::col_join(const DenseMatrix &B)
740{
741 this->row_insert(B, row_);
742}
743
744void DenseMatrix::row_insert(const DenseMatrix &B, unsigned pos)
745{
746 SYMENGINE_ASSERT(col_ == B.col_ and pos <= row_)
747
748 unsigned row = row_, col = col_;
749 this->resize(row_ + B.row_, col_);
750
751 for (unsigned i = row; i-- > pos;) {
752 for (unsigned j = col; j-- > 0;) {
753 this->m_[(i + B.row_) * col + j] = this->m_[i * col + j];
754 }
755 }
756
757 for (unsigned i = 0; i < B.row_; i++) {
758 for (unsigned j = 0; j < col; j++) {
759 this->m_[(i + pos) * col + j] = B.m_[i * col + j];
760 }
761 }
762}
763
764void DenseMatrix::col_insert(const DenseMatrix &B, unsigned pos)
765{
766 SYMENGINE_ASSERT(row_ == B.row_ and pos <= col_)
767
768 unsigned row = row_, col = col_;
769 this->resize(row_, col_ + B.col_);
770
771 for (unsigned i = row; i-- > 0;) {
772 for (unsigned j = col; j-- > 0;) {
773 if (j >= pos) {
774 this->m_[i * (col + B.col_) + j + B.col_]
775 = this->m_[i * col + j];
776 } else {
777 this->m_[i * (col + B.col_) + j] = this->m_[i * col + j];
778 }
779 }
780 }
781
782 for (unsigned i = 0; i < row; i++) {
783 for (unsigned j = 0; j < B.col_; j++) {
784 this->m_[i * (col + B.col_) + j + pos] = B.m_[i * B.col_ + j];
785 }
786 }
787}
788
789void DenseMatrix::row_del(unsigned k)
790{
791 SYMENGINE_ASSERT(k < row_)
792
793 if (row_ == 1)
794 this->resize(0, 0);
795 else {
796 for (unsigned i = k; i < row_ - 1; i++) {
797 row_exchange_dense(*this, i, i + 1);
798 }
799 this->resize(row_ - 1, col_);
800 }
801}
802
803void DenseMatrix::col_del(unsigned k)
804{
805 SYMENGINE_ASSERT(k < col_)
806
807 if (col_ == 1)
808 this->resize(0, 0);
809 else {
810 unsigned row = row_, col = col_, m = 0;
811
812 for (unsigned i = 0; i < row; i++) {
813 for (unsigned j = 0; j < col; j++) {
814 if (j != k) {
815 this->m_[m] = this->m_[i * col + j];
816 m++;
817 }
818 }
819 }
820 this->resize(row_, col_ - 1);
821 }
822}
823
824// -------------------------------- Row Operations ---------------------------//
825void row_exchange_dense(DenseMatrix &A, unsigned i, unsigned j)
826{
827 SYMENGINE_ASSERT(i != j and i < A.row_ and j < A.row_);
828
829 unsigned col = A.col_;
830
831 for (unsigned k = 0; k < A.col_; k++)
832 std::swap(A.m_[i * col + k], A.m_[j * col + k]);
833}
834
835void row_mul_scalar_dense(DenseMatrix &A, unsigned i, RCP<const Basic> &c)
836{
837 SYMENGINE_ASSERT(i < A.row_);
838
839 unsigned col = A.col_;
840
841 for (unsigned j = 0; j < A.col_; j++)
842 A.m_[i * col + j] = mul(c, A.m_[i * col + j]);
843}
844
845void row_add_row_dense(DenseMatrix &A, unsigned i, unsigned j,
846 RCP<const Basic> &c)
847{
848 SYMENGINE_ASSERT(i != j and i < A.row_ and j < A.row_);
849
850 unsigned col = A.col_;
851
852 for (unsigned k = 0; k < A.col_; k++)
853 A.m_[i * col + k] = add(A.m_[i * col + k], mul(c, A.m_[j * col + k]));
854}
855
856void permuteFwd(DenseMatrix &A, permutelist &pl)
857{
858 for (auto &p : pl) {
859 row_exchange_dense(A, p.first, p.second);
860 }
861}
862
863// ----------------------------- Column Operations ---------------------------//
864void column_exchange_dense(DenseMatrix &A, unsigned i, unsigned j)
865{
866 SYMENGINE_ASSERT(i != j and i < A.col_ and j < A.col_);
867
868 unsigned col = A.col_;
869
870 for (unsigned k = 0; k < A.row_; k++)
871 std::swap(A.m_[k * col + i], A.m_[k * col + j]);
872}
873
874// ------------------------------ Gaussian Elimination -----------------------//
875void pivoted_gaussian_elimination(const DenseMatrix &A, DenseMatrix &B,
876 permutelist &pl)
877{
878 SYMENGINE_ASSERT(A.row_ == B.row_ and A.col_ == B.col_);
879
880 unsigned row = A.row_, col = A.col_;
881 unsigned index = 0, i, j, k;
882 B.m_ = A.m_;
883
884 RCP<const Basic> scale;
885
886 for (i = 0; i < col - 1; i++) {
887 if (index == row)
888 break;
889
890 k = pivot(B, index, i);
891 if (k == row)
892 continue;
893 if (k != index) {
894 row_exchange_dense(B, k, index);
895 pl.push_back({k, index});
896 }
897
898 scale = div(one, B.m_[index * col + i]);
899 row_mul_scalar_dense(B, index, scale);
900
901 for (j = i + 1; j < row; j++) {
902 for (k = i + 1; k < col; k++)
903 B.m_[j * col + k]
904 = sub(B.m_[j * col + k],
905 mul(B.m_[j * col + i], B.m_[i * col + k]));
906 B.m_[j * col + i] = zero;
907 }
908
909 index++;
910 }
911}
912
913// Algorithm 1, page 12, Nakos, G. C., Turner, P. R., Williams, R. M. (1997).
914// Fraction-free algorithms for linear and polynomial equations.
915// ACM SIGSAM Bulletin, 31(3), 11–19. doi:10.1145/271130.271133.
916void fraction_free_gaussian_elimination(const DenseMatrix &A, DenseMatrix &B)
917{
918 SYMENGINE_ASSERT(A.row_ == B.row_ and A.col_ == B.col_);
919
920 unsigned col = A.col_;
921 B.m_ = A.m_;
922
923 for (unsigned i = 0; i < col - 1; i++)
924 for (unsigned j = i + 1; j < A.row_; j++) {
925 for (unsigned k = i + 1; k < col; k++) {
926 B.m_[j * col + k]
927 = sub(mul(B.m_[i * col + i], B.m_[j * col + k]),
928 mul(B.m_[j * col + i], B.m_[i * col + k]));
929 if (i > 0)
930 B.m_[j * col + k]
931 = div(B.m_[j * col + k], B.m_[i * col - col + i - 1]);
932 }
933 B.m_[j * col + i] = zero;
934 }
935}
936
937// Pivoted version of `fraction_free_gaussian_elimination`
938void pivoted_fraction_free_gaussian_elimination(const DenseMatrix &A,
939 DenseMatrix &B, permutelist &pl)
940{
941 SYMENGINE_ASSERT(A.row_ == B.row_ and A.col_ == B.col_);
942
943 unsigned col = A.col_, row = A.row_;
944 unsigned index = 0, i, k, j;
945 B.m_ = A.m_;
946
947 for (i = 0; i < col - 1; i++) {
948 if (index == row)
949 break;
950
951 k = pivot(B, index, i);
952 if (k == row)
953 continue;
954 if (k != index) {
955 row_exchange_dense(B, k, index);
956 pl.push_back({k, index});
957 }
958
959 for (j = i + 1; j < row; j++) {
960 for (k = i + 1; k < col; k++) {
961 B.m_[j * col + k]
962 = sub(mul(B.m_[i * col + i], B.m_[j * col + k]),
963 mul(B.m_[j * col + i], B.m_[i * col + k]));
964 if (i > 0)
965 B.m_[j * col + k]
966 = div(B.m_[j * col + k], B.m_[i * col - col + i - 1]);
967 }
968 B.m_[j * col + i] = zero;
969 }
970
971 index++;
972 }
973}
974
975void pivoted_gauss_jordan_elimination(const DenseMatrix &A, DenseMatrix &B,
976 permutelist &pl)
977{
978 SYMENGINE_ASSERT(A.row_ == B.row_ and A.col_ == B.col_);
979
980 unsigned row = A.row_, col = A.col_;
981 unsigned index = 0, i, j, k;
982 RCP<const Basic> scale;
983 B.m_ = A.m_;
984
985 for (i = 0; i < col; i++) {
986 if (index == row)
987 break;
988
989 k = pivot(B, index, i);
990 if (k == row)
991 continue;
992 if (k != index) {
993 row_exchange_dense(B, k, index);
994 pl.push_back({k, index});
995 }
996
997 scale = div(one, B.m_[index * col + i]);
998 row_mul_scalar_dense(B, index, scale);
999
1000 for (j = 0; j < row; j++) {
1001 if (j == index)
1002 continue;
1003
1004 scale = mul(minus_one, B.m_[j * col + i]);
1005 row_add_row_dense(B, j, index, scale);
1006 }
1007
1008 index++;
1009 }
1010}
1011
1012// Algorithm 2, page 13, Nakos, G. C., Turner, P. R., Williams, R. M. (1997).
1013// Fraction-free algorithms for linear and polynomial equations.
1014// ACM SIGSAM Bulletin, 31(3), 11–19. doi:10.1145/271130.271133.
1015void fraction_free_gauss_jordan_elimination(const DenseMatrix &A,
1016 DenseMatrix &B)
1017{
1018 SYMENGINE_ASSERT(A.row_ == B.row_ and A.col_ == B.col_);
1019
1020 unsigned row = A.row_, col = A.col_;
1021 unsigned i, j, k;
1022 RCP<const Basic> d;
1023
1024 B.m_ = A.m_;
1025
1026 for (i = 0; i < col; i++) {
1027 if (i > 0)
1028 d = B.m_[i * col - col + i - 1];
1029 for (j = 0; j < row; j++)
1030 if (j != i)
1031 for (k = 0; k < col; k++) {
1032 if (k != i) {
1033 B.m_[j * col + k]
1034 = sub(mul(B.m_[i * col + i], B.m_[j * col + k]),
1035 mul(B.m_[j * col + i], B.m_[i * col + k]));
1036 if (i > 0)
1037 B.m_[j * col + k] = div(B.m_[j * col + k], d);
1038 }
1039 }
1040 for (j = 0; j < row; j++)
1041 if (j != i)
1042 B.m_[j * col + i] = zero;
1043 }
1044}
1045
1046void pivoted_fraction_free_gauss_jordan_elimination(const DenseMatrix &A,
1047 DenseMatrix &B,
1048 permutelist &pl)
1049{
1050 SYMENGINE_ASSERT(A.row_ == B.row_ and A.col_ == B.col_);
1051
1052 unsigned row = A.row_, col = A.col_;
1053 unsigned index = 0, i, k, j;
1054 RCP<const Basic> d;
1055
1056 B.m_ = A.m_;
1057
1058 for (i = 0; i < col; i++) {
1059 if (index == row)
1060 break;
1061
1062 k = pivot(B, index, i);
1063 if (k == row)
1064 continue;
1065 if (k != index) {
1066 row_exchange_dense(B, k, index);
1067 pl.push_back({k, index});
1068 }
1069
1070 for (j = 0; j < row; j++) {
1071 if (j == index)
1072 continue;
1073 for (k = 0; k < col; k++) {
1074 if (k == i)
1075 continue;
1076 B.m_[j * col + k]
1077 = sub(mul(B.m_[index * col + i], B.m_[j * col + k]),
1078 mul(B.m_[j * col + i], B.m_[index * col + k]));
1079 if (index == 0)
1080 continue;
1081 B.m_[j * col + k] = div(B.m_[j * col + k], d);
1082 }
1083 }
1084
1085 d = B.m_[index * col + i];
1086
1087 for (j = 0; j < row; j++) {
1088 if (j == index)
1089 continue;
1090 B.m_[j * col + i] = zero;
1091 }
1092 index++;
1093 }
1094}
1095
1096unsigned pivot(DenseMatrix &B, unsigned r, unsigned c)
1097{
1098 for (unsigned k = r; k < B.row_; k++) {
1099 if (!is_true(is_zero(*(B.m_[k * B.col_ + c])))) {
1100 return k;
1101 }
1102 }
1103 return B.row_;
1104}
1105
1106void reduced_row_echelon_form(const DenseMatrix &A, DenseMatrix &b,
1107 vec_uint &pivot_cols, bool normalize_last)
1108{
1109 permutelist pl;
1110 if (normalize_last) {
1111 pivoted_fraction_free_gauss_jordan_elimination(A, b, pl);
1112 } else {
1113 pivoted_gauss_jordan_elimination(A, b, pl);
1114 }
1115 unsigned row = 0;
1116 for (unsigned col = 0; col < b.col_ && row < b.row_; col++) {
1117 if (is_true(is_zero(*b.get(row, col))))
1118 continue;
1119 pivot_cols.push_back(col);
1120 if (row == 0 and normalize_last) {
1121 RCP<const Basic> m = div(one, b.get(row, col));
1122 b.mul_scalar(m, b);
1123 }
1124 row++;
1125 }
1126}
1127
1128// --------------------------- Solve Ax = b ---------------------------------//
1129// Assuming A is a diagonal square matrix
1130void diagonal_solve(const DenseMatrix &A, const DenseMatrix &b, DenseMatrix &x)
1131{
1132 SYMENGINE_ASSERT(A.row_ == A.col_);
1133 SYMENGINE_ASSERT(x.row_ == A.col_ and x.col_ == b.col_);
1134
1135 const unsigned sys = b.col_;
1136
1137 // No checks are done to see if the diagonal entries are zero
1138 for (unsigned k = 0; k < sys; k++) {
1139 for (unsigned i = 0; i < A.col_; i++) {
1140 x.m_[i * sys + k] = div(b.m_[i * sys + k], A.m_[i * A.col_ + i]);
1141 }
1142 }
1143}
1144
1145void back_substitution(const DenseMatrix &U, const DenseMatrix &b,
1146 DenseMatrix &x)
1147{
1148 SYMENGINE_ASSERT(U.row_ == U.col_);
1149 SYMENGINE_ASSERT(b.row_ == U.row_);
1150 SYMENGINE_ASSERT(x.row_ == U.col_ and x.col_ == b.col_);
1151
1152 const unsigned col = U.col_;
1153 const unsigned sys = b.col_;
1154 x.m_ = b.m_;
1155
1156 for (unsigned k = 0; k < sys; k++) {
1157 for (int i = col - 1; i >= 0; i--) {
1158 for (unsigned j = i + 1; j < col; j++)
1159 x.m_[i * sys + k]
1160 = sub(x.m_[i * sys + k],
1161 mul(U.m_[i * col + j], x.m_[j * sys + k]));
1162 x.m_[i * sys + k] = div(x.m_[i * sys + k], U.m_[i * col + i]);
1163 }
1164 }
1165}
1166
1167void forward_substitution(const DenseMatrix &A, const DenseMatrix &b,
1168 DenseMatrix &x)
1169{
1170 SYMENGINE_ASSERT(A.row_ == A.col_);
1171 SYMENGINE_ASSERT(b.row_ == A.row_);
1172 SYMENGINE_ASSERT(x.row_ == A.col_ and x.col_ == b.col_);
1173
1174 unsigned col = A.col_;
1175 const unsigned sys = b.col_;
1176 x.m_ = b.m_;
1177
1178 for (unsigned k = 0; k < b.col_; k++) {
1179 for (unsigned i = 0; i < col - 1; i++) {
1180 for (unsigned j = i + 1; j < col; j++) {
1181 x.m_[j * sys + k]
1182 = sub(mul(A.m_[i * col + i], x.m_[j * sys + k]),
1183 mul(A.m_[j * col + i], x.m_[i * sys + k]));
1184 if (i > 0)
1185 x.m_[j * sys + k]
1186 = div(x.m_[j * sys + k], A.m_[(i - 1) * col + i - 1]);
1187 }
1188 }
1189 }
1190}
1191
1192void fraction_free_gaussian_elimination_solve(const DenseMatrix &A,
1193 const DenseMatrix &b,
1194 DenseMatrix &x)
1195{
1196 SYMENGINE_ASSERT(A.row_ == A.col_);
1197 SYMENGINE_ASSERT(b.row_ == A.row_ and x.row_ == A.row_);
1198 SYMENGINE_ASSERT(x.col_ == b.col_);
1199
1200 int i, j, k, col = A.col_, bcol = b.col_;
1201 DenseMatrix A_ = DenseMatrix(A.row_, A.col_, A.m_);
1202 DenseMatrix b_ = DenseMatrix(b.row_, b.col_, b.m_);
1203
1204 for (i = 0; i < col - 1; i++)
1205 for (j = i + 1; j < col; j++) {
1206 for (k = 0; k < bcol; k++) {
1207 b_.m_[j * bcol + k]
1208 = sub(mul(A_.m_[i * col + i], b_.m_[j * bcol + k]),
1209 mul(A_.m_[j * col + i], b_.m_[i * bcol + k]));
1210 if (i > 0)
1211 b_.m_[j * bcol + k] = div(b_.m_[j * bcol + k],
1212 A_.m_[i * col - col + i - 1]);
1213 }
1214
1215 for (k = i + 1; k < col; k++) {
1216 A_.m_[j * col + k]
1217 = sub(mul(A_.m_[i * col + i], A_.m_[j * col + k]),
1218 mul(A_.m_[j * col + i], A_.m_[i * col + k]));
1219 if (i > 0)
1220 A_.m_[j * col + k]
1221 = div(A_.m_[j * col + k], A_.m_[i * col - col + i - 1]);
1222 }
1223 A_.m_[j * col + i] = zero;
1224 }
1225
1226 for (i = 0; i < col * bcol; i++)
1227 x.m_[i] = zero; // Integer zero;
1228
1229 for (k = 0; k < bcol; k++) {
1230 for (i = col - 1; i >= 0; i--) {
1231 for (j = i + 1; j < col; j++)
1232 b_.m_[i * bcol + k]
1233 = sub(b_.m_[i * bcol + k],
1234 mul(A_.m_[i * col + j], x.m_[j * bcol + k]));
1235 x.m_[i * bcol + k] = div(b_.m_[i * bcol + k], A_.m_[i * col + i]);
1236 }
1237 }
1238}
1239
1240void fraction_free_gauss_jordan_solve(const DenseMatrix &A,
1241 const DenseMatrix &b, DenseMatrix &x,
1242 bool pivot)
1243{
1244 SYMENGINE_ASSERT(A.row_ == A.col_);
1245 SYMENGINE_ASSERT(b.row_ == A.row_ and x.row_ == A.row_);
1246 SYMENGINE_ASSERT(x.col_ == b.col_);
1247
1248 unsigned i, j, k, p, col = A.col_, bcol = b.col_;
1249 RCP<const Basic> d;
1250 DenseMatrix A_ = DenseMatrix(A.row_, A.col_, A.m_);
1251 DenseMatrix b_ = DenseMatrix(b.row_, b.col_, b.m_);
1252
1253 for (i = 0; i < col; i++) {
1254 if (i > 0)
1255 d = A_.m_[i * col - col + i - 1];
1256 if (pivot) {
1257 p = i;
1258 while (p < col and eq(*A_.m_[p * col + i], *zero)) {
1259 p++;
1260 }
1261 SYMENGINE_ASSERT(p != col);
1262 if (p != i) {
1263 // pivot A
1264 for (k = i; k < col; k++) {
1265 std::swap(A_.m_[p * col + k], A_.m_[i * col + k]);
1266 }
1267 for (k = 0; k < bcol; k++) {
1268 std::swap(b_.m_[p * bcol + k], b_.m_[i * bcol + k]);
1269 }
1270 }
1271 }
1272 for (j = 0; j < col; j++)
1273 if (j != i) {
1274 for (k = 0; k < bcol; k++) {
1275 b_.m_[j * bcol + k]
1276 = sub(mul(A_.m_[i * col + i], b_.m_[j * bcol + k]),
1277 mul(A_.m_[j * col + i], b_.m_[i * bcol + k]));
1278 if (i > 0)
1279 b_.m_[j * bcol + k] = div(b_.m_[j * bcol + k], d);
1280 }
1281
1282 for (k = 0; k < col; k++) {
1283 if (k != i) {
1284 A_.m_[j * col + k]
1285 = sub(mul(A_.m_[i * col + i], A_.m_[j * col + k]),
1286 mul(A_.m_[j * col + i], A_.m_[i * col + k]));
1287 if (i > 0)
1288 A_.m_[j * col + k] = div(A_.m_[j * col + k], d);
1289 }
1290 }
1291 }
1292
1293 for (j = 0; j < col; j++)
1294 if (j != i)
1295 A_.m_[j * col + i] = zero;
1296 }
1297
1298 // No checks are done to see if the diagonal entries are zero
1299 for (k = 0; k < bcol; k++)
1300 for (i = 0; i < col; i++)
1301 x.m_[i * bcol + k] = div(b_.m_[i * bcol + k], A_.m_[i * col + i]);
1302}
1303
1304void fraction_free_LU_solve(const DenseMatrix &A, const DenseMatrix &b,
1305 DenseMatrix &x)
1306{
1307 DenseMatrix LU = DenseMatrix(A.nrows(), A.ncols());
1308 DenseMatrix x_ = DenseMatrix(b.nrows(), b.ncols());
1309
1310 fraction_free_LU(A, LU);
1311 forward_substitution(LU, b, x_);
1312 back_substitution(LU, x_, x);
1313}
1314
1315void LU_solve(const DenseMatrix &A, const DenseMatrix &b, DenseMatrix &x)
1316{
1317 DenseMatrix L = DenseMatrix(A.nrows(), A.ncols());
1318 DenseMatrix U = DenseMatrix(A.nrows(), A.ncols());
1319 DenseMatrix x_ = DenseMatrix(b.nrows(), b.ncols());
1320
1321 LU(A, L, U);
1322 forward_substitution(L, b, x_);
1323 back_substitution(U, x_, x);
1324}
1325
1326void pivoted_LU_solve(const DenseMatrix &A, const DenseMatrix &b,
1327 DenseMatrix &x)
1328{
1329 DenseMatrix L = DenseMatrix(A.nrows(), A.ncols());
1330 DenseMatrix U = DenseMatrix(A.nrows(), A.ncols());
1331 DenseMatrix x_ = DenseMatrix(b);
1332 permutelist pl;
1333
1334 pivoted_LU(A, L, U, pl);
1335 permuteFwd(x_, pl);
1336 forward_substitution(L, x_, x_);
1337 back_substitution(U, x_, x);
1338}
1339
1340void LDL_solve(const DenseMatrix &A, const DenseMatrix &b, DenseMatrix &x)
1341{
1342 DenseMatrix L = DenseMatrix(A.nrows(), A.ncols());
1343 DenseMatrix D = DenseMatrix(A.nrows(), A.ncols());
1344 DenseMatrix x_ = DenseMatrix(b.nrows(), b.ncols());
1345
1346 if (not is_symmetric_dense(A))
1347 throw SymEngineException("Matrix must be symmetric");
1348
1349 LDL(A, L, D);
1350 forward_substitution(L, b, x);
1351 diagonal_solve(D, x, x_);
1352 transpose_dense(L, D);
1353 back_substitution(D, x_, x);
1354}
1355
1356// --------------------------- Matrix Decomposition --------------------------//
1357
1358// Algorithm 3, page 14, Nakos, G. C., Turner, P. R., Williams, R. M. (1997).
1359// Fraction-free algorithms for linear and polynomial equations.
1360// ACM SIGSAM Bulletin, 31(3), 11–19. doi:10.1145/271130.271133.
1361// This algorithms is not a true factorization of the matrix A(i.e. A != LU))
1362// but can be used to solve linear systems by applying forward and backward
1363// substitutions respectively.
1364void fraction_free_LU(const DenseMatrix &A, DenseMatrix &LU)
1365{
1366 SYMENGINE_ASSERT(A.row_ == A.col_ and LU.row_ == LU.col_
1367 and A.row_ == LU.row_);
1368
1369 unsigned n = A.row_;
1370 unsigned i, j, k;
1371
1372 LU.m_ = A.m_;
1373
1374 for (i = 0; i < n - 1; i++)
1375 for (j = i + 1; j < n; j++)
1376 for (k = i + 1; k < n; k++) {
1377 LU.m_[j * n + k] = sub(mul(LU.m_[i * n + i], LU.m_[j * n + k]),
1378 mul(LU.m_[j * n + i], LU.m_[i * n + k]));
1379 if (i)
1380 LU.m_[j * n + k]
1381 = div(LU.m_[j * n + k], LU.m_[i * n - n + i - 1]);
1382 }
1383}
1384
1385// SymPy LUDecomposition algorithm, in
1386// sympy.matrices.matrices.Matrix.LUdecomposition
1387// with no pivoting
1388void LU(const DenseMatrix &A, DenseMatrix &L, DenseMatrix &U)
1389{
1390 SYMENGINE_ASSERT(A.row_ == A.col_ and L.row_ == L.col_
1391 and U.row_ == U.col_);
1392 SYMENGINE_ASSERT(A.row_ == L.row_ and A.row_ == U.row_);
1393
1394 unsigned n = A.row_;
1395 unsigned i, j, k;
1396 RCP<const Basic> scale;
1397
1398 U.m_ = A.m_;
1399
1400 for (j = 0; j < n; j++) {
1401 for (i = 0; i < j; i++)
1402 for (k = 0; k < i; k++)
1403 U.m_[i * n + j] = sub(U.m_[i * n + j],
1404 mul(U.m_[i * n + k], U.m_[k * n + j]));
1405
1406 for (i = j; i < n; i++) {
1407 for (k = 0; k < j; k++)
1408 U.m_[i * n + j] = sub(U.m_[i * n + j],
1409 mul(U.m_[i * n + k], U.m_[k * n + j]));
1410 }
1411
1412 scale = div(one, U.m_[j * n + j]);
1413
1414 for (i = j + 1; i < n; i++)
1415 U.m_[i * n + j] = mul(U.m_[i * n + j], scale);
1416 }
1417
1418 for (i = 0; i < n; i++) {
1419 for (j = 0; j < i; j++) {
1420 L.m_[i * n + j] = U.m_[i * n + j];
1421 U.m_[i * n + j] = zero; // Integer zero
1422 }
1423 L.m_[i * n + i] = one; // Integer one
1424 for (j = i + 1; j < n; j++)
1425 L.m_[i * n + j] = zero; // Integer zero
1426 }
1427}
1428
1429// SymPy LUDecomposition_Simple algorithm, in
1430// sympy.matrices.matrices.Matrix.LUdecomposition_Simple with pivoting.
1431// P must be an initialized matrix and will be permuted.
1432void pivoted_LU(const DenseMatrix &A, DenseMatrix &LU, permutelist &pl)
1433{
1434 SYMENGINE_ASSERT(A.row_ == A.col_ and LU.row_ == LU.col_);
1435 SYMENGINE_ASSERT(A.row_ == LU.row_);
1436
1437 unsigned n = A.row_;
1438 unsigned i, j, k;
1439 RCP<const Basic> scale;
1440 int pivot;
1441
1442 LU.m_ = A.m_;
1443
1444 for (j = 0; j < n; j++) {
1445 for (i = 0; i < j; i++)
1446 for (k = 0; k < i; k++)
1447 LU.m_[i * n + j] = sub(LU.m_[i * n + j],
1448 mul(LU.m_[i * n + k], LU.m_[k * n + j]));
1449 pivot = -1;
1450 for (i = j; i < n; i++) {
1451 for (k = 0; k < j; k++) {
1452 LU.m_[i * n + j] = sub(LU.m_[i * n + j],
1453 mul(LU.m_[i * n + k], LU.m_[k * n + j]));
1454 }
1455 if (pivot == -1 and !is_true(is_zero(*LU.m_[i * n + j])))
1456 pivot = i;
1457 }
1458 if (pivot == -1)
1459 throw SymEngineException("Matrix is rank deficient");
1460 if (pivot - j != 0) { // row must be swapped
1461 row_exchange_dense(LU, pivot, j);
1462 pl.push_back({pivot, j});
1463 }
1464 scale = div(one, LU.m_[j * n + j]);
1465
1466 for (i = j + 1; i < n; i++)
1467 LU.m_[i * n + j] = mul(LU.m_[i * n + j], scale);
1468 }
1469}
1470
1471// SymPy LUDecomposition algorithm, in
1472// sympy.matrices.matrices.Matrix.LUdecomposition_Simple with pivoting.
1473// P must be an initialized matrix and will be permuted.
1474void pivoted_LU(const DenseMatrix &A, DenseMatrix &L, DenseMatrix &U,
1475 permutelist &pl)
1476{
1477 SYMENGINE_ASSERT(A.row_ == A.col_ and L.row_ == L.col_
1478 and U.row_ == U.col_);
1479 SYMENGINE_ASSERT(A.row_ == L.row_ and A.row_ == U.row_);
1480
1481 pivoted_LU(A, U, pl);
1482 unsigned n = A.col_;
1483 for (unsigned i = 0; i < n; i++) {
1484 for (unsigned j = 0; j < i; j++) {
1485 L.m_[i * n + j] = U.m_[i * n + j];
1486 U.m_[i * n + j] = zero; // Integer zero
1487 }
1488 L.m_[i * n + i] = one; // Integer one
1489 for (unsigned j = i + 1; j < n; j++)
1490 L.m_[i * n + j] = zero; // Integer zero
1491 }
1492}
1493
1494// SymPy's fraction free LU decomposition, without pivoting
1495// sympy.matrices.matrices.MatrixBase.LUdecompositionFF
1496// W. Zhou & D.J. Jeffrey, "Fraction-free matrix factors: new forms for LU and
1497// QR factors".
1498// Frontiers in Computer Science in China, Vol 2, no. 1, pp. 67-80, 2008.
1499void fraction_free_LDU(const DenseMatrix &A, DenseMatrix &L, DenseMatrix &D,
1500 DenseMatrix &U)
1501{
1502 SYMENGINE_ASSERT(A.row_ == L.row_ and A.row_ == U.row_);
1503 SYMENGINE_ASSERT(A.col_ == L.col_ and A.col_ == U.col_);
1504
1505 unsigned row = A.row_, col = A.col_;
1506 unsigned i, j, k;
1507 RCP<const Basic> old = integer(1);
1508
1509 U.m_ = A.m_;
1510
1511 // Initialize L
1512 for (i = 0; i < row; i++)
1513 for (j = 0; j < row; j++)
1514 if (i != j)
1515 L.m_[i * col + j] = zero;
1516 else
1517 L.m_[i * col + i] = one;
1518
1519 // Initialize D
1520 for (i = 0; i < row * row; i++)
1521 D.m_[i] = zero; // Integer zero
1522
1523 for (k = 0; k < row - 1; k++) {
1524 L.m_[k * col + k] = U.m_[k * col + k];
1525 D.m_[k * col + k] = mul(old, U.m_[k * col + k]);
1526
1527 for (i = k + 1; i < row; i++) {
1528 L.m_[i * col + k] = U.m_[i * col + k];
1529 for (j = k + 1; j < col; j++)
1530 U.m_[i * col + j]
1531 = div(sub(mul(U.m_[k * col + k], U.m_[i * col + j]),
1532 mul(U.m_[k * col + j], U.m_[i * col + k])),
1533 old);
1534 U.m_[i * col + k] = zero; // Integer zero
1535 }
1536
1537 old = U.m_[k * col + k];
1538 }
1539
1540 D.m_[row * col - col + row - 1] = old;
1541}
1542
1543// SymPy's QRecomposition in sympy.matrices.matrices.MatrixBase.QRdecomposition
1544// Rank check is not performed
1545void QR(const DenseMatrix &A, DenseMatrix &Q, DenseMatrix &R)
1546{
1547 unsigned row = A.row_;
1548 unsigned col = A.col_;
1549
1550 SYMENGINE_ASSERT(Q.row_ == row and Q.col_ == col and R.row_ == col
1551 and R.col_ == col);
1552
1553 unsigned i, j, k;
1554 RCP<const Basic> t;
1556
1557 // Initialize Q
1558 for (i = 0; i < row * col; i++)
1559 Q.m_[i] = zero;
1560
1561 // Initialize R
1562 for (i = 0; i < col * col; i++)
1563 R.m_[i] = zero;
1564
1565 for (j = 0; j < col; j++) {
1566 // Use submatrix for this
1567 for (k = 0; k < row; k++)
1568 tmp[k] = A.m_[k * col + j];
1569
1570 for (i = 0; i < j; i++) {
1571 t = zero;
1572 for (k = 0; k < row; k++)
1573 t = add(t, mul(A.m_[k * col + j], Q.m_[k * col + i]));
1574 for (k = 0; k < row; k++)
1575 tmp[k] = expand(sub(tmp[k], mul(Q.m_[k * col + i], t)));
1576 }
1577
1578 // calculate norm
1579 t = zero;
1580 for (k = 0; k < row; k++)
1581 t = add(t, pow(tmp[k], integer(2)));
1582
1583 t = pow(t, div(one, integer(2)));
1584
1585 R.m_[j * col + j] = t;
1586 for (k = 0; k < row; k++)
1587 Q.m_[k * col + j] = div(tmp[k], t);
1588
1589 for (i = 0; i < j; i++) {
1590 t = zero;
1591 for (k = 0; k < row; k++)
1592 t = add(t, mul(Q.m_[k * col + i], A.m_[k * col + j]));
1593 R.m_[i * col + j] = t;
1594 }
1595 }
1596}
1597
1598// SymPy's LDL decomposition, Assuming A is a symmetric, square, positive
1599// definite non singular matrix
1600void LDL(const DenseMatrix &A, DenseMatrix &L, DenseMatrix &D)
1601{
1602 SYMENGINE_ASSERT(A.row_ == A.col_);
1603 SYMENGINE_ASSERT(L.row_ == A.row_ and L.col_ == A.row_);
1604 SYMENGINE_ASSERT(D.row_ == A.row_ and D.col_ == A.row_);
1605
1606 unsigned col = A.col_;
1607 unsigned i, k, j;
1608 RCP<const Basic> sum;
1609 RCP<const Basic> i2 = integer(2);
1610
1611 // Initialize D
1612 for (i = 0; i < col; i++)
1613 for (j = 0; j < col; j++)
1614 D.m_[i * col + j] = zero; // Integer zero
1615
1616 // Initialize L
1617 for (i = 0; i < col; i++)
1618 for (j = 0; j < col; j++)
1619 L.m_[i * col + j] = (i != j) ? zero : one;
1620
1621 for (i = 0; i < col; i++) {
1622 for (j = 0; j < i; j++) {
1623 sum = zero;
1624 for (k = 0; k < j; k++)
1625 sum = add(sum, mul(mul(L.m_[i * col + k], L.m_[j * col + k]),
1626 D.m_[k * col + k]));
1627 L.m_[i * col + j]
1628 = mul(div(one, D.m_[j * col + j]), sub(A.m_[i * col + j], sum));
1629 }
1630 sum = zero;
1631 for (k = 0; k < i; k++)
1632 sum = add(sum, mul(pow(L.m_[i * col + k], i2), D.m_[k * col + k]));
1633 D.m_[i * col + i] = sub(A.m_[i * col + i], sum);
1634 }
1635}
1636
1637// SymPy's cholesky decomposition
1638void cholesky(const DenseMatrix &A, DenseMatrix &L)
1639{
1640 SYMENGINE_ASSERT(A.row_ == A.col_);
1641 SYMENGINE_ASSERT(L.row_ == A.row_ and L.col_ == A.row_);
1642
1643 unsigned col = A.col_;
1644 unsigned i, j, k;
1645 RCP<const Basic> sum;
1646 RCP<const Basic> i2 = integer(2);
1647 RCP<const Basic> half = div(one, i2);
1648
1649 // Initialize L
1650 for (i = 0; i < col; i++)
1651 for (j = 0; j < col; j++)
1652 L.m_[i * col + j] = zero;
1653
1654 for (i = 0; i < col; i++) {
1655 for (j = 0; j < i; j++) {
1656 sum = zero;
1657 for (k = 0; k < j; k++)
1658 sum = add(sum, mul(L.m_[i * col + k], L.m_[j * col + k]));
1659
1660 L.m_[i * col + j]
1661 = mul(div(one, L.m_[j * col + j]), sub(A.m_[i * col + j], sum));
1662 }
1663 sum = zero;
1664 for (k = 0; k < i; k++)
1665 sum = add(sum, pow(L.m_[i * col + k], i2));
1666 L.m_[i * col + i] = pow(sub(A.m_[i * col + i], sum), half);
1667 }
1668}
1669
1670// Matrix Queries
1671bool is_symmetric_dense(const DenseMatrix &A)
1672{
1673 if (A.col_ != A.row_)
1674 return false;
1675
1676 unsigned col = A.col_;
1677 bool sym = true;
1678
1679 for (unsigned i = 0; i < col; i++)
1680 for (unsigned j = i + 1; j < col; j++)
1681 if (not eq(*(A.m_[j * col + i]), *(A.m_[i * col + j]))) {
1682 sym = false;
1683 break;
1684 }
1685
1686 return sym;
1687}
1688
1689// ----------------------------- Determinant ---------------------------------//
1690RCP<const Basic> det_bareis(const DenseMatrix &A)
1691{
1692 SYMENGINE_ASSERT(A.row_ == A.col_);
1693
1694 unsigned n = A.row_;
1695
1696 if (n == 1) {
1697 return A.m_[0];
1698 } else if (n == 2) {
1699 // If A = [[a, b], [c, d]] then det(A) = ad - bc
1700 return sub(mul(A.m_[0], A.m_[3]), mul(A.m_[1], A.m_[2]));
1701 } else if (n == 3) {
1702 // if A = [[a, b, c], [d, e, f], [g, h, i]] then
1703 // det(A) = (aei + bfg + cdh) - (ceg + bdi + afh)
1704 return sub(add(add(mul(mul(A.m_[0], A.m_[4]), A.m_[8]),
1705 mul(mul(A.m_[1], A.m_[5]), A.m_[6])),
1706 mul(mul(A.m_[2], A.m_[3]), A.m_[7])),
1707 add(add(mul(mul(A.m_[2], A.m_[4]), A.m_[6]),
1708 mul(mul(A.m_[1], A.m_[3]), A.m_[8])),
1709 mul(mul(A.m_[0], A.m_[5]), A.m_[7])));
1710 } else {
1711
1712 if (A.is_lower() or A.is_upper()) {
1713 RCP<const Basic> det = A.m_[0];
1714 for (unsigned i = 1; i < n; ++i) {
1715 det = mul(det, A.m_[i * n + i]);
1716 }
1717 return det;
1718 }
1719
1720 DenseMatrix B = DenseMatrix(n, n, A.m_);
1721 unsigned i, sign = 1;
1722 RCP<const Basic> d;
1723
1724 for (unsigned k = 0; k < n - 1; k++) {
1725 if (is_true(is_zero(*B.m_[k * n + k]))) {
1726 for (i = k + 1; i < n; i++)
1727 if (!is_true(is_zero(*B.m_[i * n + k]))) {
1728 row_exchange_dense(B, i, k);
1729 sign *= -1;
1730 break;
1731 }
1732 if (i == n)
1733 return zero;
1734 }
1735
1736 for (i = k + 1; i < n; i++) {
1737 for (unsigned j = k + 1; j < n; j++) {
1738 d = sub(mul(B.m_[k * n + k], B.m_[i * n + j]),
1739 mul(B.m_[i * n + k], B.m_[k * n + j]));
1740 if (k > 0)
1741 d = div(d, B.m_[(k - 1) * n + k - 1]);
1742 B.m_[i * n + j] = d;
1743 }
1744 }
1745 }
1746
1747 return (sign == 1) ? B.m_[n * n - 1] : mul(minus_one, B.m_[n * n - 1]);
1748 }
1749}
1750
1751// Returns the coefficients of characterisitc polynomials of leading principal
1752// minors of Matrix A as elements in `polys`. Principal leading minor of kth
1753// order is the submatrix of A obtained by deleting last n-k rows and columns
1754// from A. Here `n` is the dimension of the square matrix A.
1755void berkowitz(const DenseMatrix &A, std::vector<DenseMatrix> &polys)
1756{
1757 SYMENGINE_ASSERT(A.row_ == A.col_);
1758
1759 unsigned col = A.col_;
1760 unsigned i, k, l, m;
1761
1763 std::vector<DenseMatrix> transforms;
1765
1766 for (unsigned n = col; n > 1; n--) {
1767 items.clear();
1768 k = n - 1;
1769 DenseMatrix T = DenseMatrix(n + 1, n);
1770 DenseMatrix C = DenseMatrix(k, 1);
1771
1772 // Initialize T and C
1773 for (i = 0; i < n * (n + 1); i++)
1774 T.m_[i] = zero;
1775 for (i = 0; i < k; i++)
1776 C.m_[i] = A.m_[i * col + k];
1777 items.push_back(C);
1778
1779 for (i = 0; i < n - 2; i++) {
1780 DenseMatrix B = DenseMatrix(k, 1);
1781 for (l = 0; l < k; l++) {
1782 B.m_[l] = zero;
1783 for (m = 0; m < k; m++)
1784 B.m_[l]
1785 = add(B.m_[l], mul(A.m_[l * col + m], items[i].m_[m]));
1786 }
1787 items.push_back(B);
1788 }
1789
1790 items_.clear();
1791 for (i = 0; i < n - 1; i++) {
1792 RCP<const Basic> element = zero;
1793 for (l = 0; l < k; l++)
1794 element = add(element, mul(A.m_[k * col + l], items[i].m_[l]));
1795 items_.push_back(mul(minus_one, element));
1796 }
1797 items_.insert(items_.begin(), mul(minus_one, A.m_[k * col + k]));
1798 items_.insert(items_.begin(), one);
1799
1800 for (i = 0; i < n; i++) {
1801 for (l = 0; l < n - i + 1; l++)
1802 T.m_[(i + l) * n + i] = items_[l];
1803 }
1804
1805 transforms.push_back(T);
1806 }
1807
1808 polys.push_back(DenseMatrix(2, 1, {one, mul(A.m_[0], minus_one)}));
1809
1810 for (i = 0; i < col - 1; i++) {
1811 unsigned t_row = transforms[col - 2 - i].nrows();
1812 unsigned t_col = transforms[col - 2 - i].ncols();
1813 DenseMatrix B = DenseMatrix(t_row, 1);
1814
1815 for (l = 0; l < t_row; l++) {
1816 B.m_[l] = zero;
1817 for (m = 0; m < t_col; m++) {
1818 B.m_[l] = add(B.m_[l],
1819 mul(transforms[col - 2 - i].m_[l * t_col + m],
1820 polys[i].m_[m]));
1821 B.m_[l] = expand(B.m_[l]);
1822 }
1823 }
1824 polys.push_back(B);
1825 }
1826}
1827
1828RCP<const Basic> det_berkowitz(const DenseMatrix &A)
1829{
1831
1832 berkowitz(A, polys);
1833 DenseMatrix poly = polys[polys.size() - 1];
1834
1835 if (polys.size() % 2 == 1)
1836 return mul(minus_one, poly.get(poly.nrows() - 1, 0));
1837
1838 return poly.get(poly.nrows() - 1, 0);
1839}
1840
1841void char_poly(const DenseMatrix &A, DenseMatrix &B)
1842{
1843 SYMENGINE_ASSERT(B.ncols() == 1 and B.nrows() == A.nrows() + 1);
1844 SYMENGINE_ASSERT(A.nrows() == A.ncols());
1845
1847
1848 berkowitz(A, polys);
1849 B = polys[polys.size() - 1];
1850}
1851
1852void inverse_fraction_free_LU(const DenseMatrix &A, DenseMatrix &B)
1853{
1854 SYMENGINE_ASSERT(A.row_ == A.col_ and B.row_ == B.col_
1855 and B.row_ == A.row_);
1856
1857 unsigned n = A.row_, i;
1858 DenseMatrix LU = DenseMatrix(n, n);
1859 DenseMatrix e = DenseMatrix(n, 1);
1860 DenseMatrix x = DenseMatrix(n, 1);
1861 DenseMatrix x_ = DenseMatrix(n, 1);
1862
1863 // Initialize matrices
1864 for (i = 0; i < n * n; i++) {
1865 LU.m_[i] = zero;
1866 B.m_[i] = zero;
1867 }
1868 for (i = 0; i < n; i++) {
1869 e.m_[i] = zero;
1870 x.m_[i] = zero;
1871 x_.m_[i] = zero;
1872 }
1873
1874 fraction_free_LU(A, LU);
1875
1876 // We solve AX_{i} = e_{i} for i = 1, 2, .. n and combine the row vectors
1877 // X_{1}, X_{2}, ... X_{n} to form the inverse of A. Here, e_{i}'s are the
1878 // elements of the standard basis.
1879 for (unsigned j = 0; j < n; j++) {
1880 e.m_[j] = one;
1881
1882 forward_substitution(LU, e, x_);
1883 back_substitution(LU, x_, x);
1884
1885 for (i = 0; i < n; i++)
1886 B.m_[i * n + j] = x.m_[i];
1887
1888 e.m_[j] = zero;
1889 }
1890}
1891
1892void inverse_LU(const DenseMatrix &A, DenseMatrix &B)
1893{
1894 SYMENGINE_ASSERT(A.row_ == A.col_ and B.row_ == B.col_
1895 and B.row_ == A.row_);
1896 DenseMatrix e = DenseMatrix(A.row_, A.col_);
1897 eye(e);
1898 LU_solve(A, e, B);
1899}
1900
1901void inverse_pivoted_LU(const DenseMatrix &A, DenseMatrix &B)
1902{
1903 SYMENGINE_ASSERT(A.row_ == A.col_ and B.row_ == B.col_
1904 and B.row_ == A.row_);
1905 DenseMatrix e = DenseMatrix(A.row_, A.col_);
1906 eye(e);
1907 pivoted_LU_solve(A, e, B);
1908}
1909
1910void inverse_gauss_jordan(const DenseMatrix &A, DenseMatrix &B)
1911{
1912 SYMENGINE_ASSERT(A.row_ == A.col_ and B.row_ == B.col_
1913 and B.row_ == A.row_);
1914
1915 unsigned n = A.row_;
1916 DenseMatrix e = DenseMatrix(n, n);
1917
1918 // Initialize matrices
1919 for (unsigned i = 0; i < n; i++)
1920 for (unsigned j = 0; j < n; j++) {
1921 if (i != j) {
1922 e.m_[i * n + j] = zero;
1923 } else {
1924 e.m_[i * n + i] = one;
1925 }
1926 B.m_[i * n + j] = zero;
1927 }
1928
1929 fraction_free_gauss_jordan_solve(A, e, B);
1930}
1931
1932// ----------------------- Vector-specific Methods --------------------------//
1933
1934void dot(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &C)
1935{
1936 if (A.col_ == B.row_) {
1937 if (B.col_ != 1) {
1938 DenseMatrix tmp1 = DenseMatrix(A.col_, A.row_);
1939 A.transpose(tmp1);
1940 DenseMatrix tmp2 = DenseMatrix(B.col_, B.row_);
1941 B.transpose(tmp2);
1942 C.resize(tmp1.row_, tmp2.col_);
1943 mul_dense_dense(tmp1, tmp2, C);
1944 } else {
1945 C.resize(A.row_, B.col_);
1946 mul_dense_dense(A, B, C);
1947 }
1948 C.resize(1, C.row_ * C.col_);
1949 } else if (A.col_ == B.col_) {
1950 DenseMatrix tmp2 = DenseMatrix(B.col_, B.row_);
1951 B.transpose(tmp2);
1952 dot(A, tmp2, C);
1953 } else if (A.row_ == B.row_) {
1954 DenseMatrix tmp1 = DenseMatrix(A.col_, A.row_);
1955 A.transpose(tmp1);
1956 dot(tmp1, B, C);
1957 } else {
1958 throw SymEngineException("Dimensions incorrect for dot product");
1959 }
1960}
1961
1962void cross(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &C)
1963{
1964 SYMENGINE_ASSERT((A.row_ * A.col_ == 3 and B.row_ * B.col_ == 3)
1965 and (A.row_ == C.row_ and A.col_ == C.col_));
1966
1967 C.m_[0] = sub(mul(A.m_[1], B.m_[2]), mul(A.m_[2], B.m_[1]));
1968 C.m_[1] = sub(mul(A.m_[2], B.m_[0]), mul(A.m_[0], B.m_[2]));
1969 C.m_[2] = sub(mul(A.m_[0], B.m_[1]), mul(A.m_[1], B.m_[0]));
1970}
1971
1972RCP<const Set> eigen_values(const DenseMatrix &A)
1973{
1974 unsigned n = A.nrows();
1975 if (A.is_lower() or A.is_upper()) {
1976 RCP<const Set> eigenvals = emptyset();
1977 set_basic x;
1978 for (unsigned i = 0; i < n; ++i) {
1979 x.insert(A.get(i, i));
1980 }
1981 eigenvals = finiteset(x);
1982 return eigenvals;
1983 }
1984
1985 DenseMatrix B = DenseMatrix(A.nrows() + 1, 1);
1986 char_poly(A, B);
1987 map_int_Expr coeffs;
1988 auto nr = A.nrows();
1989 for (unsigned i = 0; i <= nr; i++)
1990 insert(coeffs, nr - i, B.get(i, 0));
1991 auto lambda = symbol("lambda");
1992 return solve_poly(uexpr_poly(lambda, coeffs), lambda);
1993}
1994
1995// ------------------------- NumPy-like functions ----------------------------//
1996
1997// Mimic `eye` function in NumPy
1998void eye(DenseMatrix &A, int k)
1999{
2000 if ((k >= 0 and (unsigned) k >= A.col_) or k + A.row_ <= 0) {
2001 zeros(A);
2002 }
2003
2004 vec_basic v = vec_basic(k > 0 ? A.col_ - k : A.row_ + k, one);
2005
2006 diag(A, v, k);
2007}
2008
2009// Create diagonal matrices directly
2010void diag(DenseMatrix &A, vec_basic &v, int k)
2011{
2012 SYMENGINE_ASSERT(v.size() > 0);
2013
2014 unsigned k_ = std::abs(k);
2015
2016 if (k >= 0) {
2017 for (unsigned i = 0; i < A.row_; i++) {
2018 for (unsigned j = 0; j < A.col_; j++) {
2019 if (j != (unsigned)k) {
2020 A.m_[i * A.col_ + j] = zero;
2021 } else {
2022 A.m_[i * A.col_ + j] = v[k - k_];
2023 }
2024 }
2025 k++;
2026 }
2027 } else {
2028 k = -k;
2029
2030 for (unsigned j = 0; j < A.col_; j++) {
2031 for (unsigned i = 0; i < A.row_; i++) {
2032 if (i != (unsigned)k) {
2033 A.m_[i * A.col_ + j] = zero;
2034 } else {
2035 A.m_[i * A.col_ + j] = v[k - k_];
2036 }
2037 }
2038 k++;
2039 }
2040 }
2041}
2042
2043// Create a matrix filled with ones
2044void ones(DenseMatrix &A)
2045{
2046 for (unsigned i = 0; i < A.row_ * A.col_; i++) {
2047 A.m_[i] = one;
2048 }
2049}
2050
2051// Create a matrix filled with zeros
2052void zeros(DenseMatrix &A)
2053{
2054 for (unsigned i = 0; i < A.row_ * A.col_; i++) {
2055 A.m_[i] = zero;
2056 }
2057}
2058
2059} // namespace SymEngine
Classes and functions relating to the binary operation of addition.
T begin(T... args)
T clear(T... args)
T div(T... args)
T get(T... args)
T insert(T... args)
Main namespace for SymEngine package.
Definition: add.cpp:19
RCP< const EmptySet > emptyset()
Definition: sets.h:562
bool is_number_and_zero(const Basic &b)
Definition: number.h:139
RCP< const Set > finiteset(const set_basic &container)
Definition: sets.h:574
bool eq(const Basic &a, const Basic &b)
Checks equality for a and b
Definition: basic-inl.h:21
RCP< const Basic > sign(const RCP< const Basic > &arg)
Canonicalize Sign.
Definition: functions.cpp:527
RCP< const Basic > abs(const RCP< const Basic > &arg)
Canonicalize Abs:
Definition: functions.cpp:3492
RCP< const Basic > sub(const RCP< const Basic > &a, const RCP< const Basic > &b)
Substracts b from a.
Definition: add.cpp:495
RCP< const Basic > mul(const RCP< const Basic > &a, const RCP< const Basic > &b)
Multiplication.
Definition: mul.cpp:352
void insert(T1 &m, const T2 &first, const T3 &second)
Definition: dict.h:83
RCP< const Symbol > symbol(const std::string &name)
inline version to return Symbol
Definition: symbol.h:82
RCP< const Basic > add(const RCP< const Basic > &a, const RCP< const Basic > &b)
Adds two objects (safely).
Definition: add.cpp:425
tribool is_zero(const Basic &b, const Assumptions *assumptions=nullptr)
Check if a number is zero.
RCP< const Basic > expand(const RCP< const Basic > &self, bool deep=true)
Expands self
Definition: expand.cpp:369
std::enable_if< std::is_integral< T >::value, RCP< constInteger > >::type integer(T i)
Definition: integer.h:200
RCP< const Basic > ssubs(const RCP< const Basic > &x, const map_basic_basic &subs_dict, bool cache=true)
SymPy compatible subs.
Definition: subs.h:550
RCP< const Basic > conjugate(const RCP< const Basic > &arg)
Canonicalize Conjugate.
Definition: functions.cpp:149
T pow(T... args)
T push_back(T... args)
T size(T... args)
T swap(T... args)