1 #ifndef VIENNACL_LINALG_DETAIL_AMG_AMG_BASE_HPP_
2 #define VIENNACL_LINALG_DETAIL_AMG_AMG_BASE_HPP_
26 #include <boost/numeric/ublas/vector.hpp>
39 #define VIENNACL_AMG_COARSE_RS 1
40 #define VIENNACL_AMG_COARSE_ONEPASS 2
41 #define VIENNACL_AMG_COARSE_RS0 3
42 #define VIENNACL_AMG_COARSE_RS3 4
43 #define VIENNACL_AMG_COARSE_AG 5
44 #define VIENNACL_AMG_INTERPOL_DIRECT 1
45 #define VIENNACL_AMG_INTERPOL_CLASSIC 2
46 #define VIENNACL_AMG_INTERPOL_AG 3
47 #define VIENNACL_AMG_INTERPOL_SA 4
75 unsigned int interpol = 1,
76 double threshold = 0.25,
77 double interpolweight = 0.2,
78 double jacobiweight = 1,
79 unsigned int presmooth = 1,
80 unsigned int postsmooth = 1,
81 unsigned int coarselevels = 0)
82 : _coarse(coarse), _interpol(interpol),
83 _threshold(threshold), _interpolweight(interpolweight), _jacobiweight(jacobiweight),
84 _presmooth(presmooth), _postsmooth(postsmooth), _coarselevels(coarselevels) {};
87 void set_coarse(
unsigned int coarse) {
if (coarse > 0) _coarse = coarse; }
90 void set_interpol(
unsigned int interpol) {
if (interpol > 0) _interpol = interpol; }
93 void set_threshold(
double threshold) {
if (threshold > 0 && threshold <= 1) _threshold = threshold; }
96 void set_as(
double jacobiweight) {
if (jacobiweight > 0 && jacobiweight <= 2) _jacobiweight = jacobiweight; }
99 void set_interpolweight(
double interpolweight) {
if (interpolweight > 0 && interpolweight <= 2) _interpolweight = interpolweight; }
102 void set_presmooth(
int presmooth) {
if (presmooth >= 0) _presmooth = presmooth; }
105 void set_postsmooth(
int postsmooth) {
if (postsmooth >= 0) _postsmooth = postsmooth; }
108 void set_coarselevels(
int coarselevels) {
if (coarselevels >= 0) _coarselevels = coarselevels; }
112 unsigned int _coarse, _interpol;
113 double _threshold, _interpolweight, _jacobiweight;
114 unsigned int _presmooth, _postsmooth, _coarselevels;
121 template <
typename InternalType,
typename IteratorType,
typename ScalarType>
144 ScalarType s = 0): _m(m), _iter(iter), _i(i), _j(j), _s(s) {}
153 if (_s == 0)
return _s;
155 _m->addscalar (_iter,_i,_j,_s);
172 _m->removescalar(_iter,_i);
176 _m->addscalar (_iter,_i,_j,_s);
183 _m->removescalar(_iter,_i);
184 _m->addscalar (_iter,_i,_j,_s);
191 _m->removescalar(_iter,_i);
192 _m->addscalar (_iter,_i,_j,_s);
195 operator ScalarType (
void) {
return _s; }
200 template <
typename InternalType>
205 typedef typename InternalType::mapped_type ScalarType;
207 InternalType & internal_vec;
208 typename InternalType::iterator iter;
219 iter = internal_vec.begin();
221 iter = internal_vec.end();
226 if (iter == other.iter)
233 if (iter != other.iter)
245 unsigned int index()
const {
return (*iter).first; }
246 unsigned int index() {
return (*iter).first; }
251 template <
typename ScalarType>
259 typedef std::map<unsigned int,ScalarType> InternalType;
265 InternalType internal_vector;
277 internal_vector = InternalType();
281 unsigned int size()
const {
return _size;}
286 void clear() { internal_vector.clear(); }
288 void remove(
unsigned int i) { internal_vector.erase(i); }
291 void add (
unsigned int i, ScalarType s)
293 typename InternalType::iterator iter = internal_vector.find(i);
295 if (iter == internal_vector.end())
304 internal_vector.erase(iter);
309 template <
typename IteratorType>
310 void addscalar(IteratorType & iter,
unsigned int i,
unsigned int j, ScalarType s)
317 if (iter != internal_vector.end())
320 internal_vector[i] = s;
324 template <
typename IteratorType>
325 void removescalar(IteratorType & iter,
unsigned int i) { internal_vector.erase(iter); }
330 typename InternalType::iterator it = internal_vector.find(i);
332 if (it != internal_vector.end())
343 if (it != internal_vector.end())
356 bool isnonzero(
unsigned int i)
const {
return internal_vector.find(i) != internal_vector.end(); }
359 operator boost::numeric::ublas::vector<ScalarType> (void)
361 boost::numeric::ublas::vector<ScalarType> vec (_size);
363 vec [iter.index()] = *iter;
373 template <
typename ScalarType>
379 typedef std::map<unsigned int,ScalarType> RowType;
380 typedef std::vector<RowType> InternalType;
391 InternalType internal_mat;
394 InternalType internal_mat_trans;
399 bool transposed_mode;
412 transposed_mode =
false;
425 a_trans.
resize (j,i,
false);
430 transposed_mode =
false;
442 a.
resize(mat.size(), mat.size());
443 a_trans.
resize(mat.size(), mat.size());
446 s1 = s2 = mat.size();
448 transposed_mode =
false;
456 template <
typename MatrixType>
461 a.
resize(mat.size1(), mat.size2());
462 a_trans.
resize (mat.size2(), mat.size1());
468 for (
typename MatrixType::const_iterator1 row_iter = mat.begin1(); row_iter != mat.end1(); ++row_iter)
470 for (
typename MatrixType::const_iterator2 col_iter = row_iter.begin(); col_iter != row_iter.end(); ++col_iter)
474 unsigned int x = col_iter.index1();
475 unsigned int y = col_iter.index2();
477 a_trans (y,x) = *col_iter;
481 transposed_mode =
false;
497 bool save_mode = transposed_mode;
498 transposed_mode =
false;
501 for (
iterator2 col_iter = row_iter.
begin(); col_iter != row_iter.end(); ++col_iter)
502 internal_mat_trans[col_iter.index2()][col_iter.index1()] = *col_iter;
504 transposed_mode = save_mode;
513 transposed_mode = mode;
523 if (!transposed_mode)
525 if (internal_mat[i].find(j) != internal_mat[i].end())
532 if (internal_mat[j].find(i) != internal_mat[j].end())
540 void add (
unsigned int i,
unsigned int j, ScalarType s)
546 typename RowType::iterator col_iter = internal_mat[i].find(j);
548 if (col_iter == internal_mat[i].end())
552 s += (*col_iter).second;
555 (*col_iter).second = s;
557 internal_mat[i].erase(col_iter);
563 template <
typename IteratorType>
564 void addscalar(IteratorType & iter,
unsigned int i,
unsigned int j, ScalarType s)
570 if (iter != internal_mat[i].end())
573 internal_mat[i][j] = s;
579 template <
typename IteratorType>
582 internal_mat[i].erase(iter);
589 typename RowType::iterator iter;
591 if (!transposed_mode)
593 iter = internal_mat[i].find(j);
594 if (iter != internal_mat[i].end())
601 iter = internal_mat[j].find(i);
602 if (iter != internal_mat[j].end())
612 typename RowType::const_iterator iter;
614 if (!transposed_mode)
616 iter = internal_mat[i].find(j);
617 if (iter != internal_mat[i].end())
618 return (*iter).second;
624 iter = internal_mat[j].find(i);
625 if (iter != internal_mat[j].end())
626 return (*iter).second;
632 void resize(
unsigned int i,
unsigned int j,
bool preserve =
true)
637 a_trans.
resize(j,i,preserve);
653 if (!transposed_mode)
661 if (!transposed_mode)
670 if (!transposed_mode)
677 if (!transposed_mode)
685 if (!
trans && !transposed_mode)
700 if (!
trans && !transposed_mode)
709 return a_trans.
end1();
715 if (!
trans && !transposed_mode)
730 if (!
trans && !transposed_mode)
739 return a_trans.
end2();
746 assert((!transposed_mode || (transposed_mode && transposed)) &&
"Error: Cannot build const_iterator when transposed has not been built yet!");
753 assert((!transposed_mode || (transposed_mode && transposed)) &&
"Error: Cannot build const_iterator when transposed has not been built yet!");
755 return a_const.
end1();
760 assert((!transposed_mode || (transposed_mode && transposed)) &&
"Error: Cannot build const_iterator when transposed has not been built yet!");
767 assert((!transposed_mode || (transposed_mode && transposed)) &&
"Error: Cannot build const_iterator when transposed has not been built yet!");
769 return a_const.
end2();
775 if (!transposed_mode)
776 return &internal_mat;
780 return &internal_mat_trans;
782 operator boost::numeric::ublas::compressed_matrix<ScalarType> (void)
784 boost::numeric::ublas::compressed_matrix<ScalarType> mat;
789 for (
iterator2 col_iter = row_iter.
begin(); col_iter != row_iter.end(); ++col_iter)
790 mat (col_iter.index1(), col_iter.index2()) = *col_iter;
794 operator boost::numeric::ublas::matrix<ScalarType> (void)
796 boost::numeric::ublas::matrix<ScalarType> mat;
801 for (
iterator2 col_iter = row_iter.
begin(); col_iter != row_iter.end(); ++col_iter)
802 mat (col_iter.index1(), col_iter.index2()) = *col_iter;
819 unsigned int _influence;
824 unsigned int _coarse_index;
826 unsigned int _offset;
828 unsigned int _aggregate;
841 amg_point (
unsigned int index,
unsigned int size): _index(index), _influence(0), _undecided(true), _cpoint(false), _coarse_index(0), _offset(0), _aggregate(0)
843 influencing_points =
ListType(size);
849 void set_index(
unsigned int index) { _index = index+_offset; }
850 unsigned int get_index()
const {
return _index-_offset; }
855 bool is_cpoint()
const {
return _cpoint && !_undecided; }
856 bool is_fpoint()
const {
return !_cpoint && !_undecided; }
936 typedef std::set<amg_point*,classcomp> ListType;
938 typedef std::vector<amg_point*> VectorType;
939 VectorType pointvector;
942 unsigned int c_points, f_points;
953 pointvector = VectorType(
size);
954 c_points = f_points = 0;
960 for (
unsigned int i=0; i<
size(); ++i)
966 for (
unsigned int i=0; i<
size(); ++i)
967 delete pointvector[i];
990 for (
iterator iter = pointvector.begin(); iter != pointvector.end(); ++iter)
992 (*iter)->clear_influencing();
993 (*iter)->clear_influenced();
1006 pointvector = VectorType(size);
1008 unsigned int size()
const {
return _size; }
1019 pointlist.insert(*iter);
1025 if (pointlist.size() == 0)
1028 if ((*(--pointlist.end()))->get_influence() == 0)
1032 return *(--pointlist.end());
1037 ListType::iterator iter = pointlist.find(point);
1039 if (iter == pointlist.end())
return;
1042 ListType::iterator iter2 = iter;
1046 pointlist.erase(iter);
1050 pointlist.insert(iter2,point);
1055 pointlist.erase(point);
1062 pointlist.erase(point);
1077 unsigned int count = 0;
1079 for (
iterator iter = pointvector.begin(); iter != pointvector.end(); ++iter)
1082 if ((*iter)->is_cpoint())
1084 (*iter)->set_coarse_index(count);
1091 template <
typename MatrixType>
1098 for (
amg_point::iterator col_iter = (*row_iter)->begin_influencing(); col_iter != (*row_iter)->end_influencing(); ++col_iter)
1099 mat((*row_iter)->get_index(),(*col_iter)->get_index()) =
true;
1101 template <
typename VectorType>
1104 vec = VectorType(_size);
1108 vec[(*iter)->get_index()] = (*iter)->get_influence();
1110 template <
typename VectorType>
1113 vec = VectorType(pointlist.size());
1117 for (ListType::const_iterator iter = pointlist.begin(); iter != pointlist.end(); ++iter)
1119 vec[i] = (*iter)->get_index();
1123 template <
typename VectorType>
1126 vec = VectorType(_size);
1131 if ((*iter)->is_cpoint())
1132 vec[(*iter)->get_index()] =
true;
1135 template <
typename VectorType>
1138 vec = VectorType(_size);
1143 if ((*iter)->is_fpoint())
1144 vec[(*iter)->get_index()] =
true;
1147 template <
typename MatrixType>
1150 mat = MatrixType(_size,_size);
1155 if (!(*iter)->is_undecided())
1156 mat((*iter)->get_aggregate(),(*iter)->get_index()) =
true;
1164 template <
typename InternalType1,
typename InternalType2>
1167 typedef typename InternalType1::value_type SparseMatrixType;
1168 typedef typename InternalType2::value_type PointVectorType;
1172 boost::numeric::ublas::vector<InternalType1>
A_slice;
1175 boost::numeric::ublas::vector<boost::numeric::ublas::vector<unsigned int> >
Offset;
1180 void init(
unsigned int levels,
unsigned int threads = 0)
1199 for (
unsigned int i=0; i<
_threads; ++i)
1210 void slice (
unsigned int level, InternalType1
const & A, InternalType2
const & Pointvector)
1214 slice_new (level, A);
1219 slice_build (level, A, Pointvector);
1223 void join (
unsigned int level, InternalType2 & Pointvector)
const
1225 typedef typename InternalType2::value_type PointVectorType;
1228 Pointvector[level].clear_cf();
1229 for (
typename PointVectorType::iterator iter = Pointvector[level].begin(); iter != Pointvector[level].end(); ++iter)
1231 (*iter)->set_offset(0);
1232 Pointvector[level].update_cf(*iter);
1241 void slice_new (
unsigned int level, InternalType1
const & A)
1243 typedef typename SparseMatrixType::const_iterator1 ConstRowIterator;
1244 typedef typename SparseMatrixType::const_iterator2 ConstColIterator;
1248 #pragma omp parallel for
1250 for (
unsigned int i=0; i<=
_threads; ++i)
1253 if (i == 0)
Offset[i][level] = 0;
1264 void slice_build (
unsigned int level, InternalType1
const & A, InternalType2
const & Pointvector)
1266 typedef typename SparseMatrixType::const_iterator1 ConstRowIterator;
1267 typedef typename SparseMatrixType::const_iterator2 ConstColIterator;
1273 #pragma omp parallel for private (x,y,point)
1275 for (
unsigned int i=0; i<
_threads; ++i)
1282 ConstRowIterator row_iter = A[level].begin1();
1283 row_iter +=
Offset[i][level];
1284 x = row_iter.index1();
1286 while (x <
Offset[i+1][level] && row_iter != A[level].end1())
1289 point = Pointvector[level][x];
1290 point->set_offset(
Offset[i][level]);
1293 ConstColIterator col_iter = row_iter.begin();
1294 y = col_iter.index2();
1297 while (y <
Offset[i+1][level] && col_iter != row_iter.end())
1299 if (y >=
Offset[i][level])
1300 A_slice[i][level](x-
Offset[i][level],y-Offset[i][level]) = *col_iter;
1303 y = col_iter.index2();
1307 x = row_iter.index1();
1318 template <
typename SparseMatrixType>
1319 void amg_mat_prod (SparseMatrixType & A, SparseMatrixType & B, SparseMatrixType & RES)
1321 typedef typename SparseMatrixType::value_type ScalarType;
1322 typedef typename SparseMatrixType::iterator1 InternalRowIterator;
1323 typedef typename SparseMatrixType::iterator2 InternalColIterator;
1327 RES = SparseMatrixType(A.size1(), B.size2());
1331 #pragma omp parallel for private (x,y,z,prod) shared (A,B,RES)
1333 for (x=0; x<A.size1(); ++x)
1335 InternalRowIterator row_iter = A.begin1();
1337 for (InternalColIterator col_iter = row_iter.begin(); col_iter != row_iter.end(); ++col_iter)
1339 y = col_iter.index2();
1340 InternalRowIterator row_iter2 = B.begin1();
1343 for(InternalColIterator col_iter2 = row_iter2.begin(); col_iter2 != row_iter2.end(); ++col_iter2)
1345 z = col_iter2.index2();
1346 prod = *col_iter * *col_iter2;
1358 template <
typename SparseMatrixType>
1361 typedef typename SparseMatrixType::value_type ScalarType;
1362 typedef typename SparseMatrixType::iterator1 InternalRowIterator;
1363 typedef typename SparseMatrixType::iterator2 InternalColIterator;
1365 unsigned int x,y1,y2,z;
1367 RES = SparseMatrixType(P.size2(), P.size2());
1371 #pragma omp parallel for private (x,y1,y2,z,row) shared (A,P,RES)
1373 for (x=0; x<P.size2(); ++x)
1376 InternalRowIterator row_iter = P.begin1(
true);
1379 for (InternalColIterator col_iter = row_iter.
begin(); col_iter != row_iter.
end(); ++col_iter)
1381 y1 = col_iter.index2();
1382 InternalRowIterator row_iter2 = A.begin1();
1385 for(InternalColIterator col_iter2 = row_iter2.begin(); col_iter2 != row_iter2.end(); ++col_iter2)
1387 y2 = col_iter2.index2();
1388 row.
add (y2, *col_iter * *col_iter2);
1394 InternalRowIterator row_iter3 = P.begin1();
1397 for (InternalColIterator col_iter3 = row_iter3.begin(); col_iter3 != row_iter3.end(); ++col_iter3)
1399 z = col_iter3.index2();
1400 RES.add (x, z, *col_iter3 * *iter);
1406 std::cout <<
"Galerkin Operator: " << std::endl;
1416 template <
typename SparseMatrixType>
1419 typedef typename SparseMatrixType::value_type ScalarType;
1421 boost::numeric::ublas::compressed_matrix<ScalarType> A_temp (A.size1(), A.size2());
1423 boost::numeric::ublas::compressed_matrix<ScalarType> P_temp (P.size1(), P.size2());
1426 boost::numeric::ublas::compressed_matrix<ScalarType> R_temp (P.size1(), P.size2());
1430 boost::numeric::ublas::compressed_matrix<ScalarType> RA (R_temp.size1(),A_temp.size2());
1432 boost::numeric::ublas::compressed_matrix<ScalarType> RAP (RA.size1(),P_temp.size2());
1435 for (
unsigned int x=0; x<RAP.size1(); ++x)
1437 for (
unsigned int y=0; y<RAP.size2(); ++y)
1439 if (abs((ScalarType)RAP(x,y) - (ScalarType)A_i1(x,y)) > 0.0001)
1440 std::cout << x <<
" " << y <<
" " << RAP(x,y) <<
" " << A_i1(x,y) << std::endl;
1450 template <
typename SparseMatrixType,
typename Po
intVectorType>
1453 for (
unsigned int i=0; i<P.size1(); ++i)
1455 if (Pointvector.is_cpoint(i))
1458 for (
unsigned int j=0; j<P.size2(); ++j)
1460 if (P.isnonzero(i,j))
1463 std::cout <<
"Error 1 in row " << i << std::endl;
1464 if (P(i,j) == 1 &&
set)
1465 std::cout <<
"Error 2 in row " << i << std::endl;
1466 if (P(i,j) == 1 && !
set)
1472 if (Pointvector.is_fpoint(i))
1473 for (
unsigned int j=0; j<P.size2(); ++j)
1475 if (P.isnonzero(i,j) && j> Pointvector.get_cpoints()-1)
1476 std::cout <<
"Error 3 in row " << i << std::endl;
1477 if (P.isnonzero(i,j))
1480 for (
unsigned int k=0; k<P.size1(); ++k)
1482 if (P.isnonzero(k,j))
1484 if (Pointvector.is_cpoint(k) && P(k,j) == 1 && A.isnonzero(i,k))
1489 std::cout <<
"Error 4 in row " << i << std::endl;