1 #ifndef VIENNACL_CG_HPP_
2 #define VIENNACL_CG_HPP_
59 unsigned int iters()
const {
return iters_taken_; }
60 void iters(
unsigned int i)
const { iters_taken_ = i; }
63 double error()
const {
return last_error_; }
65 void error(
double e)
const { last_error_ = e; }
70 unsigned int _iterations;
73 mutable unsigned int iters_taken_;
74 mutable double last_error_;
87 template <
typename MatrixType,
typename VectorType>
95 VectorType result(problem_size);
98 VectorType residual = rhs;
100 VectorType tmp(problem_size);
103 ScalarType residual_norm_squared;
105 CPU_ScalarType alpha;
106 CPU_ScalarType new_ip_rr = 0;
108 CPU_ScalarType norm_rhs_squared = ip_rr;
119 alpha = ip_rr /
static_cast<CPU_ScalarType
>(tmp_in_p);
122 residual -= alpha * tmp;
125 new_ip_rr =
static_cast<CPU_ScalarType
>(residual_norm_squared);
129 beta = new_ip_rr / ip_rr;
138 tag.
error(sqrt(new_ip_rr / norm_rhs_squared));
143 template <
typename MatrixType,
typename VectorType>
146 return solve(matrix, rhs, tag);
159 template <
typename MatrixType,
typename VectorType,
typename PreconditionerType>
160 VectorType
solve(
const MatrixType &
matrix, VectorType
const & rhs,
cg_tag const & tag, PreconditionerType
const & precond)
166 VectorType result(problem_size);
169 VectorType residual = rhs;
170 VectorType tmp(problem_size);
177 ScalarType residual_in_z;
179 CPU_ScalarType alpha;
180 CPU_ScalarType new_ip_rr = 0;
182 CPU_ScalarType norm_rhs_squared = ip_rr;
183 CPU_ScalarType new_ipp_rr_over_norm_rhs;
191 alpha = ip_rr /
static_cast<CPU_ScalarType
>(tmp_in_p);
194 residual -= alpha * tmp;
199 new_ip_rr =
static_cast<CPU_ScalarType
>(residual_in_z);
200 new_ipp_rr_over_norm_rhs = new_ip_rr / norm_rhs_squared;
204 beta = new_ip_rr / ip_rr;
213 tag.
error(sqrt(std::fabs(new_ip_rr / norm_rhs_squared)));