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;
}