iterative-custom.cpp

This tutorial explains the use of iterative solvers in ViennaCL with custom monitors and initial guesses.

Note

iterative-custom.cpp and iterative-custom.cu are identical, the latter being required for compilation using CUDA nvcc

We start with including the necessary system headers:

//
// 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"

Defining Custom Monitors Functions for Iterative Solvers

Custom monitors for the iterative solvers require two ingredients: First, a structure holding all the auxiliary data we want to reuse in the monitor. Second, a callback function called by the solver in each iteration.

In this example we define a callback-routine for printing the current estimate for the residual and compare it with the true residual. To do so, we need to pass the system matrix, the right hand side, and the initial guess to the monitor routine, which we achieve by packing pointers to these objects into a struct:

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

The actual callback-routine takes the current approximation to the result as the first parameter, and the current estimate for the relative residual norm as second argument. The third argument is a pointer to our user-data, which in a first step cast to the correct type. If the monitor returns true, the iterative solver stops. This is handy for defining custom termination criteria, e.g. one-norms for the result change. Since we do not want to terminate the iterative solver with a custom criterion here, we always return 'false' at the end of the function.

Note to type-safety evangelists: This void*-interface is designed to be later exposed through a shared library ('libviennacl'). Thus, user types may not be known at the point of compilation, requiring a void*-approach.

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
}

The main Program

In the main() routine we create matrices and vectors and fill them with data.

int main()
{
  typedef float       ScalarType;

Read system matrix from a matrix-market file

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

Set up right hand side and reference solution consisting of all ones:

std::size_t size = A.size1();
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);

As initial guess we take a vector consisting of all 0.9s, except for the first entry, which we set to zero:

viennacl::vector<ScalarType> init_guess = viennacl::scalar_vector<ScalarType>(ref_result.size(), ScalarType(0.9));
init_guess[0] = 0;

Set up the monitor data, holding the system matrix, the right hand side, and the initial guess:

monitor_user_data<viennacl::compressed_matrix<ScalarType>, viennacl::vector<ScalarType> > my_monitor_data(A, b, init_guess);

set up Jacobi preconditioners (just for demonstration purposes, can be any other preconditioner here):

viennacl::linalg::jacobi_precond< viennacl::compressed_matrix<ScalarType> > jacobi(A, viennacl::linalg::jacobi_tag());

Conjugate Gradient Solver

std::cout << "----- CG Method -----" << std::endl;

Run the CG method with a relative tolerance of 1e-5 and a maximum of 20 iterations. We instantiate the CG solver object, register the monitor callback (with data), set the initial guess, and launch the solver.

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

Stabilized BiConjugate Gradient Solver

std::cout << "----- BiCGStab Method -----" << std::endl;

Run the Jacobi-preconditioned BiCGStab method with a relative tolerance of 1e-5 and a maximum of 20 iterations. We instantiate the BiCGStab solver object, register the monitor callback (with data), set the initial guess, and launch the solver.

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

GMRES Solver

std::cout << "----- GMRES Method -----" << std::endl;

Run the unpreconditioned GMRES method with a relative tolerance of 1e-5 and a maximum of 30 iterations for a Krylov size of 10 (i.e. restart every 10 iterations). We instantiate the GMRES solver object, register the monitor callback (with data), set the initial guess, and launch the solver.

Note that the monitor in the GMRES method is only called after each restart, but not in every (inner) iteration.

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

That's it, the tutorial is completed.

  std::cout << "!!!! TUTORIAL COMPLETED SUCCESSFULLY !!!!" << std::endl;

  return EXIT_SUCCESS;
}

Full Example Code

/* =========================================================================
   Copyright (c) 2010-2015, 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);

  std::size_t size = A.size1();
  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;
}