1 #ifndef VIENNACL_LINALG_AMG_HPP_
2 #define VIENNACL_LINALG_AMG_HPP_
27 #include <boost/numeric/ublas/matrix.hpp>
28 #include <boost/numeric/ublas/lu.hpp>
29 #include <boost/numeric/ublas/operation.hpp>
30 #include <boost/numeric/ublas/vector_proxy.hpp>
31 #include <boost/numeric/ublas/matrix_proxy.hpp>
32 #include <boost/numeric/ublas/vector.hpp>
33 #include <boost/numeric/ublas/triangular.hpp>
47 #ifdef VIENNACL_WITH_OPENMP
53 #define VIENNACL_AMG_COARSE_LIMIT 50
54 #define VIENNACL_AMG_MAX_LEVELS 100
71 template <
typename InternalType1,
typename InternalType2>
72 void amg_setup(InternalType1 & A, InternalType1 & P, InternalType2 & Pointvector,
amg_tag & tag)
74 typedef typename InternalType2::value_type PointVectorType;
76 unsigned int i, iterations, c_points, f_points;
86 Slicing.
init(iterations);
88 for (i=0; i<iterations; ++i)
91 Pointvector[i] = PointVectorType(static_cast<unsigned int>(A[i].
size1()));
92 Pointvector[i].init_points();
98 c_points = Pointvector[i].get_cpoints();
99 f_points = Pointvector[i].get_fpoints();
101 #if defined (VIENNACL_AMG_DEBUG) //or defined(VIENNACL_AMG_DEBUGBENCH)
102 std::cout <<
"Level " << i <<
": ";
103 std::cout <<
"No of C points = " << c_points <<
", ";
104 std::cout <<
"No of F points = " << f_points << std::endl;
108 if (c_points == 0 || f_points == 0)
120 Pointvector[i].delete_points();
122 #ifdef VIENNACL_AMG_DEBUG
123 std::cout <<
"Coarse Grid Operator Matrix:" << std::endl;
145 template <
typename MatrixType,
typename InternalType1,
typename InternalType2>
146 void amg_init(MatrixType
const & mat, InternalType1 & A, InternalType1 & P, InternalType2 & Pointvector,
amg_tag & tag)
149 typedef typename InternalType1::value_type SparseMatrixType;
165 SparseMatrixType A0 (mat);
166 A.insert_element (0, A0);
178 template <
typename InternalType1,
typename InternalType2>
179 void amg_transform_cpu (InternalType1 & A, InternalType1 & P, InternalType1 & R, InternalType2 & A_setup, InternalType2 & P_setup,
amg_tag & tag)
191 A[i].resize(A_setup[i].
size1(),A_setup[i].
size2(),
false);
196 P[i].resize(P_setup[i].
size1(),P_setup[i].
size2(),
false);
201 R[i].resize(P_setup[i].
size2(),P_setup[i].
size1(),
false);
202 P_setup[i].set_trans(
true);
204 P_setup[i].set_trans(
false);
218 template <
typename InternalType1,
typename InternalType2>
244 P_setup[i].set_trans(
true);
246 P_setup[i].set_trans(
false);
258 template <
typename InternalVectorType,
typename SparseMatrixType>
259 void amg_setup_apply (InternalVectorType & result, InternalVectorType & rhs, InternalVectorType & residual, SparseMatrixType
const & A,
amg_tag const & tag)
261 typedef typename InternalVectorType::value_type VectorType;
269 result[level] = VectorType(A[level].
size1());
270 result[level].clear();
271 rhs[level] = VectorType(A[level].
size1());
276 residual[level] = VectorType(A[level].
size1());
277 residual[level].clear();
291 template <
typename InternalVectorType,
typename SparseMatrixType>
294 typedef typename InternalVectorType::value_type VectorType;
302 result[level] = VectorType(A[level].
size1(), ctx);
303 rhs[level] = VectorType(A[level].
size1(), ctx);
307 residual[level] = VectorType(A[level].
size1(), ctx);
319 template <
typename ScalarType,
typename SparseMatrixType>
320 void amg_lu(boost::numeric::ublas::compressed_matrix<ScalarType> & op, boost::numeric::ublas::permutation_matrix<> & Permutation, SparseMatrixType
const & A)
322 typedef typename SparseMatrixType::const_iterator1 ConstRowIterator;
323 typedef typename SparseMatrixType::const_iterator2 ConstColIterator;
326 op.resize(A.size1(),A.size2(),
false);
327 for (ConstRowIterator row_iter = A.begin1(); row_iter != A.end1(); ++row_iter)
328 for (ConstColIterator col_iter = row_iter.begin(); col_iter != row_iter.end(); ++col_iter)
329 op (col_iter.index1(), col_iter.index2()) = *col_iter;
332 Permutation = boost::numeric::ublas::permutation_matrix<> (op.size1());
338 template <
typename MatrixType>
341 typedef typename MatrixType::value_type ScalarType;
342 typedef boost::numeric::ublas::vector<ScalarType> VectorType;
351 boost::numeric::ublas::vector <SparseMatrixType> A_setup;
352 boost::numeric::ublas::vector <SparseMatrixType> P_setup;
353 boost::numeric::ublas::vector <MatrixType> A;
354 boost::numeric::ublas::vector <MatrixType> P;
355 boost::numeric::ublas::vector <MatrixType> R;
356 boost::numeric::ublas::vector <PointVectorType> Pointvector;
358 mutable boost::numeric::ublas::compressed_matrix<ScalarType> op;
359 mutable boost::numeric::ublas::permutation_matrix<> Permutation;
361 mutable boost::numeric::ublas::vector <VectorType> result;
362 mutable boost::numeric::ublas::vector <VectorType> rhs;
363 mutable boost::numeric::ublas::vector <VectorType> residual;
365 mutable bool done_init_apply;
380 amg_init (mat,A_setup,P_setup,Pointvector,tag_);
382 done_init_apply =
false;
390 amg_setup(A_setup,P_setup,Pointvector,tag_);
394 done_init_apply =
false;
408 done_init_apply =
true;
416 template <
typename VectorType>
420 unsigned int nonzero=0, systemmat_nonzero=0, level_coefficients=0;
424 level_coefficients = 0;
425 for (
InternalRowIterator row_iter = A_setup[level].begin1(); row_iter != A_setup[level].end1(); ++row_iter)
432 level_coefficients++;
435 avgstencil[level] = level_coefficients/
static_cast<ScalarType
>(A_setup[level].size1());
437 return nonzero/
static_cast<ScalarType
>(systemmat_nonzero);
444 template <
typename VectorType>
448 if (!done_init_apply)
457 result[level].clear();
462 #ifdef VIENNACL_AMG_DEBUG
463 std::cout <<
"After presmooth:" << std::endl;
470 #ifdef VIENNACL_AMG_DEBUG
471 std::cout <<
"Residual:" << std::endl;
478 #ifdef VIENNACL_AMG_DEBUG
479 std::cout <<
"Restricted Residual: " << std::endl;
485 result[level] = rhs[level];
488 #ifdef VIENNACL_AMG_DEBUG
489 std::cout <<
"After direct solve: " << std::endl;
495 #ifdef VIENNACL_AMG_DEBUG
496 std::cout <<
"Coarse Error: " << std::endl;
503 #ifdef VIENNACL_AMG_DEBUG
504 std::cout <<
"Corrected Result: " << std::endl;
511 #ifdef VIENNACL_AMG_DEBUG
512 std::cout <<
"After postsmooth: " << std::endl;
525 template <
typename VectorType>
526 void smooth_jacobi(
int level,
int const iterations, VectorType & x, VectorType
const & rhs)
const
528 VectorType old_result (x.size());
530 ScalarType sum = 0,
diag = 1;
532 for (
int i=0; i<iterations; ++i)
536 #ifdef VIENNACL_WITH_OPENMP
537 #pragma omp parallel for private (sum,diag) shared (rhs,x)
539 for (index=0; index < static_cast<long>(A_setup[level].size1()); ++index)
547 if (col_iter.index1() == col_iter.index2())
550 sum += *col_iter * old_result[col_iter.index2()];
564 template <
typename ScalarType,
unsigned int MAT_ALIGNMENT>
577 boost::numeric::ublas::vector <SparseMatrixType> A_setup;
578 boost::numeric::ublas::vector <SparseMatrixType> P_setup;
579 boost::numeric::ublas::vector <MatrixType> A;
580 boost::numeric::ublas::vector <MatrixType> P;
581 boost::numeric::ublas::vector <MatrixType> R;
582 boost::numeric::ublas::vector <PointVectorType> Pointvector;
584 mutable boost::numeric::ublas::compressed_matrix<ScalarType> op;
585 mutable boost::numeric::ublas::permutation_matrix<> Permutation;
587 mutable boost::numeric::ublas::vector <VectorType> result;
588 mutable boost::numeric::ublas::vector <VectorType> rhs;
589 mutable boost::numeric::ublas::vector <VectorType> residual;
593 mutable bool done_init_apply;
611 std::vector<std::map<unsigned int, ScalarType> > mat2 = std::vector<std::map<unsigned int, ScalarType> >(mat.
size1());
615 amg_init (mat2,A_setup,P_setup,Pointvector,tag_);
617 done_init_apply =
false;
625 amg_setup(A_setup,P_setup,Pointvector, tag_);
629 done_init_apply =
false;
643 done_init_apply =
true;
651 template <
typename VectorType>
655 unsigned int nonzero=0, systemmat_nonzero=0, level_coefficients=0;
659 level_coefficients = 0;
660 for (
InternalRowIterator row_iter = A_setup[level].begin1(); row_iter != A_setup[level].end1(); ++row_iter)
667 level_coefficients++;
670 avgstencil[level] = level_coefficients/(double)A[level].
size1();
672 return nonzero/
static_cast<double>(systemmat_nonzero);
679 template <
typename VectorType>
682 if (!done_init_apply)
691 result[level].clear();
696 #ifdef VIENNACL_AMG_DEBUG
697 std::cout <<
"After presmooth: " << std::endl;
704 residual[level] = rhs[level] - residual[level];
706 #ifdef VIENNACL_AMG_DEBUG
707 std::cout <<
"Residual: " << std::endl;
715 #ifdef VIENNACL_AMG_DEBUG
716 std::cout <<
"Restricted Residual: " << std::endl;
723 result[level] = rhs[level];
724 boost::numeric::ublas::vector <ScalarType> result_cpu (result[level].
size());
726 copy (result[level],result_cpu);
728 copy (result_cpu, result[level]);
730 #ifdef VIENNACL_AMG_DEBUG
731 std::cout <<
"After direct solve: " << std::endl;
737 #ifdef VIENNACL_AMG_DEBUG
738 std::cout <<
"Coarse Error: " << std::endl;
745 #ifdef VIENNACL_AMG_DEBUG
746 std::cout <<
"Corrected Result: " << std::endl;
753 #ifdef VIENNACL_AMG_DEBUG
754 std::cout <<
"After postsmooth: " << std::endl;
767 template <
typename VectorType>
776 for (
unsigned int i=0; i<iterations; ++i)
783 viennacl::traits::opencl_handle(old_result),
784 viennacl::traits::opencl_handle(x),
785 viennacl::traits::opencl_handle(rhs),
786 static_cast<cl_uint
>(rhs.
size())));
void amg_interpol(unsigned int level, InternalType1 &A, InternalType1 &P, InternalType2 &Pointvector, amg_tag &tag)
Calls the right function to build interpolation matrix.
Definition: amg_interpol.hpp:53
viennacl::ocl::kernel & get_kernel(std::string const &program_name, std::string const &kernel_name)
Convenience function for retrieving the kernel of a program directly from the context.
Definition: context.hpp:470
Debug functionality for AMG. To be removed.
double get_jacobiweight() const
Definition: amg_base.hpp:102
void set_coarselevels(int coarselevels)
Definition: amg_base.hpp:110
size_type size() const
Returns the length of the vector (cf. std::vector)
Definition: vector.hpp:859
amg_precond()
Definition: amg.hpp:599
void amg_lu(boost::numeric::ublas::compressed_matrix< ScalarType > &op, boost::numeric::ublas::permutation_matrix<> &Permutation, SparseMatrixType const &A)
Pre-compute LU factorization for direct solve (ublas library).
Definition: amg.hpp:320
Generic interface for matrix-vector and matrix-matrix products. See viennacl/linalg/vector_operations...
Represents an OpenCL kernel within ViennaCL.
Definition: kernel.hpp:59
amg_precond(compressed_matrix< ScalarType, MAT_ALIGNMENT > const &mat, amg_tag const &tag)
The constructor. Builds data structures.
Definition: amg.hpp:606
vcl_size_t size1(MatrixType const &mat)
Generic routine for obtaining the number of rows of a matrix (ViennaCL, uBLAS, etc.)
Definition: size.hpp:216
#define VIENNACL_AMG_COARSE_RS3
Definition: amg_base.hpp:44
Manages an OpenCL context and provides the respective convenience functions for creating buffers...
Definition: context.hpp:51
void prod(const T1 &A, bool transposed_A, const T2 &B, bool transposed_B, T3 &C, ScalarType alpha, ScalarType beta)
Definition: matrix_operations.hpp:2305
void init(unsigned int levels, unsigned int threads=0)
Definition: amg_base.hpp:1183
unsigned int get_coarse() const
Definition: amg_base.hpp:90
void amg_init(MatrixType const &mat, InternalType1 &A, InternalType1 &P, InternalType2 &Pointvector, amg_tag &tag)
Initialize AMG preconditioner.
Definition: amg.hpp:146
This file provides the forward declarations for the main types used within ViennaCL.
unsigned int get_postsmooth() const
Definition: amg_base.hpp:108
void printmatrix(MatrixType &, int)
Definition: amg_debug.hpp:77
void init_apply() const
Prepare data structures for preconditioning: Build data structures for precondition phase...
Definition: amg.hpp:401
void amg_transform_cpu(InternalType1 &A, InternalType1 &P, InternalType1 &R, InternalType2 &A_setup, InternalType2 &P_setup, amg_tag &tag)
Save operators after setup phase for CPU computation.
Definition: amg.hpp:179
#define VIENNACL_AMG_COARSE_LIMIT
Definition: amg.hpp:53
result_of::size_type< MatrixType >::type size2(MatrixType const &mat)
Generic routine for obtaining the number of columns of a matrix (ViennaCL, uBLAS, etc...
Definition: size.hpp:245
void apply(VectorType &vec) const
Precondition Operation.
Definition: amg.hpp:680
#define VIENNACL_AMG_COARSE_RS0
Definition: amg_base.hpp:43
amg_tag & tag()
Definition: amg.hpp:791
ScalarType calc_complexity(VectorType &avgstencil)
Returns complexity measures.
Definition: amg.hpp:417
void amg_setup(InternalType1 &A, InternalType1 &P, InternalType2 &Pointvector, amg_tag &tag)
Setup AMG preconditioner.
Definition: amg.hpp:72
Represents a generic 'context' similar to an OpenCL context, but is backend-agnostic and thus also su...
Definition: context.hpp:39
#define VIENNACL_AMG_MAX_LEVELS
Definition: amg.hpp:54
VectorT prod(std::vector< std::vector< T, A1 >, A2 > const &matrix, VectorT const &vector)
Definition: prod.hpp:91
void smooth_jacobi(int level, int const iterations, VectorType &x, VectorType const &rhs) const
(Weighted) Jacobi Smoother (CPU version)
Definition: amg.hpp:526
void clear()
Resets all entries to zero. Does not change the size of the vector.
Definition: vector.hpp:885
void lu_factorize(matrix< SCALARTYPE, viennacl::row_major > &A)
LU factorization of a row-major dense matrix.
Definition: lu.hpp:42
void setup()
Start setup phase for this class and copy data structures.
Definition: amg.hpp:622
vcl_size_t size(VectorType const &vec)
Generic routine for obtaining the size of a vector (ViennaCL, uBLAS, etc.)
Definition: size.hpp:144
A tag for algebraic multigrid (AMG). Used to transport information from the user to the implementatio...
Definition: amg_base.hpp:61
void enqueue(KernelType &k, viennacl::ocl::command_queue const &queue)
Enqueues a kernel in the provided queue.
Definition: enqueue.hpp:48
Main kernel class for generating OpenCL kernels for compressed_matrix.
Definition: compressed_matrix.hpp:1039
Implementations of several variants of the AMG coarsening procedure (setup phase). Experimental.
void amg_galerkin_prod(SparseMatrixType &A, SparseMatrixType &P, SparseMatrixType &RES)
Sparse Galerkin product: Calculates RES = trans(P)*A*P.
Definition: amg_base.hpp:1359
void copy(std::vector< SCALARTYPE > &cpu_vec, circulant_matrix< SCALARTYPE, ALIGNMENT > &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
unsigned int get_coarselevels() const
Definition: amg_base.hpp:111
static void init(viennacl::ocl::context &ctx)
Definition: compressed_matrix.hpp:1046
AMG preconditioner class, can be supplied to solve()-routines.
Definition: amg.hpp:339
const vcl_size_t & size1() const
Returns the number of rows.
Definition: compressed_matrix.hpp:692
void lu_substitute(matrix< SCALARTYPE, F1, ALIGNMENT_A > const &A, matrix< SCALARTYPE, F2, ALIGNMENT_B > &B)
LU substitution for the system LU = rhs.
Definition: lu.hpp:201
amg_precond(MatrixType const &mat, amg_tag const &tag)
The constructor. Saves system matrix, tag and builds data structures for setup.
Definition: amg.hpp:376
A class for the sparse matrix type. Uses vector of maps as data structure for higher performance and ...
Definition: amg_base.hpp:376
Implementations of dense direct solvers are found here.
void amg_setup_apply(InternalVectorType &result, InternalVectorType &rhs, InternalVectorType &residual, SparseMatrixType const &A, amg_tag const &tag)
Setup data structures for precondition phase.
Definition: amg.hpp:259
amg_precond()
Definition: amg.hpp:370
void printvector(VectorType const &)
Definition: amg_debug.hpp:80
vector_expression< const matrix_base< NumericT, F >, const int, op_matrix_diag > diag(const matrix_base< NumericT, F > &A, int k=0)
Definition: matrix.hpp:895
void amg_transform_gpu(InternalType1 &A, InternalType1 &P, InternalType1 &R, InternalType2 &A_setup, InternalType2 &P_setup, amg_tag &tag, viennacl::context ctx)
Save operators after setup phase for GPU computation.
Definition: amg.hpp:219
void init_apply() const
Prepare data structures for preconditioning: Build data structures for precondition phase...
Definition: amg.hpp:636
void apply(VectorType &vec) const
Precondition Operation.
Definition: amg.hpp:445
amg_tag & tag()
Definition: amg.hpp:557
void smooth_jacobi(int level, unsigned int iterations, VectorType &x, VectorType const &rhs) const
Jacobi Smoother (GPU version)
Definition: amg.hpp:768
A sparse square matrix in compressed sparse rows format.
Definition: compressed_matrix.hpp:428
Helper classes and functions for the AMG preconditioner. Experimental.
void amg_coarse(unsigned int level, InternalType1 &A, InternalType2 &Pointvector, InternalType3 &Slicing, amg_tag &tag)
Calls the right coarsening procedure.
Definition: amg_coarse.hpp:52
viennacl::backend::mem_handle & handle(T &obj)
Returns the generic memory handle of an object. Non-const version.
Definition: handle.hpp:41
ScalarType calc_complexity(VectorType &avgstencil)
Returns complexity measures.
Definition: amg.hpp:652
void setup()
Start setup phase for this class and copy data structures.
Definition: amg.hpp:387
A class for the matrix slicing for parallel coarsening schemes (RS0/RS3).
Definition: amg_base.hpp:1168
unsigned int get_presmooth() const
Definition: amg_base.hpp:105
A class for the AMG points. Holds pointers of type amg_point in a vector that can be accessed using [...
Definition: amg_base.hpp:935
Implementations of several variants of the AMG interpolation operators (setup phase). Experimental.
detail::amg::amg_tag amg_tag
Definition: amg.hpp:60
void switch_memory_context(T &obj, viennacl::context new_ctx)
Generic convenience routine for migrating data of an object to a new memory domain.
Definition: memory.hpp:624