23 template <
class Matrix,
class Vector,
class Preconditioner,
class Real >
25 BiCG(
const Matrix &A, Vector &x,
const Vector &b,
26 const Preconditioner &M,
int &max_iter, Real &tol)
29 Vector rho_1(1), rho_2(1), alpha(1), beta(1);
30 Vector z, ztilde, p, ptilde, q, qtilde;
39 if ((resid =
norm(r) / normb) <= tol) {
45 for (
int i = 1; i <= max_iter; i++) {
47 ztilde = M.trans_solve(rtilde);
48 rho_1(0) =
dot(z, rtilde);
50 tol =
norm(r) / normb;
58 beta(0) = rho_1(0) / rho_2(0);
60 ptilde = ztilde + beta(0) * ptilde;
63 qtilde = A.trans_mult(ptilde);
64 alpha(0) = rho_1(0) /
dot(ptilde, q);
67 rtilde -= alpha(0) * qtilde;
70 if ((resid =
norm(r) / normb) < tol) {
int BiCG(const Matrix &A, Vector &x, const Vector &b, const Preconditioner &M, int &max_iter, Real &tol)
A::value_type dot(const A &a, const B &b)
Dot product of array or expression.
A::value_type norm(const A &a)
Euklidean norm of array or expression.