All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
iterative-custom.cpp
Go to the documentation of this file.
/* =========================================================================
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 PDF manual)
License: MIT (X11), see file LICENSE in the base directory
============================================================================= */
//
// include necessary system headers
//
#include <iostream>
//
// ViennaCL includes
//
#include "viennacl/scalar.hpp"
#include "viennacl/vector.hpp"
#include "viennacl/compressed_matrix.hpp"
#include "viennacl/linalg/prod.hpp"
#include "viennacl/linalg/jacobi_precond.hpp"
#include "viennacl/linalg/cg.hpp"
#include "viennacl/linalg/bicgstab.hpp"
#include "viennacl/linalg/gmres.hpp"
#include "viennacl/io/matrix_market.hpp"
template<typename MatrixT, typename VectorT>
struct monitor_user_data
{
monitor_user_data(MatrixT const & A, VectorT const & b, VectorT const & guess) : A_ptr(&A), b_ptr(&b), guess_ptr(&guess) {}
MatrixT const *A_ptr;
VectorT const *b_ptr;
VectorT const *guess_ptr;
};
template<typename VectorT, typename NumericT, typename MatrixT>
bool my_custom_monitor(VectorT const & current_approx, NumericT residual_estimate, void *user_data)
{
// Extract residual:
monitor_user_data<MatrixT, VectorT> const *data = reinterpret_cast<monitor_user_data<MatrixT, VectorT> const*>(user_data);
// Form residual r = b - A*x, taking an initial guess into account: r = b - A * (current_approx + x_initial)
VectorT x = current_approx + *data->guess_ptr;
VectorT residual = *data->b_ptr - viennacl::linalg::prod(*data->A_ptr, x);
VectorT initial_residual = *data->b_ptr - viennacl::linalg::prod(*data->A_ptr, *data->guess_ptr);
std::cout << "Residual estimate vs. true residual: " << residual_estimate << " vs. " << viennacl::linalg::norm_2(residual) / viennacl::linalg::norm_2(initial_residual) << std::endl;
return false; // no termination of iteration
}
int main()
{
typedef float ScalarType;
std::vector<std::map<unsigned int, ScalarType> > stl_A;
if (!viennacl::io::read_matrix_market_file(stl_A, "../examples/testdata/mat65k.mtx"))
{
std::cout << "Error reading Matrix file" << std::endl;
return EXIT_FAILURE;
}
viennacl::compressed_matrix<ScalarType> A;
viennacl::copy(stl_A, A);
viennacl::vector<ScalarType> ref_result = viennacl::scalar_vector<ScalarType>(A.size2(), ScalarType(1.0));
viennacl::vector<ScalarType> result(A.size2());
viennacl::vector<ScalarType> b = viennacl::linalg::prod(A, ref_result);
viennacl::vector<ScalarType> init_guess = viennacl::scalar_vector<ScalarType>(ref_result.size(), ScalarType(0.9));
init_guess[0] = 0;
monitor_user_data<viennacl::compressed_matrix<ScalarType>, viennacl::vector<ScalarType> > my_monitor_data(A, b, init_guess);
viennacl::linalg::jacobi_precond< viennacl::compressed_matrix<ScalarType> > jacobi(A, viennacl::linalg::jacobi_tag());
std::cout << "----- CG Method -----" << std::endl;
viennacl::linalg::cg_tag my_cg_tag(1e-5, 20);
viennacl::linalg::cg_solver<viennacl::vector<ScalarType> > my_cg_solver(my_cg_tag);
my_cg_solver.set_monitor(my_custom_monitor<viennacl::vector<ScalarType>, ScalarType, viennacl::compressed_matrix<ScalarType> >, &my_monitor_data);
my_cg_solver.set_initial_guess(init_guess);
my_cg_solver(A, b); // without preconditioner
std::cout << "----- BiCGStab Method -----" << std::endl;
viennacl::linalg::bicgstab_tag my_bicgstab_tag(1e-5, 20);
viennacl::linalg::bicgstab_solver<viennacl::vector<ScalarType> > my_bicgstab_solver(my_bicgstab_tag);
my_bicgstab_solver.set_monitor(my_custom_monitor<viennacl::vector<ScalarType>, ScalarType, viennacl::compressed_matrix<ScalarType> >, &my_monitor_data);
my_bicgstab_solver.set_initial_guess(init_guess);
my_bicgstab_solver(A, b, jacobi); // with Jacobi preconditioner
std::cout << "----- GMRES Method -----" << std::endl;
viennacl::linalg::gmres_tag my_gmres_tag(1e-5, 30, 10);
viennacl::linalg::gmres_solver<viennacl::vector<ScalarType> > my_gmres_solver(my_gmres_tag);
my_gmres_solver.set_monitor(my_custom_monitor<viennacl::vector<ScalarType>, ScalarType, viennacl::compressed_matrix<ScalarType> >, &my_monitor_data);
my_gmres_solver.set_initial_guess(init_guess);
my_gmres_solver(A, b);
std::cout << "!!!! TUTORIAL COMPLETED SUCCESSFULLY !!!!" << std::endl;
return EXIT_SUCCESS;
}
viennacl::compressed_matrix::size2
const vcl_size_t & size2() const
Returns the number of columns.
Definition: compressed_matrix.hpp:929
viennacl::linalg::norm_2
T norm_2(std::vector< T, A > const &v1)
Definition: norm_2.hpp:96
matrix_market.hpp
A reader and writer for the matrix market format is implemented here.
viennacl::linalg::jacobi_precond
Jacobi preconditioner class, can be supplied to solve()-routines. Generic version for non-ViennaCL ma...
Definition: jacobi_precond.hpp:50
prod.hpp
Generic interface for matrix-vector and matrix-matrix products. See viennacl/linalg/vector_operations...
main
int main()
Definition: bisect.cpp:91
bicgstab.hpp
The stabilized bi-conjugate gradient method is implemented here.
viennacl::linalg::cg_solver
Definition: cg.hpp:380
viennacl::linalg::bicgstab_solver
Definition: bicgstab.hpp:538
viennacl::linalg::jacobi_tag
A tag for a jacobi preconditioner.
Definition: jacobi_precond.hpp:43
jacobi_precond.hpp
Implementation of a simple Jacobi preconditioner.
viennacl::linalg::gmres_solver
Definition: gmres.hpp:678
viennacl::linalg::gmres_tag
A tag for the solver GMRES. Used for supplying solver parameters and for dispatching the solve() func...
Definition: gmres.hpp:49
NumericT
float NumericT
Definition: bisect.cpp:40
viennacl::linalg::prod
VectorT prod(std::vector< std::vector< T, A1 >, A2 > const &matrix, VectorT const &vector)
Definition: prod.hpp:102
gmres.hpp
Implementations of the generalized minimum residual method are in this file.
compressed_matrix.hpp
Implementation of the compressed_matrix class.
cg.hpp
The conjugate gradient method is implemented here.
viennacl::vector< ScalarType >
vector.hpp
The vector type with operator-overloads and proxy classes is defined here. Linear algebra operations ...
viennacl::scalar_vector
Represents a vector consisting of scalars 's' only, i.e. v[i] = s for all i. To be used as an initial...
Definition: vector_def.hpp:87
viennacl::copy
void copy(std::vector< NumericT > &cpu_vec, circulant_matrix< NumericT, AlignmentV > &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
viennacl::vector_base::size
size_type size() const
Returns the length of the vector (cf. std::vector)
Definition: vector_def.hpp:118
ScalarType
float ScalarType
Definition: fft_1d.cpp:42
viennacl::linalg::cg_tag
A tag for the conjugate gradient Used for supplying solver parameters and for dispatching the solve()...
Definition: cg.hpp:48
viennacl::compressed_matrix
A sparse square matrix in compressed sparse rows format.
Definition: compressed_matrix.hpp:559
viennacl::linalg::bicgstab_tag
A tag for the stabilized Bi-conjugate gradient solver. Used for supplying solver parameters and for d...
Definition: bicgstab.hpp:47
viennacl::io::read_matrix_market_file
long read_matrix_market_file(MatrixT &mat, const char *file, long index_base=1)
Reads a sparse matrix from a file (MatrixMarket format)
Definition: matrix_market.hpp:339
scalar.hpp
Implementation of the ViennaCL scalar class.