1 #ifndef VIENNACL_LINALG_DETAIL_AMG_AMG_INTERPOL_HPP
2 #define VIENNACL_LINALG_DETAIL_AMG_AMG_INTERPOL_HPP
24 #include <boost/numeric/ublas/vector.hpp>
51 template <
typename InternalType1,
typename InternalType2>
52 void amg_interpol(
unsigned int level, InternalType1 & A, InternalType1 & P, InternalType2 & Pointvector,
amg_tag & tag)
69 template <
typename InternalType1,
typename InternalType2>
72 typedef typename InternalType1::value_type SparseMatrixType;
73 typedef typename InternalType2::value_type PointVectorType;
74 typedef typename SparseMatrixType::value_type ScalarType;
75 typedef typename SparseMatrixType::iterator1 InternalRowIterator;
76 typedef typename SparseMatrixType::iterator2 InternalColIterator;
79 ScalarType row_sum, c_sum, diag;
83 unsigned int c_points = Pointvector[level].get_cpoints();
86 P[level] = SparseMatrixType(A[level].
size1(),c_points);
90 Pointvector[level].build_index();
94 #pragma omp parallel for private (pointx,pointy,row_sum,c_sum,temp_res,y,x,diag) shared (P,A,Pointvector,tag)
96 for (x=0; x < Pointvector[level].size(); ++x)
98 pointx = Pointvector[level][x];
112 InternalRowIterator row_iter = A[level].begin1();
116 row_sum = c_sum = diag = 0;
117 for (InternalColIterator col_iter = row_iter.begin(); col_iter != row_iter.end(); ++col_iter)
119 y = col_iter.index2();
127 row_sum += *col_iter;
129 pointy = Pointvector[level][y];
135 temp_res = -row_sum/(c_sum*diag);
159 std::cout <<
"Prolongation Matrix:" << std::endl;
170 template <
typename InternalType1,
typename InternalType2>
173 typedef typename InternalType1::value_type SparseMatrixType;
174 typedef typename InternalType2::value_type PointVectorType;
175 typedef typename SparseMatrixType::value_type ScalarType;
176 typedef typename SparseMatrixType::iterator1 InternalRowIterator;
177 typedef typename SparseMatrixType::iterator2 InternalColIterator;
180 ScalarType weak_sum, strong_sum;
183 amg_point *pointx, *pointy, *pointk, *pointm;
184 unsigned int x, y, k, m;
186 unsigned int c_points = Pointvector[level].get_cpoints();
189 P[level] = SparseMatrixType(A[level].
size1(), c_points);
193 Pointvector[level].build_index();
197 #pragma omp parallel for private (pointx,pointy,pointk,pointm,weak_sum,strong_sum,c_sum_row,temp_res,x,y,k,m,diag_sign) shared (A,P,Pointvector)
199 for (x=0; x < Pointvector[level].size(); ++x)
201 pointx = Pointvector[level][x];
202 if (A[level](x,x) > 0)
215 InternalRowIterator row_iter = A[level].begin1();
221 for (InternalColIterator col_iter = row_iter.begin(); col_iter != row_iter.end(); ++col_iter)
223 k = col_iter.index2();
224 pointk = Pointvector[level][k];
229 weak_sum += *col_iter;
243 if (A[level](k,m) * diag_sign < 0)
244 c_sum_row[k] += A[level](k,m);
265 if (A[level](k,y) * diag_sign < 0)
266 strong_sum += (A[level](x,k) * A[level](k,y)) / (*iter2);
270 temp_res = - (A[level](x,y) + strong_sum) / (weak_sum);
283 std::cout <<
"Prolongation Matrix:" << std::endl;
294 template <
typename SparseMatrixType>
297 typedef typename SparseMatrixType::value_type ScalarType;
298 typedef typename SparseMatrixType::iterator1 InternalRowIterator;
299 typedef typename SparseMatrixType::iterator2 InternalColIterator;
301 ScalarType row_max, row_min, row_sum_pos, row_sum_neg, row_sum_pos_scale, row_sum_neg_scale;
303 InternalRowIterator row_iter = P.begin1();
313 for (InternalColIterator col_iter = row_iter.begin(); col_iter != row_iter.end(); ++col_iter)
315 if (*col_iter > row_max)
317 if (*col_iter < row_min)
320 row_sum_pos += *col_iter;
322 row_sum_neg += *col_iter;
325 row_sum_pos_scale = row_sum_pos;
326 row_sum_neg_scale = row_sum_neg;
329 for (InternalColIterator col_iter = row_iter.begin(); col_iter != row_iter.end(); ++col_iter)
333 row_sum_pos_scale -= *col_iter;
338 row_sum_pos_scale -= *col_iter;
344 for (InternalColIterator col_iter = row_iter.begin(); col_iter != row_iter.end(); ++col_iter)
347 *col_iter = *col_iter *(row_sum_pos/row_sum_pos_scale);
349 *col_iter = *col_iter *(row_sum_neg/row_sum_neg_scale);
360 template <
typename InternalType1,
typename InternalType2>
361 void amg_interpol_ag(
unsigned int level, InternalType1 & A, InternalType1 & P, InternalType2 & Pointvector,
amg_tag & tag)
363 typedef typename InternalType1::value_type SparseMatrixType;
364 typedef typename InternalType2::value_type PointVectorType;
365 typedef typename SparseMatrixType::value_type ScalarType;
366 typedef typename SparseMatrixType::iterator1 InternalRowIterator;
367 typedef typename SparseMatrixType::iterator2 InternalColIterator;
371 unsigned int c_points = Pointvector[level].get_cpoints();
373 P[level] = SparseMatrixType(A[level].
size1(), c_points);
377 Pointvector[level].build_index();
381 #pragma omp parallel for private (x,pointx) shared (P)
383 for (x=0; x<Pointvector[level].size(); ++x)
385 pointx = Pointvector[level][x];
392 std::cout <<
"Aggregation based Prolongation:" << std::endl;
404 template <
typename InternalType1,
typename InternalType2>
405 void amg_interpol_sa(
unsigned int level, InternalType1 & A, InternalType1 & P, InternalType2 & Pointvector,
amg_tag & tag)
407 typedef typename InternalType1::value_type SparseMatrixType;
408 typedef typename InternalType2::value_type PointVectorType;
409 typedef typename SparseMatrixType::value_type ScalarType;
410 typedef typename SparseMatrixType::iterator1 InternalRowIterator;
411 typedef typename SparseMatrixType::iterator2 InternalColIterator;
415 unsigned int c_points = Pointvector[level].get_cpoints();
417 InternalType1 P_tentative = InternalType1(P.size());
418 SparseMatrixType Jacobi = SparseMatrixType(A[level].
size1(), A[level].
size2());
420 P[level] = SparseMatrixType(A[level].
size1(), c_points);
425 #pragma omp parallel for private (x,y,diag) shared (A,Pointvector)
427 for (x=0; x<A[level].size1(); ++x)
430 InternalRowIterator row_iter = A[level].begin1();
432 for (InternalColIterator col_iter = row_iter.begin(); col_iter != row_iter.end(); ++col_iter)
434 y = col_iter.index2();
441 else if (!Pointvector[level][x]->is_influencing(Pointvector[level][y]))
444 Jacobi (x,y) = *col_iter;
446 InternalRowIterator row_iter2 = Jacobi.begin1();
449 for (InternalColIterator col_iter2 = row_iter2.begin(); col_iter2 != row_iter2.end(); ++col_iter2)
458 std::cout <<
"Jacobi Matrix:" << std::endl;
466 std::cout <<
"Tentative Prolongation:" << std::endl;
474 std::cout <<
"Prolongation Matrix:" << std::endl;