All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
mixed_precision_cg.hpp
Go to the documentation of this file.
#ifndef VIENNACL_LINALG_MIXED_PRECISION_CG_HPP_
#define VIENNACL_LINALG_MIXED_PRECISION_CG_HPP_
/* =========================================================================
Copyright (c) 2010-2016, Institute for Microelectronics,
Institute for Analysis and Scientific Computing,
TU Wien.
Portions of this software are copyright by UChicago Argonne, LLC.
-----------------
ViennaCL - The Vienna Computing Library
-----------------
Project Head: Karl Rupp rupp@iue.tuwien.ac.at
(A list of authors and contributors can be found in the manual)
License: MIT (X11), see file LICENSE in the base directory
============================================================================= */
#include <vector>
#include <map>
#include <cmath>
#include "viennacl/forwards.h"
#include "viennacl/tools/tools.hpp"
#include "viennacl/linalg/ilu.hpp"
#include "viennacl/linalg/prod.hpp"
#include "viennacl/linalg/inner_prod.hpp"
#include "viennacl/traits/clear.hpp"
#include "viennacl/traits/size.hpp"
#include "viennacl/meta/result_of.hpp"
#include "viennacl/backend/memory.hpp"
#include "viennacl/vector_proxy.hpp"
namespace viennacl
{
namespace linalg
{
class mixed_precision_cg_tag
{
public:
mixed_precision_cg_tag(double tol = 1e-8, unsigned int max_iterations = 300, float inner_tol = 1e-2f) : tol_(tol), iterations_(max_iterations), inner_tol_(inner_tol) {}
double tolerance() const { return tol_; }
float inner_tolerance() const { return inner_tol_; }
unsigned int max_iterations() const { return iterations_; }
unsigned int iters() const { return iters_taken_; }
void iters(unsigned int i) const { iters_taken_ = i; }
double error() const { return last_error_; }
void error(double e) const { last_error_ = e; }
private:
double tol_;
unsigned int iterations_;
float inner_tol_;
//return values from solver
mutable unsigned int iters_taken_;
mutable double last_error_;
};
template<typename MatrixType, typename VectorType>
VectorType solve(const MatrixType & matrix, VectorType const & rhs, mixed_precision_cg_tag const & tag)
{
//typedef typename VectorType::value_type ScalarType;
typedef typename viennacl::result_of::cpu_value_type<VectorType>::type CPU_ScalarType;
//std::cout << "Starting CG" << std::endl;
vcl_size_t problem_size = viennacl::traits::size(rhs);
VectorType result(rhs);
viennacl::traits::clear(result);
VectorType residual = rhs;
CPU_ScalarType ip_rr = viennacl::linalg::inner_prod(rhs, rhs);
CPU_ScalarType new_ip_rr = 0;
CPU_ScalarType norm_rhs_squared = ip_rr;
if (norm_rhs_squared <= 0) //solution is zero if RHS norm is zero
return result;
viennacl::vector<float> residual_low_precision(problem_size, viennacl::traits::context(rhs));
viennacl::vector<float> result_low_precision(problem_size, viennacl::traits::context(rhs));
viennacl::vector<float> p_low_precision(problem_size, viennacl::traits::context(rhs));
viennacl::vector<float> tmp_low_precision(problem_size, viennacl::traits::context(rhs));
float inner_ip_rr = static_cast<float>(ip_rr);
float new_inner_ip_rr = 0;
float initial_inner_rhs_norm_squared = static_cast<float>(ip_rr);
float alpha;
float beta;
// transfer rhs to single precision:
p_low_precision = rhs;
residual_low_precision = p_low_precision;
// transfer matrix to single precision:
viennacl::compressed_matrix<float> matrix_low_precision(matrix.size1(), matrix.size2(), matrix.nnz(), viennacl::traits::context(rhs));
viennacl::backend::memory_copy(matrix.handle1(), const_cast<viennacl::backend::mem_handle &>(matrix_low_precision.handle1()), 0, 0, matrix_low_precision.handle1().raw_size() );
viennacl::backend::memory_copy(matrix.handle2(), const_cast<viennacl::backend::mem_handle &>(matrix_low_precision.handle2()), 0, 0, matrix_low_precision.handle2().raw_size() );
viennacl::vector_base<CPU_ScalarType> matrix_elements_high_precision(const_cast<viennacl::backend::mem_handle &>(matrix.handle()), matrix.nnz(), 0, 1);
viennacl::vector_base<float> matrix_elements_low_precision(matrix_low_precision.handle(), matrix.nnz(), 0, 1);
matrix_elements_low_precision = matrix_elements_high_precision;
matrix_low_precision.generate_row_block_information();
for (unsigned int i = 0; i < tag.max_iterations(); ++i)
{
tag.iters(i+1);
// lower precision 'inner iteration'
tmp_low_precision = viennacl::linalg::prod(matrix_low_precision, p_low_precision);
alpha = inner_ip_rr / viennacl::linalg::inner_prod(tmp_low_precision, p_low_precision);
result_low_precision += alpha * p_low_precision;
residual_low_precision -= alpha * tmp_low_precision;
new_inner_ip_rr = viennacl::linalg::inner_prod(residual_low_precision, residual_low_precision);
beta = new_inner_ip_rr / inner_ip_rr;
inner_ip_rr = new_inner_ip_rr;
p_low_precision = residual_low_precision + beta * p_low_precision;
//
// If enough progress has been achieved, update current residual with high precision evaluation
// This is effectively a restart of the CG method
//
if (new_inner_ip_rr < tag.inner_tolerance() * initial_inner_rhs_norm_squared || i == tag.max_iterations()-1)
{
residual = result_low_precision; // reusing residual vector as temporary buffer for conversion. Overwritten below anyway
result += residual;
// residual = b - Ax (without introducing a temporary)
residual = viennacl::linalg::prod(matrix, result);
residual = rhs - residual;
new_ip_rr = viennacl::linalg::inner_prod(residual, residual);
if (new_ip_rr / norm_rhs_squared < tag.tolerance() * tag.tolerance())//squared norms involved here
break;
p_low_precision = residual;
result_low_precision.clear();
residual_low_precision = p_low_precision;
initial_inner_rhs_norm_squared = static_cast<float>(new_ip_rr);
inner_ip_rr = static_cast<float>(new_ip_rr);
}
}
//store last error estimate:
tag.error(std::sqrt(new_ip_rr / norm_rhs_squared));
return result;
}
template<typename MatrixType, typename VectorType>
VectorType solve(const MatrixType & matrix, VectorType const & rhs, mixed_precision_cg_tag const & tag, viennacl::linalg::no_precond)
{
return solve(matrix, rhs, tag);
}
}
}
#endif
size.hpp
Generic size and resize functionality for different vector and matrix types.
prod.hpp
Generic interface for matrix-vector and matrix-matrix products. See viennacl/linalg/vector_operations...
tools.hpp
Various little tools used here and there in ViennaCL.
viennacl::linalg::mixed_precision_cg_tag::inner_tolerance
float inner_tolerance() const
Returns the relative tolerance.
Definition: mixed_precision_cg.hpp:61
viennacl::traits::clear
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:43
forwards.h
This file provides the forward declarations for the main types used within ViennaCL.
viennacl::matrix
A dense matrix class.
Definition: forwards.h:375
viennacl::linalg::inner_prod
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:100
inner_prod.hpp
Generic interface for the computation of inner products. See viennacl/linalg/vector_operations.hpp for implementations.
viennacl::linalg::solve
VectorT solve(MatrixT const &matrix, VectorT const &rhs, bicgstab_tag const &tag, PreconditionerT const &precond)
Definition: bicgstab.hpp:496
viennacl::linalg::prod
VectorT prod(std::vector< std::vector< T, A1 >, A2 > const &matrix, VectorT const &vector)
Definition: prod.hpp:102
viennacl::traits::size
vcl_size_t size(VectorType const &vec)
Generic routine for obtaining the size of a vector (ViennaCL, uBLAS, etc.)
Definition: size.hpp:239
viennacl::linalg::no_precond
A tag class representing the use of no preconditioner.
Definition: forwards.h:873
ilu.hpp
Implementations of incomplete factorization preconditioners. Convenience header file.
viennacl::linalg::mixed_precision_cg_tag::max_iterations
unsigned int max_iterations() const
Returns the maximum number of iterations.
Definition: mixed_precision_cg.hpp:63
viennacl::linalg::mixed_precision_cg_tag::tolerance
double tolerance() const
Returns the relative tolerance.
Definition: mixed_precision_cg.hpp:59
viennacl::vector_base
Common base class for dense vectors, vector ranges, and vector slices.
Definition: vector_def.hpp:104
viennacl::vcl_size_t
std::size_t vcl_size_t
Definition: forwards.h:75
clear.hpp
Generic clear functionality for different vector and matrix types.
viennacl::vector
Definition: forwards.h:266
viennacl::linalg::mixed_precision_cg_tag::mixed_precision_cg_tag
mixed_precision_cg_tag(double tol=1e-8, unsigned int max_iterations=300, float inner_tol=1e-2f)
The constructor.
Definition: mixed_precision_cg.hpp:56
viennacl::linalg::mixed_precision_cg_tag::error
double error() const
Returns the estimated relative error at the end of the solver run.
Definition: mixed_precision_cg.hpp:70
viennacl::result_of::cpu_value_type::type
T::ERROR_CANNOT_DEDUCE_CPU_SCALAR_TYPE_FOR_T type
Definition: result_of.hpp:271
vector_proxy.hpp
Proxy classes for vectors.
viennacl::backend::memory_copy
void memory_copy(mem_handle const &src_buffer, mem_handle &dst_buffer, vcl_size_t src_offset, vcl_size_t dst_offset, vcl_size_t bytes_to_copy)
Copies 'bytes_to_copy' bytes from address 'src_buffer + src_offset' to memory starting at address 'ds...
Definition: memory.hpp:140
viennacl::vector_base< NumericT >::clear
void clear()
Resets all entries to zero. Does not change the size of the vector.
viennacl::traits::context
viennacl::context context(T const &t)
Returns an ID for the currently active memory domain of an object.
Definition: context.hpp:40
viennacl::backend::mem_handle
Main abstraction class for multiple memory domains. Represents a buffer in either main RAM...
Definition: mem_handle.hpp:89
viennacl::linalg::mixed_precision_cg_tag::iters
void iters(unsigned int i) const
Definition: mixed_precision_cg.hpp:67
viennacl::backend::mem_handle::raw_size
vcl_size_t raw_size() const
Returns the number of bytes of the currently active buffer.
Definition: mem_handle.hpp:230
viennacl::compressed_matrix< float >
viennacl::linalg::mixed_precision_cg_tag
A tag for the conjugate gradient Used for supplying solver parameters and for dispatching the solve()...
Definition: mixed_precision_cg.hpp:47
viennacl::linalg::mixed_precision_cg_tag::error
void error(double e) const
Sets the estimated relative error at the end of the solver run.
Definition: mixed_precision_cg.hpp:72
result_of.hpp
A collection of compile time type deductions.
viennacl::linalg::mixed_precision_cg_tag::iters
unsigned int iters() const
Return the number of solver iterations:
Definition: mixed_precision_cg.hpp:66
memory.hpp
Main interface routines for memory management.