1 #ifndef VIENNACL_BICGSTAB_HPP_
2 #define VIENNACL_BICGSTAB_HPP_
57 unsigned int iters()
const {
return iters_taken_; }
58 void iters(
unsigned int i)
const { iters_taken_ = i; }
61 double error()
const {
return last_error_; }
63 void error(
double e)
const { last_error_ = e; }
67 unsigned int _iterations;
70 mutable unsigned int iters_taken_;
71 mutable double last_error_;
84 template <
typename MatrixType,
typename VectorType>
90 VectorType result(problem_size);
93 VectorType residual = rhs;
95 VectorType r0star = rhs;
96 VectorType tmp0(problem_size);
97 VectorType tmp1(problem_size);
98 VectorType s(problem_size);
101 CPU_ScalarType norm_rhs_host = ip_rr0star;
103 CPU_ScalarType alpha;
104 CPU_ScalarType omega;
105 ScalarType inner_prod_temp;
106 ScalarType new_ip_rr0star = 0;
114 alpha = ip_rr0star /
static_cast<CPU_ScalarType
>(inner_prod_temp);
123 omega = inner_prod_temp;
125 omega /= inner_prod_temp;
133 residual -= omega*tmp1;
140 CPU_ScalarType cpu_temp = new_ip_rr0star;
141 beta = cpu_temp / ip_rr0star * alpha/omega;
142 ip_rr0star = cpu_temp;
158 template <
typename MatrixType,
typename VectorType>
161 return solve(matrix, rhs, tag);
174 template <
typename MatrixType,
typename VectorType,
typename PreconditionerType>
175 VectorType
solve(
const MatrixType &
matrix, VectorType
const & rhs,
bicgstab_tag const & tag, PreconditionerType
const & precond)
180 VectorType result(problem_size);
183 VectorType residual = rhs;
184 precond.apply(residual);
185 VectorType r0star = residual;
186 VectorType tmp0(problem_size);
187 VectorType tmp1(problem_size);
188 VectorType s(problem_size);
190 VectorType p = residual;
193 CPU_ScalarType norm_rhs_host = ip_rr0star;
195 CPU_ScalarType alpha;
196 CPU_ScalarType omega;
197 ScalarType new_ip_rr0star = 0;
198 ScalarType inner_prod_temp;
207 alpha = ip_rr0star /
static_cast<CPU_ScalarType
>(inner_prod_temp);
217 omega = inner_prod_temp;
219 omega /= inner_prod_temp;
226 residual -= omega*tmp1;
233 CPU_ScalarType cpu_temp = new_ip_rr0star;
234 beta = cpu_temp / ip_rr0star * alpha/omega;
235 ip_rr0star = cpu_temp;