diff options
| -rw-r--r-- | Eigen/src/IterativeLinearSolvers/BiCGSTAB.h | 5 | ||||
| -rw-r--r-- | Eigen/src/IterativeLinearSolvers/ConjugateGradient.h | 24 | ||||
| -rw-r--r-- | test/sparse_solver.h | 56 | ||||
| -rw-r--r-- | unsupported/Eigen/src/IterativeSolvers/GMRES.h | 2 |
4 files changed, 57 insertions, 30 deletions
diff --git a/Eigen/src/IterativeLinearSolvers/BiCGSTAB.h b/Eigen/src/IterativeLinearSolvers/BiCGSTAB.h index 5a822e0ea..fbefb696f 100644 --- a/Eigen/src/IterativeLinearSolvers/BiCGSTAB.h +++ b/Eigen/src/IterativeLinearSolvers/BiCGSTAB.h @@ -44,6 +44,11 @@ bool bicgstab(const MatrixType& mat, const Rhs& rhs, Dest& x, VectorType r0 = r; RealScalar r0_sqnorm = rhs.squaredNorm(); + if(r0_sqnorm == 0) + { + x.setZero(); + return true; + } Scalar rho = 1; Scalar alpha = 1; Scalar w = 1; diff --git a/Eigen/src/IterativeLinearSolvers/ConjugateGradient.h b/Eigen/src/IterativeLinearSolvers/ConjugateGradient.h index eadf711e5..00b5647c6 100644 --- a/Eigen/src/IterativeLinearSolvers/ConjugateGradient.h +++ b/Eigen/src/IterativeLinearSolvers/ConjugateGradient.h @@ -41,21 +41,29 @@ void conjugate_gradient(const MatrixType& mat, const Rhs& rhs, Dest& x, int n = mat.cols(); VectorType residual = rhs - mat * x; //initial residual - VectorType p(n); - - p = precond.solve(residual); //initial search direction - VectorType z(n), tmp(n); - 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) + if(rhsNorm2 == 0) { x.setZero(); + iters = 0; + tol_error = 0; return; } - RealScalar residualNorm2 = 0; 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 + + VectorType z(n), tmp(n); + RealScalar absNew = internal::real(residual.dot(p)); // the square of the absolute value of r scaled by invM int i = 0; while(i < maxIters) { diff --git a/test/sparse_solver.h b/test/sparse_solver.h index 5a1be67e7..6b3c48274 100644 --- a/test/sparse_solver.h +++ b/test/sparse_solver.h @@ -176,25 +176,32 @@ template<typename Solver> void check_sparse_spd_solving(Solver& solver) // generate the problem Mat A, halfA; DenseMatrix dA; - int size = generate_sparse_spd_problem(solver, A, halfA, dA); - - // generate the right hand sides - int rhsCols = internal::random<int>(1,16); - double density = (std::max)(8./(size*rhsCols), 0.1); - SpMat B(size,rhsCols); - DenseVector b = DenseVector::Random(size); - DenseMatrix dB(size,rhsCols); - initSparse<Scalar>(density, dB, B, ForceNonZeroDiag); - for (int i = 0; i < g_repeat; i++) { + int size = generate_sparse_spd_problem(solver, A, halfA, dA); + + // generate the right hand sides + int rhsCols = internal::random<int>(1,16); + double density = (std::max)(8./(size*rhsCols), 0.1); + SpMat B(size,rhsCols); + DenseVector b = DenseVector::Random(size); + DenseMatrix dB(size,rhsCols); + initSparse<Scalar>(density, dB, B, ForceNonZeroDiag); + check_sparse_solving(solver, A, b, dA, b); check_sparse_solving(solver, halfA, b, dA, b); check_sparse_solving(solver, A, dB, dA, dB); check_sparse_solving(solver, halfA, dB, dA, dB); check_sparse_solving(solver, A, 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 #ifdef TEST_REAL_CASES if (internal::is_same<Scalar, float>::value @@ -265,21 +272,28 @@ template<typename Solver> void check_sparse_square_solving(Solver& solver) Mat A; DenseMatrix dA; - int size = generate_sparse_square_problem(solver, A, dA); - - A.makeCompressed(); - DenseVector b = DenseVector::Random(size); - DenseMatrix dB(size,rhsCols); - SpMat B(size,rhsCols); - double density = (std::max)(8./(size*rhsCols), 0.1); - initSparse<Scalar>(density, dB, B, ForceNonZeroDiag); - B.makeCompressed(); for (int i = 0; i < g_repeat; i++) { + int size = generate_sparse_square_problem(solver, A, dA); + + A.makeCompressed(); + DenseVector b = DenseVector::Random(size); + DenseMatrix dB(size,rhsCols); + SpMat B(size,rhsCols); + double density = (std::max)(8./(size*rhsCols), 0.1); + initSparse<Scalar>(density, dB, B, ForceNonZeroDiag); + B.makeCompressed(); check_sparse_solving(solver, A, b, dA, b); check_sparse_solving(solver, A, dB, 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 #ifdef TEST_REAL_CASES if (internal::is_same<Scalar, float>::value diff --git a/unsupported/Eigen/src/IterativeSolvers/GMRES.h b/unsupported/Eigen/src/IterativeSolvers/GMRES.h index 2efb6ff92..f82472159 100644 --- a/unsupported/Eigen/src/IterativeSolvers/GMRES.h +++ b/unsupported/Eigen/src/IterativeSolvers/GMRES.h @@ -349,7 +349,7 @@ public: void _solve(const Rhs& b, Dest& x) const { x = b; - if(!x.squaredNorm()) return; // Check Zero right hand side + if(x.squaredNorm() == 0) return; // Check Zero right hand side _solveWithGuess(b,x); } |