1 #ifndef VIENNACL_JACOBI_PRECOND_HPP_
2 #define VIENNACL_JACOBI_PRECOND_HPP_
45 template <
typename MatrixType>
48 typedef typename MatrixType::value_type ScalarType;
53 assert(mat.size1() == mat.size2());
54 diag_A_inv.resize(mat.size1());
56 for (
typename MatrixType::const_iterator1 row_it = system_matrix.begin1();
57 row_it != system_matrix.end1();
60 bool diag_found =
false;
61 for (
typename MatrixType::const_iterator2 col_it = row_it.begin();
62 col_it != row_it.end();
65 if (col_it.index1() == col_it.index2())
67 diag_A_inv[col_it.index1()] =
static_cast<ScalarType
>(1.0) / *col_it;
72 throw "ViennaCL: Zero in diagonal encountered while setting up Jacobi preconditioner!";
78 template <
typename VectorType>
79 void apply(VectorType & vec)
const
81 assert(vec.size() == diag_A_inv.size());
82 for (
size_t i=0; i<vec.size(); ++i)
84 vec[i] *= diag_A_inv[i];
89 MatrixType
const & system_matrix;
90 std::vector<ScalarType> diag_A_inv;
98 template <
typename ScalarType,
unsigned int MAT_ALIGNMENT>
106 assert(system_matrix.size1() == system_matrix.size2());
146 viennacl::linalg::kernels::compressed_matrix<ScalarType, MAT_ALIGNMENT>::program_name(),
150 diag_A_inv,
static_cast<cl_uint
>(diag_A_inv.size())) );
154 template <
unsigned int ALIGNMENT>
170 MatrixType
const & system_matrix;