26 #ifndef EIGEN_BICGSTAB_H
27 #define EIGEN_BICGSTAB_H
43 template<
typename MatrixType,
typename Rhs,
typename Dest,
typename Preconditioner>
44 bool bicgstab(
const MatrixType& mat,
const Rhs& rhs, Dest& x,
45 const Preconditioner& precond,
int& iters,
46 typename Dest::RealScalar& tol_error)
50 typedef typename Dest::RealScalar RealScalar;
51 typedef typename Dest::Scalar Scalar;
53 RealScalar tol = tol_error;
57 VectorType r = rhs - mat * x;
60 RealScalar r0_sqnorm = r0.squaredNorm();
65 VectorType v = VectorType::Zero(n),
p = VectorType::Zero(n);
66 VectorType
y(n), z(n);
67 VectorType kt(n), ks(n);
69 VectorType s(n), t(n);
71 RealScalar tol2 = tol*tol;
74 while ( r.squaredNorm()/r0_sqnorm > tol2 && i<maxIters )
79 if (rho == Scalar(0))
return false;
80 Scalar beta = (rho/rho_old) * (alpha / w);
81 p = r + beta * (
p - w * v);
85 v.noalias() = mat *
y;
87 alpha = rho / r0.dot(v);
91 t.noalias() = mat * z;
93 w = t.dot(s) / t.squaredNorm();
94 x += alpha * y + w * z;
98 tol_error =
sqrt(r.squaredNorm()/r0_sqnorm);
105 template<
typename _MatrixType,
106 typename _Preconditioner = DiagonalPreconditioner<typename _MatrixType::Scalar> >
111 template<
typename _MatrixType,
typename _Preconditioner>
112 struct traits<BiCGSTAB<_MatrixType,_Preconditioner> >
114 typedef _MatrixType MatrixType;
115 typedef _Preconditioner Preconditioner;
166 template<
typename _MatrixType,
typename _Preconditioner>
177 typedef typename MatrixType::Scalar
Scalar;
178 typedef typename MatrixType::Index
Index;
206 template<
typename Rhs,
typename Guess>
207 inline const internal::solve_retval_with_guess<BiCGSTAB, Rhs, Guess>
212 &&
"BiCGSTAB::solve(): invalid number of rows of the right hand side matrix b");
213 return internal::solve_retval_with_guess
214 <
BiCGSTAB, Rhs, Guess>(*
this, b.derived(), x0);
218 template<
typename Rhs,
typename Dest>
222 for(
int j=0; j<b.cols(); ++j)
227 typename Dest::ColXpr xj(x,j);
238 template<
typename Rhs,
typename Dest>
252 template<
typename _MatrixType,
typename _Preconditioner,
typename Rhs>
253 struct solve_retval<BiCGSTAB<_MatrixType, _Preconditioner>, Rhs>
254 : solve_retval_base<BiCGSTAB<_MatrixType, _Preconditioner>, Rhs>
256 typedef BiCGSTAB<_MatrixType, _Preconditioner> Dec;
259 template<typename Dest>
void evalTo(Dest& dst)
const
261 dec()._solve(rhs(),dst);
269 #endif // EIGEN_BICGSTAB_H