mirror of
https://github.com/UpsilonNumworks/Upsilon.git
synced 2026-03-18 21:30:38 +01:00
[poincare/matrix_inverse] Handle memory error in MatrixInverse
This commit is contained in:
@@ -163,7 +163,8 @@ Calculation::EqualSign Calculation::exactAndApproximateDisplayedOutputsAreEqual(
|
||||
* checkpoint: if there was not enough memory on the pool to compute the equal
|
||||
* sign, just return EqualSign::Approximation.
|
||||
* We can safely use an exception checkpoint here because we are sure of not
|
||||
* modifying any pre-existing node in the pool. */
|
||||
* modifying any pre-existing node in the pool. We are sure there cannot be a
|
||||
* Store in the exactOutput. */
|
||||
Poincare::ExceptionCheckpoint ecp;
|
||||
if (ExceptionRun(ecp)) {
|
||||
constexpr int bufferSize = Constant::MaxSerializedExpressionSize;
|
||||
|
||||
@@ -80,7 +80,7 @@ public:
|
||||
Matrix createTranspose() const;
|
||||
/* createInverse can be called on any matrix, reduced or not, approximated or
|
||||
* not. */
|
||||
Expression createInverse(ExpressionNode::ReductionContext reductionContext) const;
|
||||
Expression createInverse(ExpressionNode::ReductionContext reductionContext, bool * couldComputeInverse) const;
|
||||
#if MATRIX_EXACT_REDUCING
|
||||
Expression determinant() const;
|
||||
#endif
|
||||
|
||||
@@ -1,5 +1,6 @@
|
||||
#include <poincare/matrix.h>
|
||||
#include <poincare/division.h>
|
||||
#include <poincare/exception_checkpoint.h>
|
||||
#include <poincare/matrix_complex.h>
|
||||
#include <poincare/matrix_layout.h>
|
||||
#include <poincare/multiplication_explicite.h>
|
||||
@@ -301,38 +302,52 @@ Matrix Matrix::createTranspose() const {
|
||||
return matrix;
|
||||
}
|
||||
|
||||
Expression Matrix::createInverse(ExpressionNode::ReductionContext reductionContext) const {
|
||||
Expression Matrix::createInverse(ExpressionNode::ReductionContext reductionContext, bool * couldComputeInverse) const {
|
||||
int dim = numberOfRows();
|
||||
if (dim != numberOfColumns()) {
|
||||
return Undefined::Builder();
|
||||
}
|
||||
/* Create the matrix inv = (A|I) with A is the input matrix and I the dim
|
||||
* identity matrix */
|
||||
Matrix matrixAI = Matrix::Builder();
|
||||
for (int i = 0; i < dim; i++) {
|
||||
for (int j = 0; j < dim; j++) {
|
||||
matrixAI.addChildAtIndexInPlace(const_cast<Matrix *>(this)->matrixChild(i, j).clone(), i*2*dim+j, i*2*dim+j);
|
||||
/* If the matrix is too big, the inversion might not be computed exactly
|
||||
* because of a pool allocation error, but we might still be able to compute
|
||||
* it approximately. We thus encapsulate the inversion creation in an
|
||||
* exception checkpoint.
|
||||
* We can safely use an exception checkpoint here because we are sure of not
|
||||
* modifying any pre-existing node in the pool. We are sure there is no Store
|
||||
* in the matrix. */
|
||||
Poincare::ExceptionCheckpoint ecp;
|
||||
if (ExceptionRun(ecp)) {
|
||||
/* Create the matrix inv = (A|I) with A is the input matrix and I the dim
|
||||
* identity matrix */
|
||||
Matrix matrixAI = Matrix::Builder();
|
||||
for (int i = 0; i < dim; i++) {
|
||||
for (int j = 0; j < dim; j++) {
|
||||
matrixAI.addChildAtIndexInPlace(const_cast<Matrix *>(this)->matrixChild(i, j).clone(), i*2*dim+j, i*2*dim+j);
|
||||
}
|
||||
for (int j = dim; j < 2*dim; j++) {
|
||||
matrixAI.addChildAtIndexInPlace(j-dim == i ? Rational::Builder(1) : Rational::Builder(0), i*2*dim+j, i*2*dim+j);
|
||||
}
|
||||
}
|
||||
for (int j = dim; j < 2*dim; j++) {
|
||||
matrixAI.addChildAtIndexInPlace(j-dim == i ? Rational::Builder(1) : Rational::Builder(0), i*2*dim+j, i*2*dim+j);
|
||||
matrixAI.setDimensions(dim, 2*dim);
|
||||
matrixAI = matrixAI.rowCanonize(reductionContext);
|
||||
// Check inversibility
|
||||
for (int i = 0; i < dim; i++) {
|
||||
if (!matrixAI.matrixChild(i, i).isRationalOne()) {
|
||||
return Undefined::Builder();
|
||||
}
|
||||
}
|
||||
Matrix inverse = Matrix::Builder();
|
||||
for (int i = 0; i < dim; i++) {
|
||||
for (int j = 0; j < dim; j++) {
|
||||
inverse.addChildAtIndexInPlace(matrixAI.matrixChild(i, j+dim), i*dim+j, i*dim+j); // We can steal matrixAI's children
|
||||
}
|
||||
}
|
||||
inverse.setDimensions(dim, dim);
|
||||
*couldComputeInverse = true;
|
||||
return inverse;
|
||||
} else {
|
||||
*couldComputeInverse = false;
|
||||
return Expression();
|
||||
}
|
||||
matrixAI.setDimensions(dim, 2*dim);
|
||||
matrixAI = matrixAI.rowCanonize(reductionContext);
|
||||
// Check inversibility
|
||||
for (int i = 0; i < dim; i++) {
|
||||
if (!matrixAI.matrixChild(i, i).isRationalOne()) {
|
||||
return Undefined::Builder();
|
||||
}
|
||||
}
|
||||
Matrix inverse = Matrix::Builder();
|
||||
for (int i = 0; i < dim; i++) {
|
||||
for (int j = 0; j < dim; j++) {
|
||||
inverse.addChildAtIndexInPlace(matrixAI.matrixChild(i, j+dim), i*dim+j, i*dim+j); // We can steal matrixAI's children
|
||||
}
|
||||
}
|
||||
inverse.setDimensions(dim, dim);
|
||||
return inverse;
|
||||
}
|
||||
|
||||
template int Matrix::ArrayInverse<float>(float *, int, int);
|
||||
|
||||
@@ -57,9 +57,15 @@ Expression MatrixInverse::shallowReduce(ExpressionNode::ReductionContext reducti
|
||||
}
|
||||
/* Power(matrix, -n) creates a matrixInverse, so the simplification must be
|
||||
* done here and not in power. */
|
||||
Expression result = matrixChild.createInverse(reductionContext);
|
||||
replaceWithInPlace(result);
|
||||
return result.shallowReduce(reductionContext);
|
||||
bool couldComputeInverse = false;
|
||||
Expression result = matrixChild.createInverse(reductionContext, &couldComputeInverse);
|
||||
if (couldComputeInverse) {
|
||||
replaceWithInPlace(result);
|
||||
return result.shallowReduce(reductionContext);
|
||||
}
|
||||
// The matrix could not be inverted exactly
|
||||
// TODO Poincare error?
|
||||
return *this;
|
||||
}
|
||||
Expression result = Power::Builder(c, Rational::Builder(-1));
|
||||
replaceWithInPlace(result);
|
||||
|
||||
Reference in New Issue
Block a user