All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
amg.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 "viennacl/vector.hpp"
#include "viennacl/coordinate_matrix.hpp"
#include "viennacl/compressed_matrix.hpp"
#include "viennacl/linalg/ilu.hpp"
#include "viennacl/linalg/cg.hpp"
#include "viennacl/linalg/bicgstab.hpp"
#include "viennacl/io/matrix_market.hpp"
#include "viennacl/linalg/norm_2.hpp"
#include "viennacl/tools/matrix_generation.hpp"
#include "viennacl/linalg/amg.hpp"
#include <iostream>
#include <vector>
#include <ctime>
#include "vector-io.hpp"
#include "viennacl/tools/timer.hpp"
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)
{
VectorType result(rhs);
VectorType residual(rhs);
viennacl::tools::timer timer;
timer.start();
result = viennacl::linalg::solve(matrix, rhs, solver, precond);
viennacl::backend::finish();
std::cout << " > Solver time: " << timer.get() << std::endl;
residual -= viennacl::linalg::prod(matrix, result);
std::cout << " > Relative residual: " << viennacl::linalg::norm_2(residual) / viennacl::linalg::norm_2(rhs) << 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>
void run_amg(viennacl::linalg::cg_tag & cg_solver,
viennacl::vector<ScalarType> & vcl_vec,
viennacl::vector<ScalarType> & vcl_result,
viennacl::compressed_matrix<ScalarType> & vcl_compressed_matrix,
std::string info,
viennacl::linalg::amg_tag & amg_tag)
{
std::cout << "-- CG with AMG preconditioner, " << info << " --" << std::endl;
viennacl::linalg::amg_precond<viennacl::compressed_matrix<ScalarType> > vcl_amg(vcl_compressed_matrix, amg_tag);
std::cout << " * Setup phase (ViennaCL types)..." << std::endl;
viennacl::tools::timer timer;
timer.start();
vcl_amg.setup();
viennacl::backend::finish();
std::cout << " > Setup time: " << timer.get() << std::endl;
std::cout << " * CG solver (ViennaCL types)..." << std::endl;
run_solver(vcl_compressed_matrix, vcl_vec, vcl_result, cg_solver, vcl_amg);
}
int main(int argc, char **argv)
{
std::string filename("../examples/testdata/mat65k.mtx");
if (argc == 2)
filename = argv[1];
std::cout << std::endl;
std::cout << "----------------------------------------------" << std::endl;
std::cout << " Device Info" << std::endl;
std::cout << "----------------------------------------------" << std::endl;
#ifdef VIENNACL_WITH_OPENCL
// Optional: Customize OpenCL backend
viennacl::ocl::platform pf = viennacl::ocl::get_platforms()[0];
std::vector<viennacl::ocl::device> const & devices = pf.devices();
// Optional: Set first device to first context:
viennacl::ocl::setup_context(0, devices[0]);
// Optional: 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(0));
#else
viennacl::context ctx;
#endif
typedef double ScalarType; // feel free to change this to double if supported by your device
viennacl::compressed_matrix<ScalarType> vcl_compressed_matrix(ctx);
//viennacl::tools::generate_fdm_laplace(vcl_compressed_matrix, points_per_dim, points_per_dim);
// Read matrix
std::cout << "Reading matrix..." << std::endl;
std::vector< std::map<unsigned int, ScalarType> > read_in_matrix;
if (!viennacl::io::read_matrix_market_file(read_in_matrix, filename))
{
std::cout << "Error reading Matrix file" << std::endl;
return EXIT_FAILURE;
}
viennacl::copy(read_in_matrix, vcl_compressed_matrix);
std::cout << "Reading matrix completed." << std::endl;
viennacl::vector<ScalarType> vcl_vec(vcl_compressed_matrix.size1(), ctx);
viennacl::vector<ScalarType> vcl_result(vcl_compressed_matrix.size1(), ctx);
std::vector<ScalarType> std_vec, std_result;
// rhs and result vector:
std_vec.resize(vcl_compressed_matrix.size1());
std_result.resize(vcl_compressed_matrix.size1());
for (std::size_t i=0; i<std_result.size(); ++i)
std_result[i] = ScalarType(1);
// Copy to GPU
viennacl::copy(std_vec, vcl_vec);
viennacl::copy(std_result, vcl_result);
vcl_vec = viennacl::linalg::prod(vcl_compressed_matrix, vcl_result);
viennacl::linalg::cg_tag cg_solver(1e-8, 10000);
viennacl::context host_ctx(viennacl::MAIN_MEMORY);
viennacl::context target_ctx = viennacl::traits::context(vcl_compressed_matrix);
std::cout << "-- CG solver (no preconditioner, warmup) --" << std::endl;
run_solver(vcl_compressed_matrix, vcl_vec, vcl_result, cg_solver, viennacl::linalg::no_precond());
viennacl::linalg::amg_tag amg_tag_direct;
amg_tag_direct.set_coarsening_method(viennacl::linalg::AMG_COARSENING_METHOD_ONEPASS);
amg_tag_direct.set_interpolation_method(viennacl::linalg::AMG_INTERPOLATION_METHOD_DIRECT);
amg_tag_direct.set_strong_connection_threshold(0.25);
amg_tag_direct.set_jacobi_weight(0.67);
amg_tag_direct.set_presmooth_steps(1);
amg_tag_direct.set_postsmooth_steps(1);
amg_tag_direct.set_setup_context(host_ctx); // run setup on host
amg_tag_direct.set_target_context(target_ctx); // run solver cycles on device
run_amg(cg_solver, vcl_vec, vcl_result, vcl_compressed_matrix, "ONEPASS COARSENING, DIRECT INTERPOLATION", amg_tag_direct);
viennacl::linalg::amg_tag amg_tag_agg_pmis;
amg_tag_agg_pmis.set_coarsening_method(viennacl::linalg::AMG_COARSENING_METHOD_MIS2_AGGREGATION);
amg_tag_agg_pmis.set_interpolation_method(viennacl::linalg::AMG_INTERPOLATION_METHOD_AGGREGATION);
run_amg(cg_solver, vcl_vec, vcl_result, vcl_compressed_matrix, "AG COARSENING (PMIS), AG INTERPOLATION", amg_tag_agg_pmis);
viennacl::linalg::amg_tag amg_tag_sa_pmis;
amg_tag_sa_pmis.set_coarsening_method(viennacl::linalg::AMG_COARSENING_METHOD_MIS2_AGGREGATION);
amg_tag_sa_pmis.set_interpolation_method(viennacl::linalg::AMG_INTERPOLATION_METHOD_SMOOTHED_AGGREGATION);
run_amg (cg_solver, vcl_vec, vcl_result, vcl_compressed_matrix, "AG COARSENING (PMIS), SA INTERPOLATION", amg_tag_sa_pmis);
std::cout << std::endl;
std::cout << " -------------- Benchmark runs -------------- " << std::endl;
std::cout << std::endl;
std::cout << "-- CG solver (no preconditioner) --" << std::endl;
run_solver(vcl_compressed_matrix, vcl_vec, vcl_result, cg_solver, viennacl::linalg::no_precond());
run_amg(cg_solver, vcl_vec, vcl_result, vcl_compressed_matrix, "ONEPASS COARSENING, DIRECT INTERPOLATION", amg_tag_direct);
run_amg(cg_solver, vcl_vec, vcl_result, vcl_compressed_matrix, "AG COARSENING (PMIS), AG INTERPOLATION", amg_tag_agg_pmis);
run_amg (cg_solver, vcl_vec, vcl_result, vcl_compressed_matrix, "AG COARSENING (PMIS), SA INTERPOLATION", amg_tag_sa_pmis);
std::cout << "!!!! TUTORIAL COMPLETED SUCCESSFULLY !!!!" << std::endl;
return EXIT_SUCCESS;
}
viennacl::tools::timer
Simple timer class based on gettimeofday (POSIX) or QueryPerformanceCounter (Windows).
Definition: timer.hpp:90
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
norm_2.hpp
Generic interface for the l^2-norm. See viennacl/linalg/vector_operations.hpp for implementations...
viennacl::linalg::amg_tag::set_setup_context
void set_setup_context(viennacl::context ctx)
Sets the ViennaCL context for the setup stage. Set this to a host context if you want to run the setu...
Definition: amg_base.hpp:140
viennacl::ocl::platform
Wrapper class for an OpenCL platform.
Definition: platform.hpp:45
viennacl::linalg::amg_tag::set_strong_connection_threshold
void set_strong_connection_threshold(double threshold)
Sets the strong connection threshold. Customizations by the user essential for best results! ...
Definition: amg_base.hpp:100
viennacl::linalg::amg_tag::set_jacobi_weight
void set_jacobi_weight(double w)
Sets the weight (damping) for the Jacobi smoother.
Definition: amg_base.hpp:111
viennacl::compressed_matrix::size1
const vcl_size_t & size1() const
Returns the number of rows.
Definition: compressed_matrix.hpp:927
matrix_generation.hpp
Helper routines for generating sparse matrices.
main
int main()
Definition: bisect.cpp:91
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
viennacl::linalg::amg_tag::set_target_context
void set_target_context(viennacl::context ctx)
Sets the ViennaCL context for the solver cycle stage (i.e. preconditioner applications).
Definition: amg_base.hpp:148
bicgstab.hpp
The stabilized bi-conjugate gradient method is implemented here.
vector-io.hpp
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::AMG_INTERPOLATION_METHOD_AGGREGATION
Definition: amg_base.hpp:57
viennacl::linalg::amg_tag::set_postsmooth_steps
void set_postsmooth_steps(vcl_size_t steps)
Sets the number of smoother applications on the coarse level before interpolation to the finer level...
Definition: amg_base.hpp:121
viennacl::linalg::AMG_INTERPOLATION_METHOD_DIRECT
Definition: amg_base.hpp:56
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
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
viennacl::linalg::AMG_COARSENING_METHOD_MIS2_AGGREGATION
Definition: amg_base.hpp:50
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::prod
VectorT prod(std::vector< std::vector< T, A1 >, A2 > const &matrix, VectorT const &vector)
Definition: prod.hpp:102
viennacl::linalg::amg_tag::set_interpolation_method
void set_interpolation_method(amg_interpolation_method interpol)
Sets the interpolation method to the provided method.
Definition: amg_base.hpp:91
viennacl::linalg::AMG_INTERPOLATION_METHOD_SMOOTHED_AGGREGATION
Definition: amg_base.hpp:58
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::amg_tag::set_presmooth_steps
void set_presmooth_steps(vcl_size_t steps)
Sets the number of smoother applications on the fine level before restriction to the coarser level...
Definition: amg_base.hpp:116
compressed_matrix.hpp
Implementation of the compressed_matrix class.
viennacl::vector::resize
void resize(size_type new_size, bool preserve=true)
Resizes the allocated memory for the vector. Pads the memory to be a multiple of 'AlignmentV'.
Definition: vector.hpp:1046
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
amg.hpp
Main include file for algebraic multigrid (AMG) preconditioners. Experimental.
viennacl::linalg::amg_precond
AMG preconditioner class, can be supplied to solve()-routines.
Definition: amg.hpp:253
cg.hpp
The conjugate gradient method is implemented here.
viennacl::vector< ScalarType >
viennacl::linalg::amg_tag::set_coarsening_method
void set_coarsening_method(amg_coarsening_method s)
Sets the strategy used for constructing coarse grids.
Definition: amg_base.hpp:86
viennacl::MAIN_MEMORY
Definition: forwards.h:348
timer.hpp
A simple, yet (mostly) sufficiently accurate timer for benchmarking and profiling.
viennacl::linalg::amg_tag
A tag for algebraic multigrid (AMG). Used to transport information from the user to the implementatio...
Definition: amg_base.hpp:64
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
vector.hpp
The vector type with operator-overloads and proxy classes is defined here. Linear algebra operations ...
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< NumericT >::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::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::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::AMG_COARSENING_METHOD_ONEPASS
Definition: amg_base.hpp:48
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