1 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=4 sw=2 sts=2:
3 #include "config.h"
4
5 #include "anisotropic.hh"
6 #include <dune/common/timer.hh>
7 #include <dune/common/parallel/indexset.hh>
8 #include <dune/common/parallel/communication.hh>
9 #include <dune/istl/paamg/fastamg.hh>
10 #include <dune/istl/paamg/pinfo.hh>
11 #include <dune/istl/solvers.hh>
12 #include <cstdlib>
13 #include <ctime>
14
15 namespace Dune
16 {
17 using Mat = BCRSMatrix<FieldMatrix<double,1,1>>;
18 using Vec = BlockVector<FieldVector<double,1>>;
19
20 // explicit template instantion of FastAMG preconditioner
21 template class Amg::FastAMG<MatrixAdapter<Mat,Vec,Vec>, Vec, Amg::SequentialInformation>;
22
23 } // end namespace Dune
24
25
26 template<class M, class V>
randomize(const M & mat,V & b)27 void randomize(const M& mat, V& b)
28 {
29 V x=b;
30
31 srand((unsigned)std::clock());
32
33 typedef typename V::iterator iterator;
34 for(iterator i=x.begin(); i != x.end(); ++i)
35 *i=(rand() / (RAND_MAX + 1.0));
36
37 mat.mv(static_cast<const V&>(x), b);
38 }
39
40 template <int BS>
testAMG(int N,int coarsenTarget,int ml)41 void testAMG(int N, int coarsenTarget, int ml)
42 {
43 std::cout<<"N="<<N<<" coarsenTarget="<<coarsenTarget<<" maxlevel="<<ml<<std::endl;
44
45 typedef Dune::ParallelIndexSet<int,LocalIndex,512> ParallelIndexSet;
46
47 ParallelIndexSet indices;
48 typedef Dune::FieldMatrix<double,BS,BS> MatrixBlock;
49 typedef Dune::BCRSMatrix<MatrixBlock> BCRSMat;
50 typedef Dune::FieldVector<double,BS> VectorBlock;
51 typedef Dune::BlockVector<VectorBlock> Vector;
52 typedef Dune::MatrixAdapter<BCRSMat,Vector,Vector> Operator;
53 typedef Dune::CollectiveCommunication<void*> Comm;
54 int n;
55
56 Comm c;
57 BCRSMat mat = setupAnisotropic2d<MatrixBlock>(N, indices, c, &n, 1);
58
59 Vector b(mat.N()), x(mat.M());
60
61 b=0;
62 x=100;
63
64 setBoundary(x, b, N);
65
66 x=0;
67 randomize(mat, b);
68
69 if(N<6) {
70 Dune::printmatrix(std::cout, mat, "A", "row");
71 Dune::printvector(std::cout, x, "x", "row");
72 }
73
74 Dune::Timer watch;
75
76 watch.reset();
77 Operator fop(mat);
78
79 typedef Dune::Amg::AggregationCriterion<Dune::Amg::SymmetricMatrixDependency<BCRSMat,Dune::Amg::FirstDiagonal> > CriterionBase;
80 typedef Dune::Amg::CoarsenCriterion<CriterionBase> Criterion;
81
82 Criterion criterion(15,coarsenTarget);
83 criterion.setDefaultValuesIsotropic(2);
84 criterion.setAlpha(.67);
85 criterion.setBeta(1.0e-4);
86 criterion.setMaxLevel(ml);
87 criterion.setSkipIsolated(false);
88
89 Dune::SeqScalarProduct<Vector> sp;
90 typedef Dune::Amg::FastAMG<Operator,Vector> AMG;
91 Dune::Amg::Parameters parms;
92
93 AMG amg(fop, criterion, parms);
94
95 // check if recalculation of matrix hierarchy works
96 amg.recalculateHierarchy();
97
98 double buildtime = watch.elapsed();
99
100 std::cout<<"Building hierarchy took "<<buildtime<<" seconds"<<std::endl;
101
102 Dune::GeneralizedPCGSolver<Vector> amgCG(fop,amg,1e-6,80,2);
103 //Dune::LoopSolver<Vector> amgCG(fop, amg, 1e-4, 10000, 2);
104 watch.reset();
105 Dune::InverseOperatorResult r;
106 amgCG.apply(x,b,r);
107
108 double solvetime = watch.elapsed();
109
110 std::cout<<"AMG solving took "<<solvetime<<" seconds"<<std::endl;
111
112 std::cout<<"AMG building took "<<(buildtime/r.elapsed*r.iterations)<<" iterations"<<std::endl;
113 std::cout<<"AMG building together with solving took "<<buildtime+solvetime<<std::endl;
114
115 /*
116 watch.reset();
117 cg.apply(x,b,r);
118
119 std::cout<<"CG solving took "<<watch.elapsed()<<" seconds"<<std::endl;
120 */
121 }
122
123
main(int argc,char ** argv)124 int main(int argc, char** argv)
125 try
126 {
127 int N=100;
128 int coarsenTarget=1200;
129 int ml=10;
130
131 if(argc>1)
132 N = atoi(argv[1]);
133
134 if(argc>2)
135 coarsenTarget = atoi(argv[2]);
136
137 if(argc>3)
138 ml = atoi(argv[3]);
139
140 testAMG<1>(N, coarsenTarget, ml);
141 testAMG<2>(N, coarsenTarget, ml);
142
143 return 0;
144 }
145 catch (std::exception &e)
146 {
147 std::cout << "ERROR: " << e.what() << std::endl;
148 return 1;
149 }
150 catch (...)
151 {
152 std::cerr << "Dune reported an unknown error." << std::endl;
153 exit(1);
154 }
155