From 62670c83a0ba7cb4f45a734a4817a818a7c92bba Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Mon, 10 Jun 2013 23:40:56 +0200 Subject: Fix bug #314: move remaining math functions from internal to numext namespace --- Eigen/src/IterativeLinearSolvers/ConjugateGradient.h | 4 ++-- Eigen/src/IterativeLinearSolvers/IncompleteLUT.h | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) (limited to 'Eigen/src/IterativeLinearSolvers') diff --git a/Eigen/src/IterativeLinearSolvers/ConjugateGradient.h b/Eigen/src/IterativeLinearSolvers/ConjugateGradient.h index 00b5647c6..a74a8155e 100644 --- a/Eigen/src/IterativeLinearSolvers/ConjugateGradient.h +++ b/Eigen/src/IterativeLinearSolvers/ConjugateGradient.h @@ -63,7 +63,7 @@ void conjugate_gradient(const MatrixType& mat, const Rhs& rhs, Dest& x, 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 absNew = numext::real(residual.dot(p)); // the square of the absolute value of r scaled by invM int i = 0; while(i < maxIters) { @@ -80,7 +80,7 @@ void conjugate_gradient(const MatrixType& mat, const Rhs& rhs, Dest& x, z = precond.solve(residual); // approximately solve for "A z = residual" RealScalar absOld = absNew; - absNew = internal::real(residual.dot(z)); // update the absolute value of r + absNew = numext::real(residual.dot(z)); // update the absolute value of r RealScalar beta = absNew / absOld; // calculate the Gram-Schmidt value used to create the new search direction p = z + beta * p; // update search direction i++; diff --git a/Eigen/src/IterativeLinearSolvers/IncompleteLUT.h b/Eigen/src/IterativeLinearSolvers/IncompleteLUT.h index 17d18ef58..50a870aec 100644 --- a/Eigen/src/IterativeLinearSolvers/IncompleteLUT.h +++ b/Eigen/src/IterativeLinearSolvers/IncompleteLUT.h @@ -310,7 +310,7 @@ void IncompleteLUT::factorize(const _MatrixType& amat) jr(k) = jpos; ++sizeu; } - rownorm += internal::abs2(j_it.value()); + rownorm += numext::abs2(j_it.value()); } // 2 - detect possible zero row -- cgit v1.2.3 From 22820e950e916917cebfc8aa1633162b847b53f7 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Mon, 1 Jul 2013 11:49:23 +0200 Subject: Improve BiCGSTAB robustness: fix a divide by zero and allow to restart with a new initial residual reference. --- Eigen/src/IterativeLinearSolvers/BiCGSTAB.h | 24 ++++++++++++++++++------ 1 file changed, 18 insertions(+), 6 deletions(-) (limited to 'Eigen/src/IterativeLinearSolvers') diff --git a/Eigen/src/IterativeLinearSolvers/BiCGSTAB.h b/Eigen/src/IterativeLinearSolvers/BiCGSTAB.h index fbefb696f..048d59170 100644 --- a/Eigen/src/IterativeLinearSolvers/BiCGSTAB.h +++ b/Eigen/src/IterativeLinearSolvers/BiCGSTAB.h @@ -43,8 +43,9 @@ bool bicgstab(const MatrixType& mat, const Rhs& rhs, Dest& x, VectorType r = rhs - mat * x; VectorType r0 = r; - RealScalar r0_sqnorm = rhs.squaredNorm(); - if(r0_sqnorm == 0) + RealScalar r0_sqnorm = r0.squaredNorm(); + RealScalar rhs_sqnorm = rhs.squaredNorm(); + if(rhs_sqnorm == 0) { x.setZero(); return true; @@ -61,13 +62,22 @@ bool bicgstab(const MatrixType& mat, const Rhs& rhs, Dest& x, RealScalar tol2 = tol*tol; int i = 0; + int restarts = 0; - while ( r.squaredNorm()/r0_sqnorm > tol2 && i tol2 && iRealScalar(0)) + w = t.dot(s) / t.squaredNorm(); x += alpha * y + w * z; r = s - w * t; ++i; } - tol_error = sqrt(r.squaredNorm()/r0_sqnorm); + tol_error = sqrt(r.squaredNorm()/rhs_sqnorm); iters = i; return true; } -- cgit v1.2.3 From 1caeb814f05b9477bcf385b31cc8f4cc95142bc2 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Tue, 2 Jul 2013 08:14:10 +0200 Subject: Fix bicgstab for complexes, and avoid a duplicate computation --- Eigen/src/IterativeLinearSolvers/BiCGSTAB.h | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) (limited to 'Eigen/src/IterativeLinearSolvers') diff --git a/Eigen/src/IterativeLinearSolvers/BiCGSTAB.h b/Eigen/src/IterativeLinearSolvers/BiCGSTAB.h index 048d59170..6fc6ab852 100644 --- a/Eigen/src/IterativeLinearSolvers/BiCGSTAB.h +++ b/Eigen/src/IterativeLinearSolvers/BiCGSTAB.h @@ -74,7 +74,7 @@ bool bicgstab(const MatrixType& mat, const Rhs& rhs, Dest& x, // The new residual vector became too orthogonal to the arbitrarily choosen direction r0 // Let's restart with a new r0: r0 = r; - r0_sqnorm = rho = r0.dot(r); + rho = r0_sqnorm = r.squaredNorm(); if(restarts++ == 0) i = 0; } @@ -91,9 +91,11 @@ bool bicgstab(const MatrixType& mat, const Rhs& rhs, Dest& x, z = precond.solve(s); t.noalias() = mat * z; - w = t.squaredNorm(); - if(w>RealScalar(0)) - w = t.dot(s) / t.squaredNorm(); + RealScalar tmp = t.squaredNorm(); + if(tmp>RealScalar(0)) + w = t.dot(s) / tmp; + else + w = Scalar(0); x += alpha * y + w * z; r = s - w * t; ++i; -- cgit v1.2.3 From bd689ccc2836ac0854961c7a3dc5b1151e576a9e Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Wed, 17 Jul 2013 09:21:07 +0200 Subject: IncompleteLUT should not raise an assert in compute if factorize failed. --- Eigen/src/IterativeLinearSolvers/IncompleteLUT.h | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) (limited to 'Eigen/src/IterativeLinearSolvers') diff --git a/Eigen/src/IterativeLinearSolvers/IncompleteLUT.h b/Eigen/src/IterativeLinearSolvers/IncompleteLUT.h index 50a870aec..b55afc136 100644 --- a/Eigen/src/IterativeLinearSolvers/IncompleteLUT.h +++ b/Eigen/src/IterativeLinearSolvers/IncompleteLUT.h @@ -150,8 +150,7 @@ class IncompleteLUT : internal::noncopyable { analyzePattern(amat); factorize(amat); - eigen_assert(m_factorizationIsOk == true); - m_isInitialized = true; + m_isInitialized = m_factorizationIsOk; return *this; } -- cgit v1.2.3