devtools/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2026-04-10 11:34:33 +08:00
Code Issues Packages Projects Releases Wiki Activity

Add the implementation of the Incomplete LU preconditioner with dual threshold (ILUT)

Browse Source 
Modify the BiCGSTAB function to check the residual norm of the initial guess
This commit is contained in:
Desire NUENTSA
2012-02-10 10:59:39 +01:00
parent 9ed6a267a3
commit a815d962da
3 changed files with 414 additions and 9 deletions
Download Patch File Download Diff File Expand all files Collapse all files

16
Eigen/src/IterativeLinearSolvers/BiCGSTAB.h
View File

@@ -53,6 +53,7 @@ void bicgstab(const MatrixType& mat, const Rhs& rhs, Dest& x,
int n = mat.cols();
VectorType r = rhs - mat * x;
VectorType r0 = r;
RealScalar r0_sqnorm = r0.squaredNorm();
Scalar rho = 1;
Scalar alpha = 1;
@@ -67,15 +68,17 @@ void bicgstab(const MatrixType& mat, const Rhs& rhs, Dest& x,
RealScalar tol2 = tol*tol;
int i = 0;
do
while ( r.squaredNorm()/r0_sqnorm > tol2 && i<maxIters )
{
Scalar rho_old = rho;
rho = r0.dot(r);
eigen_assert((rho != Scalar(0)) && "BiCGSTAB BROKE DOWN !!!");
Scalar beta = (rho/rho_old) * (alpha / w);
p = r + beta * (p - w * v);
y = precond.solve(p);
v.noalias() = mat * y;
alpha = rho / r0.dot(v);
@@ -84,17 +87,12 @@ void bicgstab(const MatrixType& mat, const Rhs& rhs, Dest& x,
z = precond.solve(s);
t.noalias() = mat * z;
kt = precond.solve(t);
ks = precond.solve(s);
w = kt.dot(ks) / kt.squaredNorm();
w = t.dot(s) / t.squaredNorm();
x += alpha * y + w * z;
r = s - w * t;
++i;
} while ( r.squaredNorm()/r0_sqnorm > tol2 && i<maxIters );
}
tol_error = sqrt(r.squaredNorm()/r0_sqnorm);
//tol_error = sqrt(abs(absNew / absInit));
iters = i;
}
@@ -233,7 +231,7 @@ public:
template<typename Rhs,typename Dest>
void _solve(const Rhs& b, Dest& x) const
{
x.setOnes();
x.setZero();
_solveWithGuess(b,x);
}