Bug567 : Fix iterative solvers to immediately return when the initial guess is the true solution and for trivial solution

This commit is contained in:
Desire NUENTSA 2013-03-20 16:15:18 +01:00
parent 22460edb49
commit da6219b19d
4 changed files with 59 additions and 32 deletions

View File

@ -44,6 +44,11 @@ bool bicgstab(const MatrixType& mat, const Rhs& rhs, Dest& x,
VectorType r0 = r; VectorType r0 = r;
RealScalar r0_sqnorm = rhs.squaredNorm(); RealScalar r0_sqnorm = rhs.squaredNorm();
if(r0_sqnorm == 0)
{
x.setZero();
return true;
}
Scalar rho = 1; Scalar rho = 1;
Scalar alpha = 1; Scalar alpha = 1;
Scalar w = 1; Scalar w = 1;

View File

@ -41,21 +41,29 @@ void conjugate_gradient(const MatrixType& mat, const Rhs& rhs, Dest& x,
int n = mat.cols(); int n = mat.cols();
VectorType residual = rhs - mat * x; //initial residual VectorType residual = rhs - mat * x; //initial residual
VectorType p(n);
RealScalar rhsNorm2 = rhs.squaredNorm();
if(rhsNorm2 == 0)
{
x.setZero();
iters = 0;
tol_error = 0;
return;
}
RealScalar threshold = tol*tol*rhsNorm2;
RealScalar residualNorm2 = residual.squaredNorm();
if (residualNorm2 < threshold)
{
iters = 0;
tol_error = sqrt(residualNorm2 / rhsNorm2);
return;
}
VectorType p(n);
p = precond.solve(residual); //initial search direction p = precond.solve(residual); //initial search direction
VectorType z(n), tmp(n); VectorType z(n), tmp(n);
RealScalar absNew = internal::real(residual.dot(p)); // the square of the absolute value of r scaled by invM RealScalar absNew = internal::real(residual.dot(p)); // the square of the absolute value of r scaled by invM
RealScalar rhsNorm2 = rhs.squaredNorm();
// Check Zero right hand side
if(!rhsNorm2)
{
x.setZero();
return;
}
RealScalar residualNorm2 = 0;
RealScalar threshold = tol*tol*rhsNorm2;
int i = 0; int i = 0;
while(i < maxIters) while(i < maxIters)
{ {

View File

@ -176,6 +176,7 @@ template<typename Solver> void check_sparse_spd_solving(Solver& solver)
// generate the problem // generate the problem
Mat A, halfA; Mat A, halfA;
DenseMatrix dA; DenseMatrix dA;
for (int i = 0; i < g_repeat; i++) {
int size = generate_sparse_spd_problem(solver, A, halfA, dA); int size = generate_sparse_spd_problem(solver, A, halfA, dA);
// generate the right hand sides // generate the right hand sides
@ -186,13 +187,19 @@ template<typename Solver> void check_sparse_spd_solving(Solver& solver)
DenseMatrix dB(size,rhsCols); DenseMatrix dB(size,rhsCols);
initSparse<Scalar>(density, dB, B, ForceNonZeroDiag); initSparse<Scalar>(density, dB, B, ForceNonZeroDiag);
for (int i = 0; i < g_repeat; i++) {
check_sparse_solving(solver, A, b, dA, b); check_sparse_solving(solver, A, b, dA, b);
check_sparse_solving(solver, halfA, b, dA, b); check_sparse_solving(solver, halfA, b, dA, b);
check_sparse_solving(solver, A, dB, dA, dB); check_sparse_solving(solver, A, dB, dA, dB);
check_sparse_solving(solver, halfA, dB, dA, dB); check_sparse_solving(solver, halfA, dB, dA, dB);
check_sparse_solving(solver, A, B, dA, dB); check_sparse_solving(solver, A, B, dA, dB);
check_sparse_solving(solver, halfA, B, dA, dB); check_sparse_solving(solver, halfA, B, dA, dB);
// check only once
if(i==0)
{
b = DenseVector::Zero(size);
check_sparse_solving(solver, A, b, dA, b);
}
} }
// First, get the folder // First, get the folder
@ -265,6 +272,7 @@ template<typename Solver> void check_sparse_square_solving(Solver& solver)
Mat A; Mat A;
DenseMatrix dA; DenseMatrix dA;
for (int i = 0; i < g_repeat; i++) {
int size = generate_sparse_square_problem(solver, A, dA); int size = generate_sparse_square_problem(solver, A, dA);
A.makeCompressed(); A.makeCompressed();
@ -274,10 +282,16 @@ template<typename Solver> void check_sparse_square_solving(Solver& solver)
double density = (std::max)(8./(size*rhsCols), 0.1); double density = (std::max)(8./(size*rhsCols), 0.1);
initSparse<Scalar>(density, dB, B, ForceNonZeroDiag); initSparse<Scalar>(density, dB, B, ForceNonZeroDiag);
B.makeCompressed(); B.makeCompressed();
for (int i = 0; i < g_repeat; i++) {
check_sparse_solving(solver, A, b, dA, b); check_sparse_solving(solver, A, b, dA, b);
check_sparse_solving(solver, A, dB, dA, dB); check_sparse_solving(solver, A, dB, dA, dB);
check_sparse_solving(solver, A, B, dA, dB); check_sparse_solving(solver, A, B, dA, dB);
// check only once
if(i==0)
{
b = DenseVector::Zero(size);
check_sparse_solving(solver, A, b, dA, b);
}
} }
// First, get the folder // First, get the folder

View File

@ -349,7 +349,7 @@ public:
void _solve(const Rhs& b, Dest& x) const void _solve(const Rhs& b, Dest& x) const
{ {
x = b; x = b;
if(!x.squaredNorm()) return; // Check Zero right hand side if(x.squaredNorm() == 0) return; // Check Zero right hand side
_solveWithGuess(b,x); _solveWithGuess(b,x);
} }