ViennaCL - The Vienna Computing Library  1.7.0
Free open-source GPU-accelerated linear algebra and solver library.
amg.cpp
Go to the documentation of this file.
1 /* =========================================================================
2  Copyright (c) 2010-2015, Institute for Microelectronics,
3  Institute for Analysis and Scientific Computing,
4  TU Wien.
5  Portions of this software are copyright by UChicago Argonne, LLC.
6 
7  -----------------
8  ViennaCL - The Vienna Computing Library
9  -----------------
10 
11  Project Head: Karl Rupp rupp@iue.tuwien.ac.at
12 
13  (A list of authors and contributors can be found in the PDF manual)
14 
15  License: MIT (X11), see file LICENSE in the base directory
16 ============================================================================= */
17 
18 
27 #include "viennacl/vector.hpp"
30 #include "viennacl/linalg/ilu.hpp"
31 #include "viennacl/linalg/cg.hpp"
36 
40 #include "viennacl/linalg/amg.hpp"
41 
45 #include <iostream>
46 #include <vector>
47 #include <ctime>
48 #include "vector-io.hpp"
49 #include "viennacl/tools/timer.hpp"
50 
51 
57 template<typename MatrixType, typename VectorType, typename SolverTag, typename PrecondTag>
58 void run_solver(MatrixType const & matrix, VectorType const & rhs, VectorType const & ref_result, SolverTag const & solver, PrecondTag const & precond)
59 {
60  VectorType result(rhs);
61  VectorType residual(rhs);
62 
64  timer.start();
65  result = viennacl::linalg::solve(matrix, rhs, solver, precond);
67  std::cout << " > Solver time: " << timer.get() << std::endl;
68  residual -= viennacl::linalg::prod(matrix, result);
69  std::cout << " > Relative residual: " << viennacl::linalg::norm_2(residual) / viennacl::linalg::norm_2(rhs) << std::endl;
70  std::cout << " > Iterations: " << solver.iters() << std::endl;
71  result -= ref_result;
72  std::cout << " > Relative deviation from result: " << viennacl::linalg::norm_2(result) / viennacl::linalg::norm_2(ref_result) << std::endl;
73 }
74 
80 template<typename ScalarType>
81 void run_amg(viennacl::linalg::cg_tag & cg_solver,
83  viennacl::vector<ScalarType> & vcl_result,
84  viennacl::compressed_matrix<ScalarType> & vcl_compressed_matrix,
85  std::string info,
86  viennacl::linalg::amg_tag & amg_tag)
87 {
88  std::cout << "-- CG with AMG preconditioner, " << info << " --" << std::endl;
89 
90  viennacl::linalg::amg_precond<viennacl::compressed_matrix<ScalarType> > vcl_amg(vcl_compressed_matrix, amg_tag);
91  std::cout << " * Setup phase (ViennaCL types)..." << std::endl;
93  timer.start();
94  vcl_amg.setup();
96  std::cout << " > Setup time: " << timer.get() << std::endl;
97 
98  std::cout << " * CG solver (ViennaCL types)..." << std::endl;
99  run_solver(vcl_compressed_matrix, vcl_vec, vcl_result, cg_solver, vcl_amg);
100 }
101 
107 int main(int argc, char **argv)
108 {
109  std::string filename("../examples/testdata/mat65k.mtx");
110  if (argc == 2)
111  filename = argv[1];
112 
116  std::cout << std::endl;
117  std::cout << "----------------------------------------------" << std::endl;
118  std::cout << " Device Info" << std::endl;
119  std::cout << "----------------------------------------------" << std::endl;
120 
121 #ifdef VIENNACL_WITH_OPENCL
122  // Optional: Customize OpenCL backend
124  std::vector<viennacl::ocl::device> const & devices = pf.devices();
125 
126  // Optional: Set first device to first context:
127  viennacl::ocl::setup_context(0, devices[0]);
128 
129  // Optional: Set second device for second context (use the same device for the second context if only one device available):
130  if (devices.size() > 1)
131  viennacl::ocl::setup_context(1, devices[1]);
132  else
133  viennacl::ocl::setup_context(1, devices[0]);
134 
135  std::cout << viennacl::ocl::current_device().info() << std::endl;
137 #else
138  viennacl::context ctx;
139 #endif
140 
141  typedef double ScalarType; // feel free to change this to double if supported by your device
142 
143 
147  viennacl::compressed_matrix<ScalarType> vcl_compressed_matrix(ctx);
148 
149  //viennacl::tools::generate_fdm_laplace(vcl_compressed_matrix, points_per_dim, points_per_dim);
150  // Read matrix
151  std::cout << "Reading matrix..." << std::endl;
152  std::vector< std::map<unsigned int, ScalarType> > read_in_matrix;
153  if (!viennacl::io::read_matrix_market_file(read_in_matrix, filename))
154  {
155  std::cout << "Error reading Matrix file" << std::endl;
156  return EXIT_FAILURE;
157  }
158  viennacl::copy(read_in_matrix, vcl_compressed_matrix);
159  std::cout << "Reading matrix completed." << std::endl;
160 
161  viennacl::vector<ScalarType> vcl_vec(vcl_compressed_matrix.size1(), ctx);
162  viennacl::vector<ScalarType> vcl_result(vcl_compressed_matrix.size1(), ctx);
163 
164  std::vector<ScalarType> std_vec, std_result;
165 
166 
167  // rhs and result vector:
168  std_vec.resize(vcl_compressed_matrix.size1());
169  std_result.resize(vcl_compressed_matrix.size1());
170  for (std::size_t i=0; i<std_result.size(); ++i)
171  std_result[i] = ScalarType(1);
172 
173  // Copy to GPU
174  viennacl::copy(std_vec, vcl_vec);
175  viennacl::copy(std_result, vcl_result);
176 
177  vcl_vec = viennacl::linalg::prod(vcl_compressed_matrix, vcl_result);
178 
179 
183  viennacl::linalg::cg_tag cg_solver(1e-8, 10000);
184 
186  viennacl::context target_ctx = viennacl::traits::context(vcl_compressed_matrix);
187 
192  std::cout << "-- CG solver (no preconditioner, warmup) --" << std::endl;
193  run_solver(vcl_compressed_matrix, vcl_vec, vcl_result, cg_solver, viennacl::linalg::no_precond());
194 
198  viennacl::linalg::amg_tag amg_tag_direct;
201  amg_tag_direct.set_strong_connection_threshold(0.25);
202  amg_tag_direct.set_jacobi_weight(0.67);
203  amg_tag_direct.set_presmooth_steps(1);
204  amg_tag_direct.set_postsmooth_steps(1);
205  amg_tag_direct.set_setup_context(host_ctx); // run setup on host
206  amg_tag_direct.set_target_context(target_ctx); // run solver cycles on device
207  run_amg(cg_solver, vcl_vec, vcl_result, vcl_compressed_matrix, "ONEPASS COARSENING, DIRECT INTERPOLATION", amg_tag_direct);
208 
212  viennacl::linalg::amg_tag amg_tag_agg_pmis;
215  run_amg(cg_solver, vcl_vec, vcl_result, vcl_compressed_matrix, "AG COARSENING (PMIS), AG INTERPOLATION", amg_tag_agg_pmis);
216 
220  viennacl::linalg::amg_tag amg_tag_sa_pmis;
223  run_amg (cg_solver, vcl_vec, vcl_result, vcl_compressed_matrix, "AG COARSENING (PMIS), SA INTERPOLATION", amg_tag_sa_pmis);
224 
225  std::cout << std::endl;
226  std::cout << " -------------- Benchmark runs -------------- " << std::endl;
227  std::cout << std::endl;
228 
229  std::cout << "-- CG solver (no preconditioner) --" << std::endl;
230  run_solver(vcl_compressed_matrix, vcl_vec, vcl_result, cg_solver, viennacl::linalg::no_precond());
231  run_amg(cg_solver, vcl_vec, vcl_result, vcl_compressed_matrix, "ONEPASS COARSENING, DIRECT INTERPOLATION", amg_tag_direct);
232  run_amg(cg_solver, vcl_vec, vcl_result, vcl_compressed_matrix, "AG COARSENING (PMIS), AG INTERPOLATION", amg_tag_agg_pmis);
233  run_amg (cg_solver, vcl_vec, vcl_result, vcl_compressed_matrix, "AG COARSENING (PMIS), SA INTERPOLATION", amg_tag_sa_pmis);
234 
238  std::cout << "!!!! TUTORIAL COMPLETED SUCCESSFULLY !!!!" << std::endl;
239 
240  return EXIT_SUCCESS;
241 }
242 
Simple timer class based on gettimeofday (POSIX) or QueryPerformanceCounter (Windows).
Definition: timer.hpp:90
T norm_2(std::vector< T, A > const &v1)
Definition: norm_2.hpp:96
A reader and writer for the matrix market format is implemented here.
std::vector< platform > get_platforms()
Definition: platform.hpp:124
Generic interface for the l^2-norm. See viennacl/linalg/vector_operations.hpp for implementations...
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
Wrapper class for an OpenCL platform.
Definition: platform.hpp:45
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
void set_jacobi_weight(double w)
Sets the weight (damping) for the Jacobi smoother.
Definition: amg_base.hpp:111
const vcl_size_t & size1() const
Returns the number of rows.
Helper routines for generating sparse matrices.
int main()
Definition: bisect.cpp:91
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
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
The stabilized bi-conjugate gradient method is implemented here.
void finish()
Synchronizes the execution. finish() will only return after all compute kernels (CUDA, OpenCL) have completed.
Definition: memory.hpp:54
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::ocl::device const & current_device()
Convenience function for returning the active device in the current context.
Definition: backend.hpp:351
Implementation of the coordinate_matrix class.
VectorT solve(MatrixT const &matrix, VectorT const &rhs, bicgstab_tag const &tag, PreconditionerT const &precond)
Definition: bicgstab.hpp:496
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
Represents a generic 'context' similar to an OpenCL context, but is backend-agnostic and thus also su...
Definition: context.hpp:39
VectorT prod(std::vector< std::vector< T, A1 >, A2 > const &matrix, VectorT const &vector)
Definition: prod.hpp:102
void set_interpolation_method(amg_interpolation_method interpol)
Sets the interpolation method to the provided method.
Definition: amg_base.hpp:91
A tag class representing the use of no preconditioner.
Definition: forwards.h:873
Implementations of incomplete factorization preconditioners. Convenience header file.
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
Implementation of the compressed_matrix class.
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
Main include file for algebraic multigrid (AMG) preconditioners. Experimental.
AMG preconditioner class, can be supplied to solve()-routines.
Definition: amg.hpp:253
The conjugate gradient method is implemented here.
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
void set_coarsening_method(amg_coarsening_method s)
Sets the strategy used for constructing coarse grids.
Definition: amg_base.hpp:86
A simple, yet (mostly) sufficiently accurate timer for benchmarking and profiling.
A tag for algebraic multigrid (AMG). Used to transport information from the user to the implementatio...
Definition: amg_base.hpp:64
viennacl::context context(T const &t)
Returns an ID for the currently active memory domain of an object.
Definition: context.hpp:40
The vector type with operator-overloads and proxy classes is defined here. Linear algebra operations ...
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) ...
size_type size() const
Returns the length of the vector (cf. std::vector)
Definition: vector_def.hpp:118
float ScalarType
Definition: fft_1d.cpp:42
A tag for the conjugate gradient Used for supplying solver parameters and for dispatching the solve()...
Definition: cg.hpp:48
A sparse square matrix in compressed sparse rows format.
double get() const
Definition: timer.hpp:104
long read_matrix_market_file(MatrixT &mat, const char *file, long index_base=1)
Reads a sparse matrix from a file (MatrixMarket format)
viennacl::ocl::context & get_context(long i)
Convenience function for returning the current context.
Definition: backend.hpp:225
void setup_context(long i, std::vector< cl_device_id > const &devices)
Convenience function for setting devices for a context.
Definition: backend.hpp:231