1 #ifndef VIENNACL_LINALG_ROW_SCALING_HPP_
2 #define VIENNACL_LINALG_ROW_SCALING_HPP_
50 unsigned int norm()
const {
return norm_; }
59 template <
typename MatrixType>
62 typedef typename MatrixType::value_type ScalarType;
72 assert(mat.size1() == mat.size2());
73 diag_M_inv.resize(mat.size1());
75 for (
typename MatrixType::const_iterator1 row_it = system_matrix.begin1();
76 row_it != system_matrix.end1();
79 diag_M_inv[row_it.index1()];
80 for (
typename MatrixType::const_iterator2 col_it = row_it.begin();
81 col_it != row_it.end();
85 diag_M_inv[col_it.index1()] += std::fabs(*col_it);
87 diag_M_inv[col_it.index1()] += (*col_it) * (*col_it);
89 if (diag_M_inv[row_it.index1()] == 0)
90 throw "ViennaCL: Zero row encountered while setting up row scaling preconditioner!";
93 diag_M_inv[row_it.index1()] =
static_cast<ScalarType
>(1.0) / diag_M_inv[row_it.index1()];
95 diag_M_inv[row_it.index1()] =
static_cast<ScalarType
>(1.0) / std::sqrt(diag_M_inv[row_it.index1()]);
101 template <
typename VectorType>
104 assert(vec.size() == diag_M_inv.size());
105 for (
size_t i=0; i<vec.size(); ++i)
107 vec[i] *= diag_M_inv[i];
112 MatrixType
const & system_matrix;
114 std::vector<ScalarType> diag_M_inv;
122 template <
typename ScalarType,
unsigned int MAT_ALIGNMENT>
135 assert(system_matrix.size1() == system_matrix.size2());
178 if (tag_.
norm() == 1)
181 viennacl::linalg::kernels::compressed_matrix<ScalarType, MAT_ALIGNMENT>::program_name(),
185 diag_M_inv,
static_cast<cl_uint
>(diag_M_inv.size())) );
190 viennacl::linalg::kernels::compressed_matrix<ScalarType, MAT_ALIGNMENT>::program_name(),
194 diag_M_inv,
static_cast<cl_uint
>(diag_M_inv.size())) );
198 template <
unsigned int ALIGNMENT>
215 MatrixType
const & system_matrix;