1 #ifndef VIENNACL_GMRES_HPP_
2 #define VIENNACL_GMRES_HPP_
65 unsigned int ret = iterations_ / krylov_dim_;
66 if (ret > 0 && (ret * krylov_dim_ == iterations_) )
72 unsigned int iters()
const {
return iters_taken_; }
74 void iters(
unsigned int i)
const { iters_taken_ = i; }
77 double error()
const {
return last_error_; }
79 void error(
double e)
const { last_error_ = e; }
83 unsigned int iterations_;
84 unsigned int krylov_dim_;
87 mutable unsigned int iters_taken_;
88 mutable double last_error_;
94 template <
typename SRC_VECTOR,
typename DEST_VECTOR>
101 template <
typename ScalarType,
typename DEST_VECTOR>
106 src.
begin() +
static_cast<difference_type
>(
start + len),
107 dest.begin() +
static_cast<difference_type
>(
start));
118 template <
typename VectorType,
typename ScalarType>
121 ScalarType input_j = input_vec(j);
136 mu = std::sqrt(sigma + input_j*input_j);
138 ScalarType hh_vec_0 = (input_j <= 0) ? (input_j - mu) : (-sigma / (input_j + mu));
140 beta = ScalarType(2) * hh_vec_0 * hh_vec_0 / (sigma + hh_vec_0 * hh_vec_0);
144 hh_vec[j] = ScalarType(1);
149 template <
typename VectorType,
typename ScalarType>
153 x -= (beta * hT_in_x) * h;
168 template <
typename MatrixType,
typename VectorType,
typename PreconditionerType>
169 VectorType
solve(
const MatrixType &
matrix, VectorType
const & rhs,
gmres_tag const & tag, PreconditionerType
const & precond)
174 VectorType result = rhs;
179 krylov_dim = problem_size;
181 VectorType res = rhs;
182 VectorType v_k_tilde = rhs;
183 VectorType v_k_tilde_temp = rhs;
185 std::vector< std::vector<CPU_ScalarType> > R(krylov_dim, std::vector<CPU_ScalarType>(tag.
krylov_dim()));
186 std::vector<CPU_ScalarType> projection_rhs(krylov_dim);
188 std::vector<VectorType> householder_reflectors(krylov_dim, rhs);
189 std::vector<CPU_ScalarType> betas(krylov_dim);
198 for (
unsigned int it = 0; it <= tag.
max_restarts(); ++it)
214 tag.
error(rho_0 / norm_rhs);
222 CPU_ScalarType rho =
static_cast<CPU_ScalarType
>(1.0);
229 for (k = 0; k < krylov_dim; ++k)
241 precond.apply(v_k_tilde);
246 v_k_tilde[k-1] = CPU_ScalarType(1);
249 for (
int i = k-1; i > -1; --i)
253 precond.apply(v_k_tilde_temp);
254 v_k_tilde = v_k_tilde_temp;
257 for (
unsigned int i = 0; i < k; ++i)
264 CPU_ScalarType rho_k_k = 0;
285 projection_rhs[k] = res[k];
287 rho *= std::sin( std::acos(projection_rhs[k] / rho) );
289 if (std::fabs(rho * rho_0 / norm_rhs) < tag.
tolerance())
291 tag.
error( std::fabs(rho*rho_0 / norm_rhs) );
301 for (
int i=k-1; i>-1; --i)
303 for (
unsigned int j=i+1; j<k; ++j)
304 projection_rhs[i] -= R[j][i] * projection_rhs[j];
306 projection_rhs[i] /= R[i][i];
313 res *= projection_rhs[0];
317 for (
unsigned int i = 0; i < k-1; ++i)
318 res[i] += projection_rhs[i+1];
324 for (
int i=k-1; i>=0; --i)
333 tag.
error(std::fabs(rho*rho_0 / norm_rhs));
343 template <
typename MatrixType,
typename VectorType>
double error() const
Returns the estimated relative error at the end of the solver run.
Definition: gmres.hpp:77
base_type::difference_type difference_type
Definition: vector.hpp:985
std::size_t vcl_size_t
Definition: forwards.h:58
T norm_2(std::vector< T, A > const &v1)
Definition: norm_2.hpp:86
Generic interface for the l^2-norm. See viennacl/linalg/vector_operations.hpp for implementations...
Generic size and resize functionality for different vector and matrix types.
Generic interface for matrix-vector and matrix-matrix products. See viennacl/linalg/vector_operations...
void error(double e) const
Sets the estimated relative error at the end of the solver run.
Definition: gmres.hpp:79
A dense matrix class.
Definition: forwards.h:293
void clear(VectorType &vec)
Generic routine for setting all entries of a vector to zero. This is the version for non-ViennaCL obj...
Definition: clear.hpp:57
This file provides the forward declarations for the main types used within ViennaCL.
unsigned int krylov_dim() const
Returns the maximum dimension of the Krylov space before restart.
Definition: gmres.hpp:61
viennacl::enable_if< viennacl::is_stl< typename viennacl::traits::tag_of< VectorT1 >::type >::value, typename VectorT1::value_type >::type inner_prod(VectorT1 const &v1, VectorT2 const &v2)
Definition: inner_prod.hpp:89
Generic interface for the computation of inner products. See viennacl/linalg/vector_operations.hpp for implementations.
A tag for the solver GMRES. Used for supplying solver parameters and for dispatching the solve() func...
Definition: gmres.hpp:44
double tolerance() const
Returns the relative tolerance.
Definition: gmres.hpp:57
Main namespace in ViennaCL. Holds all the basic types such as vector, matrix, etc. and defines operations upon them.
Definition: cpu_ram.hpp:29
VectorT prod(std::vector< std::vector< T, A1 >, A2 > const &matrix, VectorT const &vector)
Definition: prod.hpp:91
vcl_size_t size(VectorType const &vec)
Generic routine for obtaining the size of a vector (ViennaCL, uBLAS, etc.)
Definition: size.hpp:144
A tag class representing the use of no preconditioner.
Definition: forwards.h:727
result_of::size_type< T >::type start(T const &obj)
Definition: start.hpp:43
void copy(std::vector< SCALARTYPE > &cpu_vec, circulant_matrix< SCALARTYPE, ALIGNMENT > &gpu_mat)
Copies a circulant matrix from the std::vector to the OpenCL device (either GPU or multi-core CPU) ...
Definition: circulant_matrix.hpp:150
iterator begin()
Returns an iterator pointing to the beginning of the vector (STL like)
Definition: vector.hpp:803
T::value_type type
Definition: result_of.hpp:230
void gmres_householder_reflect(VectorType &x, VectorType const &h, ScalarType beta)
Definition: gmres.hpp:150
Generic clear functionality for different vector and matrix types.
T::ERROR_CANNOT_DEDUCE_CPU_SCALAR_TYPE_FOR_T type
Definition: result_of.hpp:276
void gmres_copy_helper(SRC_VECTOR const &src, DEST_VECTOR &dest, vcl_size_t len, vcl_size_t start=0)
Definition: gmres.hpp:95
unsigned int max_restarts() const
Returns the maximum number of GMRES restarts.
Definition: gmres.hpp:63
unsigned int max_iterations() const
Returns the maximum number of iterations.
Definition: gmres.hpp:59
VectorType solve(const MatrixType &matrix, VectorType const &rhs, bicgstab_tag const &tag)
Implementation of the stabilized Bi-conjugate gradient solver.
Definition: bicgstab.hpp:92
gmres_tag(double tol=1e-10, unsigned int max_iterations=300, unsigned int krylov_dim=20)
The constructor.
Definition: gmres.hpp:53
void iters(unsigned int i) const
Set the number of solver iterations (should only be modified by the solver)
Definition: gmres.hpp:74
unsigned int iters() const
Return the number of solver iterations:
Definition: gmres.hpp:72
A collection of compile time type deductions.
void gmres_setup_householder_vector(VectorType const &input_vec, VectorType &hh_vec, ScalarType &beta, ScalarType &mu, vcl_size_t j)
Computes the householder vector 'hh_vec' which rotates 'input_vec' such that all entries below the j-...
Definition: gmres.hpp:119