All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
solver.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
============================================================================= */
/*
*
* Benchmark: Iterative solver tests (solver.cpp and solver.cu are identical, the latter being required for compilation using CUDA nvcc)
*
*/
#ifndef BOOST_UBLAS_NDEBUG
#define BOOST_UBLAS_NDEBUG
#endif
#include <boost/numeric/ublas/matrix_sparse.hpp>
#include <boost/numeric/ublas/io.hpp>
#include <boost/numeric/ublas/operation_sparse.hpp>
#define VIENNACL_WITH_UBLAS 1
#include "viennacl/scalar.hpp"
#include "viennacl/vector.hpp"
#include "viennacl/coordinate_matrix.hpp"
#include "viennacl/compressed_matrix.hpp"
#include "viennacl/ell_matrix.hpp"
#include "viennacl/sliced_ell_matrix.hpp"
#include "viennacl/hyb_matrix.hpp"
#include "viennacl/context.hpp"
#include "viennacl/linalg/cg.hpp"
#include "viennacl/linalg/bicgstab.hpp"
#include "viennacl/linalg/gmres.hpp"
#include "viennacl/linalg/mixed_precision_cg.hpp"
#include "viennacl/linalg/ilu.hpp"
#include "viennacl/linalg/ichol.hpp"
#include "viennacl/linalg/jacobi_precond.hpp"
#include "viennacl/linalg/row_scaling.hpp"
#include "viennacl/io/matrix_market.hpp"
#include "viennacl/tools/timer.hpp"
#include <iostream>
#include <vector>
using namespace boost::numeric;
#define BENCHMARK_RUNS 1
inline void printOps(double num_ops, double exec_time)
{
std::cout << "GFLOPs: " << num_ops / (1000000 * exec_time * 1000) << std::endl;
}
template<typename ScalarType>
ScalarType diff_inf(ublas::vector<ScalarType> & v1, viennacl::vector<ScalarType> & v2)
{
ublas::vector<ScalarType> v2_cpu(v2.size());
viennacl::copy(v2.begin(), v2.end(), v2_cpu.begin());
for (unsigned int i=0;i<v1.size(); ++i)
{
if ( std::max( fabs(v2_cpu[i]), fabs(v1[i]) ) > 0 )
v2_cpu[i] = fabs(v2_cpu[i] - v1[i]) / std::max( fabs(v2_cpu[i]), fabs(v1[i]) );
else
v2_cpu[i] = 0.0;
}
return norm_inf(v2_cpu);
}
template<typename ScalarType>
ScalarType diff_2(ublas::vector<ScalarType> & v1, viennacl::vector<ScalarType> & v2)
{
ublas::vector<ScalarType> v2_cpu(v2.size());
viennacl::copy(v2.begin(), v2.end(), v2_cpu.begin());
return norm_2(v1 - v2_cpu) / norm_2(v1);
}
template<typename MatrixType, typename VectorType, typename SolverTag, typename PrecondTag>
void run_solver(MatrixType const & matrix, VectorType const & rhs, VectorType const & ref_result, SolverTag const & solver, PrecondTag const & precond, long ops)
{
viennacl::tools::timer timer;
VectorType result(rhs);
VectorType residual(rhs);
viennacl::backend::finish();
timer.start();
for (int runs=0; runs<BENCHMARK_RUNS; ++runs)
{
result = viennacl::linalg::solve(matrix, rhs, solver, precond);
}
viennacl::backend::finish();
double exec_time = timer.get();
std::cout << "Exec. time: " << exec_time << std::endl;
std::cout << "Est. "; printOps(static_cast<double>(ops), exec_time / BENCHMARK_RUNS);
residual -= viennacl::linalg::prod(matrix, result);
std::cout << "Relative residual: " << viennacl::linalg::norm_2(residual) / viennacl::linalg::norm_2(rhs) << std::endl;
std::cout << "Estimated rel. residual: " << solver.error() << std::endl;
std::cout << "Iterations: " << solver.iters() << std::endl;
result -= ref_result;
std::cout << "Relative deviation from result: " << viennacl::linalg::norm_2(result) / viennacl::linalg::norm_2(ref_result) << std::endl;
}
template<typename ScalarType>
int run_benchmark(viennacl::context ctx)
{
viennacl::tools::timer timer;
double exec_time;
ublas::vector<ScalarType> ublas_vec1;
ublas::vector<ScalarType> ublas_vec2;
ublas::vector<ScalarType> ublas_result;
unsigned int solver_iters = 100;
unsigned int solver_krylov_dim = 20;
double solver_tolerance = 1e-6;
ublas::compressed_matrix<ScalarType> ublas_matrix;
if (!viennacl::io::read_matrix_market_file(ublas_matrix, "../examples/testdata/mat65k.mtx"))
{
std::cout << "Error reading Matrix file" << std::endl;
return EXIT_FAILURE;
}
std::cout << "done reading matrix" << std::endl;
ublas_result = ublas::scalar_vector<ScalarType>(ublas_matrix.size1(), ScalarType(1.0));
ublas_vec1 = ublas::prod(ublas_matrix, ublas_result);
ublas_vec2 = ublas_vec1;
viennacl::compressed_matrix<ScalarType> vcl_compressed_matrix(ublas_vec1.size(), ublas_vec1.size(), ctx);
viennacl::coordinate_matrix<ScalarType> vcl_coordinate_matrix(ublas_vec1.size(), ublas_vec1.size(), ctx);
viennacl::ell_matrix<ScalarType> vcl_ell_matrix(ctx);
viennacl::sliced_ell_matrix<ScalarType> vcl_sliced_ell_matrix(ctx);
viennacl::hyb_matrix<ScalarType> vcl_hyb_matrix(ctx);
viennacl::vector<ScalarType> vcl_vec1(ublas_vec1.size(), ctx);
viennacl::vector<ScalarType> vcl_vec2(ublas_vec1.size(), ctx);
viennacl::vector<ScalarType> vcl_result(ublas_vec1.size(), ctx);
//cpu to gpu:
viennacl::copy(ublas_matrix, vcl_compressed_matrix);
viennacl::copy(ublas_matrix, vcl_coordinate_matrix);
viennacl::copy(ublas_matrix, vcl_ell_matrix);
viennacl::copy(ublas_matrix, vcl_sliced_ell_matrix);
viennacl::copy(ublas_matrix, vcl_hyb_matrix);
viennacl::copy(ublas_vec1, vcl_vec1);
viennacl::copy(ublas_vec2, vcl_vec2);
viennacl::copy(ublas_result, vcl_result);
std::cout << "------- Jacobi preconditioner ----------" << std::endl;
viennacl::linalg::jacobi_precond< ublas::compressed_matrix<ScalarType> > ublas_jacobi(ublas_matrix, viennacl::linalg::jacobi_tag());
viennacl::linalg::jacobi_precond< viennacl::compressed_matrix<ScalarType> > vcl_jacobi_csr(vcl_compressed_matrix, viennacl::linalg::jacobi_tag());
viennacl::linalg::jacobi_precond< viennacl::coordinate_matrix<ScalarType> > vcl_jacobi_coo(vcl_coordinate_matrix, viennacl::linalg::jacobi_tag());
std::cout << "------- Row-Scaling preconditioner ----------" << std::endl;
viennacl::linalg::row_scaling< ublas::compressed_matrix<ScalarType> > ublas_row_scaling(ublas_matrix, viennacl::linalg::row_scaling_tag(1));
viennacl::linalg::row_scaling< viennacl::compressed_matrix<ScalarType> > vcl_row_scaling_csr(vcl_compressed_matrix, viennacl::linalg::row_scaling_tag(1));
viennacl::linalg::row_scaling< viennacl::coordinate_matrix<ScalarType> > vcl_row_scaling_coo(vcl_coordinate_matrix, viennacl::linalg::row_scaling_tag(1));
std::cout << "------- ICHOL0 on CPU (ublas) ----------" << std::endl;
timer.start();
viennacl::linalg::ichol0_precond< ublas::compressed_matrix<ScalarType> > ublas_ichol0(ublas_matrix, viennacl::linalg::ichol0_tag());
exec_time = timer.get();
std::cout << "Setup time: " << exec_time << std::endl;
timer.start();
for (int runs=0; runs<BENCHMARK_RUNS; ++runs)
ublas_ichol0.apply(ublas_vec1);
exec_time = timer.get();
std::cout << "ublas time: " << exec_time << std::endl;
std::cout << "------- ICHOL0 with ViennaCL ----------" << std::endl;
timer.start();
viennacl::linalg::ichol0_precond< viennacl::compressed_matrix<ScalarType> > vcl_ichol0(vcl_compressed_matrix, viennacl::linalg::ichol0_tag());
exec_time = timer.get();
std::cout << "Setup time: " << exec_time << std::endl;
viennacl::backend::finish();
timer.start();
for (int runs=0; runs<BENCHMARK_RUNS; ++runs)
vcl_ichol0.apply(vcl_vec1);
viennacl::backend::finish();
exec_time = timer.get();
std::cout << "ViennaCL time: " << exec_time << std::endl;
std::cout << "------- Chow-Patel parallel ICC with ViennaCL ----------" << std::endl;
timer.start();
viennacl::linalg::chow_patel_icc_precond< viennacl::compressed_matrix<ScalarType> > vcl_chow_patel_icc(vcl_compressed_matrix, viennacl::linalg::chow_patel_tag());
viennacl::backend::finish();
std::cout << "Setup time: " << timer.get() << std::endl;
timer.start();
for (int runs=0; runs<BENCHMARK_RUNS; ++runs)
vcl_chow_patel_icc.apply(vcl_vec1);
viennacl::backend::finish();
std::cout << "ViennaCL Chow-Patel-ICC substitution time: " << timer.get() << std::endl;
std::cout << "------- ILU0 on with ublas ----------" << std::endl;
timer.start();
viennacl::linalg::ilu0_precond< ublas::compressed_matrix<ScalarType> > ublas_ilu0(ublas_matrix, viennacl::linalg::ilu0_tag());
exec_time = timer.get();
std::cout << "Setup time (no level scheduling): " << exec_time << std::endl;
timer.start();
for (int runs=0; runs<BENCHMARK_RUNS; ++runs)
ublas_ilu0.apply(ublas_vec1);
exec_time = timer.get();
std::cout << "ublas ILU0 substitution time (no level scheduling): " << exec_time << std::endl;
std::cout << "------- ILU0 with ViennaCL ----------" << std::endl;
timer.start();
viennacl::linalg::ilu0_precond< viennacl::compressed_matrix<ScalarType> > vcl_ilu0(vcl_compressed_matrix, viennacl::linalg::ilu0_tag());
exec_time = timer.get();
std::cout << "Setup time (no level scheduling): " << exec_time << std::endl;
viennacl::backend::finish();
timer.start();
for (int runs=0; runs<BENCHMARK_RUNS; ++runs)
vcl_ilu0.apply(vcl_vec1);
viennacl::backend::finish();
exec_time = timer.get();
std::cout << "ViennaCL ILU0 substitution time (no level scheduling): " << exec_time << std::endl;
timer.start();
viennacl::linalg::ilu0_tag ilu0_with_level_scheduling; ilu0_with_level_scheduling.use_level_scheduling(true);
viennacl::linalg::ilu0_precond< viennacl::compressed_matrix<ScalarType> > vcl_ilu0_level_scheduling(vcl_compressed_matrix, ilu0_with_level_scheduling);
exec_time = timer.get();
std::cout << "Setup time (with level scheduling): " << exec_time << std::endl;
viennacl::backend::finish();
timer.start();
for (int runs=0; runs<BENCHMARK_RUNS; ++runs)
vcl_ilu0_level_scheduling.apply(vcl_vec1);
viennacl::backend::finish();
exec_time = timer.get();
std::cout << "ViennaCL ILU0 substitution time (with level scheduling): " << exec_time << std::endl;
std::cout << "------- Block-ILU0 with ublas ----------" << std::endl;
ublas_vec1 = ublas_vec2;
viennacl::copy(ublas_vec1, vcl_vec1);
timer.start();
viennacl::linalg::block_ilu_precond< ublas::compressed_matrix<ScalarType>,
viennacl::linalg::ilu0_tag> ublas_block_ilu0(ublas_matrix, viennacl::linalg::ilu0_tag());
exec_time = timer.get();
std::cout << "Setup time: " << exec_time << std::endl;
timer.start();
for (int runs=0; runs<BENCHMARK_RUNS; ++runs)
ublas_block_ilu0.apply(ublas_vec1);
exec_time = timer.get();
std::cout << "ublas time: " << exec_time << std::endl;
std::cout << "------- Block-ILU0 with ViennaCL ----------" << std::endl;
timer.start();
viennacl::linalg::block_ilu_precond< viennacl::compressed_matrix<ScalarType>,
viennacl::linalg::ilu0_tag> vcl_block_ilu0(vcl_compressed_matrix, viennacl::linalg::ilu0_tag());
exec_time = timer.get();
std::cout << "Setup time: " << exec_time << std::endl;
//vcl_block_ilu0.apply(vcl_vec1); //warm-up
viennacl::backend::finish();
timer.start();
for (int runs=0; runs<BENCHMARK_RUNS; ++runs)
vcl_block_ilu0.apply(vcl_vec1);
viennacl::backend::finish();
exec_time = timer.get();
std::cout << "ViennaCL time: " << exec_time << std::endl;
std::cout << "------- ILUT with ublas ----------" << std::endl;
ublas_vec1 = ublas_vec2;
viennacl::copy(ublas_vec1, vcl_vec1);
timer.start();
viennacl::linalg::ilut_precond< ublas::compressed_matrix<ScalarType> > ublas_ilut(ublas_matrix, viennacl::linalg::ilut_tag());
exec_time = timer.get();
std::cout << "Setup time (no level scheduling): " << exec_time << std::endl;
timer.start();
for (int runs=0; runs<BENCHMARK_RUNS; ++runs)
ublas_ilut.apply(ublas_vec1);
exec_time = timer.get();
std::cout << "ublas ILUT substitution time (no level scheduling): " << exec_time << std::endl;
std::cout << "------- ILUT with ViennaCL ----------" << std::endl;
timer.start();
viennacl::linalg::ilut_precond< viennacl::compressed_matrix<ScalarType> > vcl_ilut(vcl_compressed_matrix, viennacl::linalg::ilut_tag());
exec_time = timer.get();
std::cout << "Setup time (no level scheduling): " << exec_time << std::endl;
viennacl::backend::finish();
timer.start();
for (int runs=0; runs<BENCHMARK_RUNS; ++runs)
vcl_ilut.apply(vcl_vec1);
viennacl::backend::finish();
exec_time = timer.get();
std::cout << "ViennaCL ILUT substitution time (no level scheduling): " << exec_time << std::endl;
timer.start();
viennacl::linalg::ilut_tag ilut_with_level_scheduling; ilut_with_level_scheduling.use_level_scheduling(true);
viennacl::linalg::ilut_precond< viennacl::compressed_matrix<ScalarType> > vcl_ilut_level_scheduling(vcl_compressed_matrix, ilut_with_level_scheduling);
exec_time = timer.get();
std::cout << "Setup time (with level scheduling): " << exec_time << std::endl;
viennacl::backend::finish();
timer.start();
for (int runs=0; runs<BENCHMARK_RUNS; ++runs)
vcl_ilut_level_scheduling.apply(vcl_vec1);
viennacl::backend::finish();
exec_time = timer.get();
std::cout << "ViennaCL ILUT substitution time (with level scheduling): " << exec_time << std::endl;
std::cout << "------- Block-ILUT with ublas ----------" << std::endl;
ublas_vec1 = ublas_vec2;
viennacl::copy(ublas_vec1, vcl_vec1);
timer.start();
viennacl::linalg::block_ilu_precond< ublas::compressed_matrix<ScalarType>,
viennacl::linalg::ilut_tag> ublas_block_ilut(ublas_matrix, viennacl::linalg::ilut_tag());
exec_time = timer.get();
std::cout << "Setup time: " << exec_time << std::endl;
//ublas_block_ilut.apply(ublas_vec1);
timer.start();
for (int runs=0; runs<BENCHMARK_RUNS; ++runs)
ublas_block_ilut.apply(ublas_vec1);
exec_time = timer.get();
std::cout << "ublas time: " << exec_time << std::endl;
std::cout << "------- Block-ILUT with ViennaCL ----------" << std::endl;
timer.start();
viennacl::linalg::block_ilu_precond< viennacl::compressed_matrix<ScalarType>,
viennacl::linalg::ilut_tag> vcl_block_ilut(vcl_compressed_matrix, viennacl::linalg::ilut_tag());
exec_time = timer.get();
std::cout << "Setup time: " << exec_time << std::endl;
//vcl_block_ilut.apply(vcl_vec1); //warm-up
viennacl::backend::finish();
timer.start();
for (int runs=0; runs<BENCHMARK_RUNS; ++runs)
vcl_block_ilut.apply(vcl_vec1);
viennacl::backend::finish();
exec_time = timer.get();
std::cout << "ViennaCL time: " << exec_time << std::endl;
std::cout << "------- Chow-Patel parallel ILU with ViennaCL ----------" << std::endl;
timer.start();
viennacl::linalg::chow_patel_ilu_precond< viennacl::compressed_matrix<ScalarType> > vcl_chow_patel_ilu(vcl_compressed_matrix, viennacl::linalg::chow_patel_tag());
viennacl::backend::finish();
std::cout << "Setup time: " << timer.get() << std::endl;
timer.start();
for (int runs=0; runs<BENCHMARK_RUNS; ++runs)
vcl_chow_patel_ilu.apply(vcl_vec1);
viennacl::backend::finish();
std::cout << "ViennaCL Chow-Patel-ILU substitution time: " << timer.get() << std::endl;
long cg_ops = static_cast<long>(solver_iters * (ublas_matrix.nnz() + 6 * ublas_vec2.size()));
viennacl::linalg::cg_tag cg_solver(solver_tolerance, solver_iters);
std::cout << "------- CG solver (no preconditioner) using ublas ----------" << std::endl;
run_solver(ublas_matrix, ublas_vec2, ublas_result, cg_solver, viennacl::linalg::no_precond(), cg_ops);
std::cout << "------- CG solver (no preconditioner) via ViennaCL, compressed_matrix ----------" << std::endl;
run_solver(vcl_compressed_matrix, vcl_vec2, vcl_result, cg_solver, viennacl::linalg::no_precond(), cg_ops);
bool is_double = (sizeof(ScalarType) == sizeof(double));
if (is_double)
{
std::cout << "------- CG solver, mixed precision (no preconditioner) via ViennaCL, compressed_matrix ----------" << std::endl;
viennacl::linalg::mixed_precision_cg_tag mixed_precision_cg_solver(solver_tolerance, solver_iters);
run_solver(vcl_compressed_matrix, vcl_vec2, vcl_result, mixed_precision_cg_solver, viennacl::linalg::no_precond(), cg_ops);
}
std::cout << "------- CG solver (no preconditioner) via ViennaCL, coordinate_matrix ----------" << std::endl;
run_solver(vcl_coordinate_matrix, vcl_vec2, vcl_result, cg_solver, viennacl::linalg::no_precond(), cg_ops);
std::cout << "------- CG solver (no preconditioner) via ViennaCL, ell_matrix ----------" << std::endl;
run_solver(vcl_ell_matrix, vcl_vec2, vcl_result, cg_solver, viennacl::linalg::no_precond(), cg_ops);
std::cout << "------- CG solver (no preconditioner) via ViennaCL, sliced_ell_matrix ----------" << std::endl;
run_solver(vcl_sliced_ell_matrix, vcl_vec2, vcl_result, cg_solver, viennacl::linalg::no_precond(), cg_ops);
std::cout << "------- CG solver (no preconditioner) via ViennaCL, hyb_matrix ----------" << std::endl;
run_solver(vcl_hyb_matrix, vcl_vec2, vcl_result, cg_solver, viennacl::linalg::no_precond(), cg_ops);
std::cout << "------- CG solver (ICHOL0 preconditioner) using ublas ----------" << std::endl;
run_solver(ublas_matrix, ublas_vec2, ublas_result, cg_solver, ublas_ichol0, cg_ops);
std::cout << "------- CG solver (ICHOL0 preconditioner) via ViennaCL, compressed_matrix ----------" << std::endl;
run_solver(vcl_compressed_matrix, vcl_vec2, vcl_result, cg_solver, vcl_ichol0, cg_ops);
std::cout << "------- CG solver (Chow-Patel ICHOL0 preconditioner) via ViennaCL, compressed_matrix ----------" << std::endl;
run_solver(vcl_compressed_matrix, vcl_vec2, vcl_result, cg_solver, vcl_chow_patel_icc, cg_ops);
std::cout << "------- CG solver (ILU0 preconditioner) using ublas ----------" << std::endl;
run_solver(ublas_matrix, ublas_vec2, ublas_result, cg_solver, ublas_ilu0, cg_ops);
std::cout << "------- CG solver (ILU0 preconditioner) via ViennaCL, compressed_matrix ----------" << std::endl;
run_solver(vcl_compressed_matrix, vcl_vec2, vcl_result, cg_solver, vcl_ilu0, cg_ops);
std::cout << "------- CG solver (Block-ILU0 preconditioner) using ublas ----------" << std::endl;
run_solver(ublas_matrix, ublas_vec2, ublas_result, cg_solver, ublas_block_ilu0, cg_ops);
std::cout << "------- CG solver (Block-ILU0 preconditioner) via ViennaCL, compressed_matrix ----------" << std::endl;
run_solver(vcl_compressed_matrix, vcl_vec2, vcl_result, cg_solver, vcl_block_ilu0, cg_ops);
std::cout << "------- CG solver (ILUT preconditioner) using ublas ----------" << std::endl;
run_solver(ublas_matrix, ublas_vec2, ublas_result, cg_solver, ublas_ilut, cg_ops);
std::cout << "------- CG solver (ILUT preconditioner) via ViennaCL, compressed_matrix ----------" << std::endl;
run_solver(vcl_compressed_matrix, vcl_vec2, vcl_result, cg_solver, vcl_ilut, cg_ops);
std::cout << "------- CG solver (ILUT preconditioner) via ViennaCL, coordinate_matrix ----------" << std::endl;
run_solver(vcl_coordinate_matrix, vcl_vec2, vcl_result, cg_solver, vcl_ilut, cg_ops);
std::cout << "------- CG solver (Block-ILUT preconditioner) using ublas ----------" << std::endl;
run_solver(ublas_matrix, ublas_vec2, ublas_result, cg_solver, ublas_block_ilut, cg_ops);
std::cout << "------- CG solver (Block-ILUT preconditioner) via ViennaCL, compressed_matrix ----------" << std::endl;
run_solver(vcl_compressed_matrix, vcl_vec2, vcl_result, cg_solver, vcl_block_ilut, cg_ops);
std::cout << "------- CG solver (Jacobi preconditioner) using ublas ----------" << std::endl;
run_solver(ublas_matrix, ublas_vec2, ublas_result, cg_solver, ublas_jacobi, cg_ops);
std::cout << "------- CG solver (Jacobi preconditioner) via ViennaCL, compressed_matrix ----------" << std::endl;
run_solver(vcl_compressed_matrix, vcl_vec2, vcl_result, cg_solver, vcl_jacobi_csr, cg_ops);
std::cout << "------- CG solver (Jacobi preconditioner) via ViennaCL, coordinate_matrix ----------" << std::endl;
run_solver(vcl_coordinate_matrix, vcl_vec2, vcl_result, cg_solver, vcl_jacobi_coo, cg_ops);
std::cout << "------- CG solver (row scaling preconditioner) using ublas ----------" << std::endl;
run_solver(ublas_matrix, ublas_vec2, ublas_result, cg_solver, ublas_row_scaling, cg_ops);
std::cout << "------- CG solver (row scaling preconditioner) via ViennaCL, compressed_matrix ----------" << std::endl;
run_solver(vcl_compressed_matrix, vcl_vec2, vcl_result, cg_solver, vcl_row_scaling_csr, cg_ops);
std::cout << "------- CG solver (row scaling preconditioner) via ViennaCL, coordinate_matrix ----------" << std::endl;
run_solver(vcl_coordinate_matrix, vcl_vec2, vcl_result, cg_solver, vcl_row_scaling_coo, cg_ops);
long bicgstab_ops = static_cast<long>(solver_iters * (2 * ublas_matrix.nnz() + 13 * ublas_vec2.size()));
viennacl::linalg::bicgstab_tag bicgstab_solver(solver_tolerance, solver_iters);
std::cout << "------- BiCGStab solver (no preconditioner) using ublas ----------" << std::endl;
run_solver(ublas_matrix, ublas_vec2, ublas_result, bicgstab_solver, viennacl::linalg::no_precond(), bicgstab_ops);
std::cout << "------- BiCGStab solver (no preconditioner) via ViennaCL, compressed_matrix ----------" << std::endl;
run_solver(vcl_compressed_matrix, vcl_vec2, vcl_result, bicgstab_solver, viennacl::linalg::no_precond(), bicgstab_ops);
std::cout << "------- BiCGStab solver (ILU0 preconditioner) via ViennaCL, compressed_matrix ----------" << std::endl;
run_solver(vcl_compressed_matrix, vcl_vec2, vcl_result, bicgstab_solver, vcl_ilu0, bicgstab_ops);
std::cout << "------- BiCGStab solver (Chow-Patel-ILU preconditioner) via ViennaCL, compressed_matrix ----------" << std::endl;
run_solver(vcl_compressed_matrix, vcl_vec2, vcl_result, bicgstab_solver, vcl_chow_patel_ilu, bicgstab_ops);
std::cout << "------- BiCGStab solver (no preconditioner) via ViennaCL, coordinate_matrix ----------" << std::endl;
run_solver(vcl_coordinate_matrix, vcl_vec2, vcl_result, bicgstab_solver, viennacl::linalg::no_precond(), bicgstab_ops);
std::cout << "------- BiCGStab solver (no preconditioner) via ViennaCL, ell_matrix ----------" << std::endl;
run_solver(vcl_ell_matrix, vcl_vec2, vcl_result, bicgstab_solver, viennacl::linalg::no_precond(), bicgstab_ops);
std::cout << "------- BiCGStab solver (no preconditioner) via ViennaCL, sliced_ell_matrix ----------" << std::endl;
run_solver(vcl_sliced_ell_matrix, vcl_vec2, vcl_result, bicgstab_solver, viennacl::linalg::no_precond(), bicgstab_ops);
std::cout << "------- BiCGStab solver (no preconditioner) via ViennaCL, hyb_matrix ----------" << std::endl;
run_solver(vcl_hyb_matrix, vcl_vec2, vcl_result, bicgstab_solver, viennacl::linalg::no_precond(), bicgstab_ops);
std::cout << "------- BiCGStab solver (ILUT preconditioner) using ublas ----------" << std::endl;
run_solver(ublas_matrix, ublas_vec2, ublas_result, bicgstab_solver, ublas_ilut, bicgstab_ops);
std::cout << "------- BiCGStab solver (ILUT preconditioner) via ViennaCL, compressed_matrix ----------" << std::endl;
run_solver(vcl_compressed_matrix, vcl_vec2, vcl_result, bicgstab_solver, vcl_ilut, bicgstab_ops);
std::cout << "------- BiCGStab solver (Block-ILUT preconditioner) using ublas ----------" << std::endl;
run_solver(ublas_matrix, ublas_vec2, ublas_result, bicgstab_solver, ublas_block_ilut, bicgstab_ops);
#ifdef VIENNACL_WITH_OPENCL
std::cout << "------- BiCGStab solver (Block-ILUT preconditioner) via ViennaCL, compressed_matrix ----------" << std::endl;
run_solver(vcl_compressed_matrix, vcl_vec2, vcl_result, bicgstab_solver, vcl_block_ilut, bicgstab_ops);
#endif
// std::cout << "------- BiCGStab solver (ILUT preconditioner) via ViennaCL, coordinate_matrix ----------" << std::endl;
// run_solver(vcl_coordinate_matrix, vcl_vec2, vcl_result, bicgstab_solver, vcl_ilut, bicgstab_ops);
std::cout << "------- BiCGStab solver (Jacobi preconditioner) using ublas ----------" << std::endl;
run_solver(ublas_matrix, ublas_vec2, ublas_result, bicgstab_solver, ublas_jacobi, bicgstab_ops);
std::cout << "------- BiCGStab solver (Jacobi preconditioner) via ViennaCL, compressed_matrix ----------" << std::endl;
run_solver(vcl_compressed_matrix, vcl_vec2, vcl_result, bicgstab_solver, vcl_jacobi_csr, bicgstab_ops);
std::cout << "------- BiCGStab solver (Jacobi preconditioner) via ViennaCL, coordinate_matrix ----------" << std::endl;
run_solver(vcl_coordinate_matrix, vcl_vec2, vcl_result, bicgstab_solver, vcl_jacobi_coo, bicgstab_ops);
std::cout << "------- BiCGStab solver (row scaling preconditioner) using ublas ----------" << std::endl;
run_solver(ublas_matrix, ublas_vec2, ublas_result, bicgstab_solver, ublas_row_scaling, bicgstab_ops);
std::cout << "------- BiCGStab solver (row scaling preconditioner) via ViennaCL, compressed_matrix ----------" << std::endl;
run_solver(vcl_compressed_matrix, vcl_vec2, vcl_result, bicgstab_solver, vcl_row_scaling_csr, bicgstab_ops);
std::cout << "------- BiCGStab solver (row scaling preconditioner) via ViennaCL, coordinate_matrix ----------" << std::endl;
run_solver(vcl_coordinate_matrix, vcl_vec2, vcl_result, bicgstab_solver, vcl_row_scaling_coo, bicgstab_ops);
long gmres_ops = static_cast<long>(solver_iters * (ublas_matrix.nnz() + (solver_iters * 2 + 7) * ublas_vec2.size()));
viennacl::linalg::gmres_tag gmres_solver(solver_tolerance, solver_iters, solver_krylov_dim);
std::cout << "------- GMRES solver (no preconditioner) using ublas ----------" << std::endl;
run_solver(ublas_matrix, ublas_vec2, ublas_result, gmres_solver, viennacl::linalg::no_precond(), gmres_ops);
std::cout << "------- GMRES solver (no preconditioner) via ViennaCL, compressed_matrix ----------" << std::endl;
run_solver(vcl_compressed_matrix, vcl_vec2, vcl_result, gmres_solver, viennacl::linalg::no_precond(), gmres_ops);
std::cout << "------- GMRES solver (no preconditioner) on GPU, coordinate_matrix ----------" << std::endl;
run_solver(vcl_coordinate_matrix, vcl_vec2, vcl_result, gmres_solver, viennacl::linalg::no_precond(), gmres_ops);
std::cout << "------- GMRES solver (no preconditioner) on GPU, ell_matrix ----------" << std::endl;
run_solver(vcl_ell_matrix, vcl_vec2, vcl_result, gmres_solver, viennacl::linalg::no_precond(), gmres_ops);
std::cout << "------- GMRES solver (no preconditioner) on GPU, sliced_ell_matrix ----------" << std::endl;
run_solver(vcl_sliced_ell_matrix, vcl_vec2, vcl_result, gmres_solver, viennacl::linalg::no_precond(), gmres_ops);
std::cout << "------- GMRES solver (no preconditioner) on GPU, hyb_matrix ----------" << std::endl;
run_solver(vcl_hyb_matrix, vcl_vec2, vcl_result, gmres_solver, viennacl::linalg::no_precond(), gmres_ops);
std::cout << "------- GMRES solver (ILUT preconditioner) using ublas ----------" << std::endl;
run_solver(ublas_matrix, ublas_vec2, ublas_result, gmres_solver, ublas_ilut, gmres_ops);
std::cout << "------- GMRES solver (ILUT preconditioner) via ViennaCL, compressed_matrix ----------" << std::endl;
run_solver(vcl_compressed_matrix, vcl_vec2, vcl_result, gmres_solver, vcl_ilut, gmres_ops);
std::cout << "------- GMRES solver (ILUT preconditioner) via ViennaCL, coordinate_matrix ----------" << std::endl;
run_solver(vcl_coordinate_matrix, vcl_vec2, vcl_result, gmres_solver, vcl_ilut, gmres_ops);
std::cout << "------- GMRES solver (Jacobi preconditioner) using ublas ----------" << std::endl;
run_solver(ublas_matrix, ublas_vec2, ublas_result, gmres_solver, ublas_jacobi, gmres_ops);
std::cout << "------- GMRES solver (Jacobi preconditioner) via ViennaCL, compressed_matrix ----------" << std::endl;
run_solver(vcl_compressed_matrix, vcl_vec2, vcl_result, gmres_solver, vcl_jacobi_csr, gmres_ops);
std::cout << "------- GMRES solver (Jacobi preconditioner) via ViennaCL, coordinate_matrix ----------" << std::endl;
run_solver(vcl_coordinate_matrix, vcl_vec2, vcl_result, gmres_solver, vcl_jacobi_coo, gmres_ops);
std::cout << "------- GMRES solver (row scaling preconditioner) using ublas ----------" << std::endl;
run_solver(ublas_matrix, ublas_vec2, ublas_result, gmres_solver, ublas_row_scaling, gmres_ops);
std::cout << "------- GMRES solver (row scaling preconditioner) via ViennaCL, compressed_matrix ----------" << std::endl;
run_solver(vcl_compressed_matrix, vcl_vec2, vcl_result, gmres_solver, vcl_row_scaling_csr, gmres_ops);
std::cout << "------- GMRES solver (row scaling preconditioner) via ViennaCL, coordinate_matrix ----------" << std::endl;
run_solver(vcl_coordinate_matrix, vcl_vec2, vcl_result, gmres_solver, vcl_row_scaling_coo, gmres_ops);
return EXIT_SUCCESS;
}
int main()
{
std::cout << std::endl;
std::cout << "----------------------------------------------" << std::endl;
std::cout << " Device Info" << std::endl;
std::cout << "----------------------------------------------" << std::endl;
#ifdef VIENNACL_WITH_OPENCL
viennacl::ocl::platform pf = viennacl::ocl::get_platforms()[0];
std::vector<viennacl::ocl::device> const & devices = pf.devices();
// Set first device to first context:
viennacl::ocl::setup_context(0, devices[0]);
// Set second device for second context (use the same device for the second context if only one device available):
if (devices.size() > 1)
viennacl::ocl::setup_context(1, devices[1]);
else
viennacl::ocl::setup_context(1, devices[0]);
std::cout << viennacl::ocl::current_device().info() << std::endl;
viennacl::context ctx(viennacl::ocl::get_context(1));
#else
viennacl::context ctx;
#endif
std::cout << "---------------------------------------------------------------------------" << std::endl;
std::cout << "---------------------------------------------------------------------------" << std::endl;
std::cout << " Benchmark for Execution Times of Iterative Solvers provided with ViennaCL " << std::endl;
std::cout << "---------------------------------------------------------------------------" << std::endl;
std::cout << " Note that the purpose of this benchmark is not to run solvers until" << std::endl;
std::cout << " convergence. Instead, only the execution times of a few iterations are" << std::endl;
std::cout << " recorded. Residual errors are only printed for information." << std::endl << std::endl;
std::cout << std::endl;
std::cout << "----------------------------------------------" << std::endl;
std::cout << "----------------------------------------------" << std::endl;
std::cout << "## Benchmark :: Solver" << std::endl;
std::cout << "----------------------------------------------" << std::endl;
std::cout << std::endl;
std::cout << " -------------------------------" << std::endl;
std::cout << " # benchmarking single-precision" << std::endl;
std::cout << " -------------------------------" << std::endl;
run_benchmark<float>(ctx);
#ifdef VIENNACL_WITH_OPENCL
if ( viennacl::ocl::current_device().double_support() )
#endif
{
std::cout << std::endl;
std::cout << " -------------------------------" << std::endl;
std::cout << " # benchmarking double-precision" << std::endl;
std::cout << " -------------------------------" << std::endl;
run_benchmark<double>(ctx);
}
return 0;
}
viennacl::hyb_matrix
Sparse matrix class using a hybrid format composed of the ELL and CSR format for storing the nonzeros...
Definition: forwards.h:406
viennacl::tools::timer
Simple timer class based on gettimeofday (POSIX) or QueryPerformanceCounter (Windows).
Definition: timer.hpp:90
viennacl::linalg::chow_patel_tag
A tag for incomplete LU and incomplete Cholesky factorization with static pattern (Parallel-ILU0...
Definition: chow_patel_ilu.hpp:44
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::ocl::get_platforms
std::vector< platform > get_platforms()
Definition: platform.hpp:124
viennacl::linalg::ilu0_precond
ILU0 preconditioner class, can be supplied to solve()-routines.
Definition: ilu0.hpp:146
viennacl::linalg::ichol0_precond
Incomplete Cholesky preconditioner class with static pattern (ICHOL0), can be supplied to solve()-rou...
Definition: ichol.hpp:127
run_benchmark
int run_benchmark(viennacl::context ctx)
Definition: solver.cpp:127
viennacl::linalg::ichol0_tag
A tag for incomplete Cholesky factorization with static pattern (ILU0)
Definition: ichol.hpp:43
viennacl::linalg::jacobi_precond
Jacobi preconditioner class, can be supplied to solve()-routines. Generic version for non-ViennaCL ma...
Definition: jacobi_precond.hpp:50
viennacl::ocl::platform
Wrapper class for an OpenCL platform.
Definition: platform.hpp:45
viennacl::linalg::ilu0_tag
A tag for incomplete LU factorization with static pattern (ILU0)
Definition: ilu0.hpp:58
viennacl::ocl::platform::devices
std::vector< device > devices(cl_device_type dtype=CL_DEVICE_TYPE_DEFAULT)
Returns the available devices of the supplied device type.
Definition: platform.hpp:91
bicgstab.hpp
The stabilized bi-conjugate gradient method is implemented here.
diff_2
ScalarType diff_2(ublas::vector< ScalarType > &v1, viennacl::vector< ScalarType > &v2)
Definition: solver.cpp:91
ichol.hpp
Implementations of incomplete Cholesky factorization preconditioners with static nonzero pattern...
diff_inf
ScalarType diff_inf(ublas::vector< ScalarType > &v1, viennacl::vector< ScalarType > &v2)
Definition: solver.cpp:74
viennacl::backend::finish
void finish()
Synchronizes the execution. finish() will only return after all compute kernels (CUDA, OpenCL) have completed.
Definition: memory.hpp:54
viennacl::linalg::ilut_tag::use_level_scheduling
bool use_level_scheduling() const
Definition: ilut.hpp:76
viennacl::linalg::jacobi_tag
A tag for a jacobi preconditioner.
Definition: jacobi_precond.hpp:43
viennacl::linalg::chow_patel_icc_precond
Parallel Chow-Patel ILU preconditioner class, can be supplied to solve()-routines.
Definition: chow_patel_ilu.hpp:158
viennacl::linalg::block_ilu_precond
A block ILU preconditioner class, can be supplied to solve()-routines.
Definition: block_ilu.hpp:132
viennacl::linalg::detail::max
T max(const T &lhs, const T &rhs)
Maximum.
Definition: util.hpp:59
jacobi_precond.hpp
Implementation of a simple Jacobi preconditioner.
viennacl::linalg::row_scaling
Jacobi-type preconditioner class, can be supplied to solve()-routines. This is a diagonal preconditio...
Definition: row_scaling.hpp:87
viennacl::ocl::current_device
viennacl::ocl::device const & current_device()
Convenience function for returning the active device in the current context.
Definition: backend.hpp:351
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
coordinate_matrix.hpp
Implementation of the coordinate_matrix class.
viennacl::linalg::solve
VectorT solve(MatrixT const &matrix, VectorT const &rhs, bicgstab_tag const &tag, PreconditionerT const &precond)
Definition: bicgstab.hpp:496
context.hpp
Implementation of a OpenCL-like context, which serves as a unification of {OpenMP, CUDA, OpenCL} at the user API.
viennacl::ocl::device::info
std::string info(vcl_size_t indent=0, char indent_char= ' ') const
Returns an info string with a few properties of the device. Use full_info() to get all details...
Definition: device.hpp:995
viennacl::context
Represents a generic 'context' similar to an OpenCL context, but is backend-agnostic and thus also su...
Definition: context.hpp:39
viennacl::linalg::ilut_precond::apply
void apply(VectorT &vec) const
Definition: ilut.hpp:366
v1
viennacl::vector< float > v1
Definition: global_variables.cpp:60
hyb_matrix.hpp
Implementation of the hyb_matrix class.
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.
BENCHMARK_RUNS
#define BENCHMARK_RUNS
Definition: solver.cpp:64
viennacl::ell_matrix
Sparse matrix class using the ELLPACK format for storing the nonzeros.
Definition: ell_matrix.hpp:53
viennacl::vector_base::begin
iterator begin()
Returns an iterator pointing to the beginning of the vector (STL like)
Definition: vector.hpp:841
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::ilut_tag
A tag for incomplete LU factorization with threshold (ILUT)
Definition: ilut.hpp:45
viennacl::sliced_ell_matrix
Sparse matrix class using the sliced ELLPACK with parameters C, .
Definition: forwards.h:403
compressed_matrix.hpp
Implementation of the compressed_matrix class.
sliced_ell_matrix.hpp
Implementation of the sliced_ell_matrix class.
viennacl::linalg::row_scaling_tag
A tag for a row scaling preconditioner which merely normalizes the equation system such that each row...
Definition: row_scaling.hpp:40
run_solver
void run_solver(MatrixType const &matrix, VectorType const &rhs, VectorType const &ref_result, SolverTag const &solver, PrecondTag const &precond, long ops)
Definition: solver.cpp:101
viennacl::ocl::device::double_support
bool double_support() const
ViennaCL convenience function: Returns true if the device supports double precision.
Definition: device.hpp:956
viennacl::linalg::ilut_precond
ILUT preconditioner class, can be supplied to solve()-routines.
Definition: ilut.hpp:352
main
int main()
Definition: solver.cpp:625
cg.hpp
The conjugate gradient method is implemented here.
printOps
void printOps(double num_ops, double exec_time)
Definition: solver.cpp:67
viennacl::linalg::ilu0_precond::apply
void apply(VectorT &vec) const
Definition: ilu0.hpp:160
ell_matrix.hpp
Implementation of the ell_matrix class.
viennacl::linalg::block_ilu_precond::apply
void apply(VectorT &vec) const
Definition: block_ilu.hpp:174
viennacl::vector< ScalarType >
row_scaling.hpp
A row normalization preconditioner is implemented here.
timer.hpp
A simple, yet (mostly) sufficiently accurate timer for benchmarking and profiling.
viennacl::linalg::ichol0_precond::apply
void apply(VectorT &vec) const
Definition: ichol.hpp:141
v2
viennacl::vector< int > v2
Definition: global_variables.cpp:61
viennacl::linalg::ilu0_tag::use_level_scheduling
bool use_level_scheduling() const
Definition: ilu0.hpp:63
prod
void prod(std::vector< std::map< IndexT, NumericT > > const &stl_A, std::vector< std::map< IndexT, NumericT > > const &stl_B, std::vector< std::map< IndexT, NumericT > > &stl_C)
Definition: sparse_prod.cpp:114
vector.hpp
The vector type with operator-overloads and proxy classes is defined here. Linear algebra operations ...
viennacl::linalg::norm_inf
T norm_inf(std::vector< T, A > const &v1)
Definition: norm_inf.hpp:60
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::tools::timer::get
double get() const
Definition: timer.hpp:104
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
viennacl::vector_base::end
iterator end()
Returns an iterator pointing to the end of the vector (STL like)
Definition: vector.hpp:848
viennacl::tools::timer::start
void start()
Definition: timer.hpp:97
viennacl::ocl::get_context
viennacl::ocl::context & get_context(long i)
Convenience function for returning the current context.
Definition: backend.hpp:225
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
scalar.hpp
Implementation of the ViennaCL scalar class.
viennacl::ocl::setup_context
void setup_context(long i, std::vector< cl_device_id > const &devices)
Convenience function for setting devices for a context.
Definition: backend.hpp:231
mixed_precision_cg.hpp
The conjugate gradient method using mixed precision is implemented here. Experimental.
viennacl::linalg::chow_patel_ilu_precond
Parallel Chow-Patel ILU preconditioner class, can be supplied to solve()-routines.
Definition: chow_patel_ilu.hpp:236
viennacl::coordinate_matrix
A sparse square matrix, where entries are stored as triplets (i,j, val), where i and j are the row an...
Definition: coordinate_matrix.hpp:186