1 #ifndef VIENNACL_LINALG_AMG_HPP_
2 #define VIENNACL_LINALG_AMG_HPP_
26 #include <boost/numeric/ublas/matrix.hpp>
27 #include <boost/numeric/ublas/lu.hpp>
28 #include <boost/numeric/ublas/operation.hpp>
29 #include <boost/numeric/ublas/vector_proxy.hpp>
30 #include <boost/numeric/ublas/matrix_proxy.hpp>
31 #include <boost/numeric/ublas/vector.hpp>
32 #include <boost/numeric/ublas/triangular.hpp>
52 #define VIENNACL_AMG_COARSE_LIMIT 50
53 #define VIENNACL_AMG_MAX_LEVELS 100
70 template <
typename InternalType1,
typename InternalType2>
71 void amg_setup(InternalType1 & A, InternalType1 & P, InternalType2 & Pointvector,
amg_tag & tag)
73 typedef typename InternalType1::value_type SparseMatrixType;
74 typedef typename InternalType2::value_type PointVectorType;
75 typedef typename SparseMatrixType::value_type ScalarType;
76 typedef typename SparseMatrixType::iterator1 InternalRowIterator;
77 typedef typename SparseMatrixType::iterator2 InternalColIterator;
79 unsigned int i, iterations, c_points, f_points;
89 Slicing.
init(iterations);
91 for (i=0; i<iterations; ++i)
94 Pointvector[i] = PointVectorType(A[i].
size1());
95 Pointvector[i].init_points();
101 c_points = Pointvector[i].get_cpoints();
102 f_points = Pointvector[i].get_fpoints();
104 #if defined (DEBUG) //or defined(DEBUGBENCH)
105 std::cout <<
"Level " << i <<
": ";
106 std::cout <<
"No of C points = " << c_points <<
", ";
107 std::cout <<
"No of F points = " << f_points << std::endl;
111 if (c_points == 0 || f_points == 0)
123 Pointvector[i].delete_points();
126 std::cout <<
"Coarse Grid Operator Matrix:" << std::endl;
148 template <
typename MatrixType,
typename InternalType1,
typename InternalType2>
149 void amg_init(MatrixType
const & mat, InternalType1 & A, InternalType1 & P, InternalType2 & Pointvector,
amg_tag & tag)
151 typedef typename MatrixType::value_type ScalarType;
152 typedef typename InternalType1::value_type SparseMatrixType;
168 SparseMatrixType A0 (mat);
169 A.insert_element (0, A0);
181 template <
typename InternalType1,
typename InternalType2>
182 void amg_transform_cpu (InternalType1 & A, InternalType1 & P, InternalType1 & R, InternalType2 & A_setup, InternalType2 & P_setup,
amg_tag & tag)
184 typedef typename InternalType1::value_type MatrixType;
194 A[i].resize(A_setup[i].
size1(),A_setup[i].
size2(),
false);
199 P[i].resize(P_setup[i].
size1(),P_setup[i].
size2(),
false);
204 R[i].resize(P_setup[i].
size2(),P_setup[i].
size1(),
false);
205 P_setup[i].set_trans(
true);
207 P_setup[i].set_trans(
false);
220 template <
typename InternalType1,
typename InternalType2>
221 void amg_transform_gpu (InternalType1 & A, InternalType1 & P, InternalType1 & R, InternalType2 & A_setup, InternalType2 & P_setup,
amg_tag & tag)
223 typedef typename InternalType1::value_type MatrixType;
224 typedef typename InternalType2::value_type::value_type ScalarType;
234 A[i].resize(A_setup[i].
size1(),A_setup[i].
size2(),
false);
239 P[i].resize(P_setup[i].
size1(),P_setup[i].
size2(),
false);
245 R[i].resize(P_setup[i].
size2(),P_setup[i].
size1(),
false);
246 P_setup[i].set_trans(
true);
248 P_setup[i].set_trans(
false);
260 template <
typename InternalVectorType,
typename SparseMatrixType>
261 void amg_setup_apply (InternalVectorType & result, InternalVectorType & rhs, InternalVectorType & residual, SparseMatrixType
const & A,
amg_tag const & tag)
263 typedef typename InternalVectorType::value_type VectorType;
271 result[level] = VectorType(A[level].
size1());
272 result[level].clear();
273 rhs[level] = VectorType(A[level].
size1());
278 residual[level] = VectorType(A[level].
size1());
279 residual[level].clear();
289 template <
typename ScalarType,
typename SparseMatrixType>
290 void amg_lu(boost::numeric::ublas::compressed_matrix<ScalarType> & op, boost::numeric::ublas::permutation_matrix<ScalarType> & Permutation, SparseMatrixType
const & A)
292 typedef typename SparseMatrixType::const_iterator1 ConstRowIterator;
293 typedef typename SparseMatrixType::const_iterator2 ConstColIterator;
296 op.resize(A.size1(),A.size2(),
false);
297 for (ConstRowIterator row_iter = A.begin1(); row_iter != A.end1(); ++row_iter)
298 for (ConstColIterator col_iter = row_iter.begin(); col_iter != row_iter.end(); ++col_iter)
299 op (col_iter.index1(), col_iter.index2()) = *col_iter;
302 Permutation = boost::numeric::ublas::permutation_matrix<ScalarType> (op.size1());
308 template <
typename MatrixType>
311 typedef typename MatrixType::value_type ScalarType;
312 typedef boost::numeric::ublas::vector<ScalarType> VectorType;
321 boost::numeric::ublas::vector <SparseMatrixType> A_setup;
322 boost::numeric::ublas::vector <SparseMatrixType> P_setup;
323 boost::numeric::ublas::vector <MatrixType> A;
324 boost::numeric::ublas::vector <MatrixType> P;
325 boost::numeric::ublas::vector <MatrixType> R;
326 boost::numeric::ublas::vector <PointVectorType> Pointvector;
328 mutable boost::numeric::ublas::compressed_matrix<ScalarType> op;
329 mutable boost::numeric::ublas::permutation_matrix<ScalarType> Permutation;
331 mutable boost::numeric::ublas::vector <VectorType> result;
332 mutable boost::numeric::ublas::vector <VectorType> rhs;
333 mutable boost::numeric::ublas::vector <VectorType> residual;
335 mutable bool done_init_apply;
350 amg_init (mat,A_setup,P_setup,Pointvector,_tag);
352 done_init_apply =
false;
360 amg_setup(A_setup,P_setup,Pointvector,_tag);
364 done_init_apply =
false;
378 done_init_apply =
true;
386 template <
typename VectorType>
390 unsigned int nonzero=0, systemmat_nonzero=0, level_coefficients=0;
394 level_coefficients = 0;
395 for (
InternalRowIterator row_iter = A_setup[level].begin1(); row_iter != A_setup[level].end1(); ++row_iter)
402 level_coefficients++;
405 avgstencil[level] = level_coefficients/(double)A_setup[level].
size1();
407 return nonzero/
static_cast<double>(systemmat_nonzero);
414 template <
typename VectorType>
418 if (!done_init_apply)
427 result[level].clear();
433 std::cout <<
"After presmooth:" << std::endl;
441 std::cout <<
"Residual:" << std::endl;
449 std::cout <<
"Restricted Residual: " << std::endl;
455 result[level] = rhs[level];
459 std::cout <<
"After direct solve: " << std::endl;
466 std::cout <<
"Coarse Error: " << std::endl;
474 std::cout <<
"Corrected Result: " << std::endl;
482 std::cout <<
"After postsmooth: " << std::endl;
495 template <
typename VectorType>
496 void smooth_jacobi(
int level,
int const iterations, VectorType & x, VectorType
const & rhs)
const
498 VectorType old_result (x.size());
500 ScalarType sum = 0, diag = 1;
502 for (
int i=0; i<iterations; ++i)
507 #pragma omp parallel for private (sum,diag) shared (rhs,x)
509 for (index=0; index<A_setup[level].size1(); ++index)
517 if (col_iter.index1() == col_iter.index2())
520 sum += *col_iter * old_result[col_iter.index2()];
534 template <
typename ScalarType,
unsigned int MAT_ALIGNMENT>
547 boost::numeric::ublas::vector <SparseMatrixType> A_setup;
548 boost::numeric::ublas::vector <SparseMatrixType> P_setup;
549 boost::numeric::ublas::vector <MatrixType> A;
550 boost::numeric::ublas::vector <MatrixType> P;
551 boost::numeric::ublas::vector <MatrixType> R;
552 boost::numeric::ublas::vector <PointVectorType> Pointvector;
554 mutable boost::numeric::ublas::compressed_matrix<ScalarType> op;
555 mutable boost::numeric::ublas::permutation_matrix<ScalarType> Permutation;
557 mutable boost::numeric::ublas::vector <VectorType> result;
558 mutable boost::numeric::ublas::vector <VectorType> rhs;
559 mutable boost::numeric::ublas::vector <VectorType> residual;
561 mutable bool done_init_apply;
579 std::vector<std::map<unsigned int, ScalarType> > mat2 = std::vector<std::map<unsigned int, ScalarType> >(mat.
size1());
583 amg_init (mat2,A_setup,P_setup,Pointvector,_tag);
585 done_init_apply =
false;
593 amg_setup(A_setup,P_setup,Pointvector,_tag);
597 done_init_apply =
false;
611 done_init_apply =
true;
619 template <
typename VectorType>
623 unsigned int nonzero=0, systemmat_nonzero=0, level_coefficients=0;
627 level_coefficients = 0;
628 for (
InternalRowIterator row_iter = A_setup[level].begin1(); row_iter != A_setup[level].end1(); ++row_iter)
635 level_coefficients++;
638 avgstencil[level] = level_coefficients/(double)A[level].
size1();
640 return nonzero/
static_cast<double>(systemmat_nonzero);
647 template <
typename VectorType>
650 if (!done_init_apply)
659 result[level].clear();
665 std::cout <<
"After presmooth: " << std::endl;
673 std::cout <<
"Residual: " << std::endl;
682 std::cout <<
"Restricted Residual: " << std::endl;
689 result[level] = rhs[level];
690 boost::numeric::ublas::vector <ScalarType> result_cpu (result[level].
size());
692 copy (result[level],result_cpu);
694 copy (result_cpu, result[level]);
697 std::cout <<
"After direct solve: " << std::endl;
704 std::cout <<
"Coarse Error: " << std::endl;
712 std::cout <<
"Corrected Result: " << std::endl;
720 std::cout <<
"After postsmooth: " << std::endl;
733 template <
typename VectorType>
745 for (
unsigned int i=0; i<iterations; ++i)
754 static_cast<cl_uint>(rhs.
size())));