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 #include "anisotropic.hh"
5 #include <dune/common/timer.hh>
6 #include <dune/common/parallel/indexset.hh>
7 #include <dune/common/parallel/communication.hh>
8 #include <dune/istl/paamg/kamg.hh>
9 #include <dune/istl/paamg/pinfo.hh>
10 #include <cstdlib>
11 #include <ctime>
12
13 namespace Dune
14 {
15 using Mat = BCRSMatrix<FieldMatrix<double,1,1>>;
16 using Vec = BlockVector<FieldVector<double,1>>;
17 using LinOp = MatrixAdapter<Mat,Vec,Vec>;
18 using Comm = Amg::SequentialInformation;
19
20 // explicit template instantion of FastAMG preconditioner
21 template class Amg::KAMG<LinOp, Vec, Richardson<Vec,Vec>, Comm, GeneralizedPCGSolver<Vec>>;
22 template class Amg::KAMG<LinOp, Vec, SeqJac<Mat,Vec,Vec>, Comm, GeneralizedPCGSolver<Vec>>;
23 template class Amg::KAMG<LinOp, Vec, SeqSOR<Mat,Vec,Vec>, Comm, GeneralizedPCGSolver<Vec>>;
24 template class Amg::KAMG<LinOp, Vec, SeqSSOR<Mat,Vec,Vec>, Comm, GeneralizedPCGSolver<Vec>>;
25
26 template class Amg::KAMG<LinOp, Vec, SeqJac<Mat,Vec,Vec>, Comm, RestartedFCGSolver<Vec>>;
27 template class Amg::KAMG<LinOp, Vec, SeqJac<Mat,Vec,Vec>, Comm, CompleteFCGSolver<Vec>>;
28 template class Amg::KAMG<LinOp, Vec, SeqJac<Mat,Vec,Vec>, Comm, BiCGSTABSolver<Vec>>;
29 template class Amg::KAMG<LinOp, Vec, SeqJac<Mat,Vec,Vec>, Comm, MINRESSolver<Vec>>;
30
31 /**
32 * Note: RestartedGMResSolver can not be used, since additional constructor argument required.
33 * Needs something like ConstructionTraits for solvers.
34 **/
35 // template class Amg::KAMG<LinOp, Vec, SeqJac<Matrix,Vec,Vec>, Comm, RestartedGMResSolver<Vec>>;
36
37 } // end namespace Dune
38
39
40 template<class M, class V>
randomize(const M & mat,V & b)41 void randomize(const M& mat, V& b)
42 {
43 V x=b;
44
45 srand((unsigned)std::clock());
46
47 typedef typename V::iterator iterator;
48 for(iterator i=x.begin(); i != x.end(); ++i)
49 *i=(rand() / (RAND_MAX + 1.0));
50
51 mat.mv(static_cast<const V&>(x), b);
52 }
53
54
55 template <int BS>
testAMG(int N,int coarsenTarget,int ml)56 void testAMG(int N, int coarsenTarget, int ml)
57 {
58
59 std::cout<<"N="<<N<<" coarsenTarget="<<coarsenTarget<<" maxlevel="<<ml<<std::endl;
60
61
62 typedef Dune::ParallelIndexSet<int,LocalIndex,512> ParallelIndexSet;
63
64 ParallelIndexSet indices;
65 typedef Dune::FieldMatrix<double,BS,BS> MatrixBlock;
66 typedef Dune::BCRSMatrix<MatrixBlock> BCRSMat;
67 typedef Dune::FieldVector<double,BS> VectorBlock;
68 typedef Dune::BlockVector<VectorBlock> Vector;
69 typedef Dune::MatrixAdapter<BCRSMat,Vector,Vector> Operator;
70 typedef Dune::CollectiveCommunication<void*> Comm;
71 int n;
72
73 Comm c;
74 BCRSMat mat = setupAnisotropic2d<MatrixBlock>(N, indices, c, &n, .001);
75
76 Vector b(mat.N()), x(mat.M());
77
78 b=0;
79 x=100;
80
81 setBoundary(x, b, N);
82
83 x=0;
84 randomize(mat, b);
85
86 if(N<6) {
87 Dune::printmatrix(std::cout, mat, "A", "row");
88 Dune::printvector(std::cout, x, "x", "row");
89 }
90
91 Dune::Timer watch;
92
93 watch.reset();
94 Operator fop(mat);
95
96 typedef Dune::Amg::CoarsenCriterion<Dune::Amg::UnSymmetricCriterion<BCRSMat,Dune::Amg::FirstDiagonal> >
97 Criterion;
98 typedef Dune::SeqSSOR<BCRSMat,Vector,Vector> Smoother;
99 //typedef Dune::SeqSOR<BCRSMat,Vector,Vector> Smoother;
100 //typedef Dune::SeqJac<BCRSMat,Vector,Vector> Smoother;
101 //typedef Dune::SeqOverlappingSchwarz<BCRSMat,Vector,Dune::MultiplicativeSchwarzMode> Smoother;
102 //typedef Dune::SeqOverlappingSchwarz<BCRSMat,Vector,Dune::SymmetricMultiplicativeSchwarzMode> Smoother;
103 //typedef Dune::SeqOverlappingSchwarz<BCRSMat,Vector> Smoother;
104 typedef typename Dune::Amg::SmootherTraits<Smoother>::Arguments SmootherArgs;
105
106 SmootherArgs smootherArgs;
107
108 smootherArgs.iterations = 1;
109
110 //smootherArgs.overlap=SmootherArgs::vertex;
111 //smootherArgs.overlap=SmootherArgs::none;
112 //smootherArgs.overlap=SmootherArgs::aggregate;
113
114 smootherArgs.relaxationFactor = 1;
115
116 Criterion criterion(15,coarsenTarget);
117 criterion.setDefaultValuesIsotropic(2);
118 criterion.setAlpha(.67);
119 criterion.setBeta(1.0e-4);
120 criterion.setMaxLevel(ml);
121 criterion.setSkipIsolated(false);
122 // specify pre/post smoother steps
123 criterion.setNoPreSmoothSteps(1);
124 criterion.setNoPostSmoothSteps(1);
125
126 Dune::SeqScalarProduct<Vector> sp;
127 typedef Dune::Amg::KAMG<Operator,Vector,Smoother,Dune::Amg::SequentialInformation> AMG;
128
129 AMG amg(fop, criterion, smootherArgs);
130
131
132 double buildtime = watch.elapsed();
133
134 std::cout<<"Building hierarchy took "<<buildtime<<" seconds"<<std::endl;
135
136 //Dune::BiCGSTABSolver<Vector> amgCG(fop,amg,1e-6,80,2);
137 Dune::CGSolver<Vector> amgCG(fop,amg,1e-6,80,2);
138 //Dune::LoopSolver<Vector> amgCG(fop, amg, 1e-4, 10000, 2);
139 watch.reset();
140 Dune::InverseOperatorResult r;
141 amgCG.apply(x,b,r);
142
143 double solvetime = watch.elapsed();
144
145 std::cout<<"AMG solving took "<<solvetime<<" seconds"<<std::endl;
146
147 std::cout<<"AMG building took "<<(buildtime/r.elapsed*r.iterations)<<" iterations"<<std::endl;
148 std::cout<<"AMG building together with solving took "<<buildtime+solvetime<<std::endl;
149
150 /*
151 watch.reset();
152 cg.apply(x,b,r);
153
154 std::cout<<"CG solving took "<<watch.elapsed()<<" seconds"<<std::endl;
155 */
156 }
157
158
main(int argc,char ** argv)159 int main(int argc, char** argv)
160 try
161 {
162 int N=100;
163 int coarsenTarget=1200;
164 int ml=10;
165
166 if(argc>1)
167 N = atoi(argv[1]);
168
169 if(argc>2)
170 coarsenTarget = atoi(argv[2]);
171
172 if(argc>3)
173 ml = atoi(argv[3]);
174
175 testAMG<1>(N, coarsenTarget, ml);
176 testAMG<2>(N, coarsenTarget, ml);
177
178 return 0;
179 }
180 catch (std::exception &e)
181 {
182 std::cout << "ERROR: " << e.what() << std::endl;
183 return 1;
184 }
185 catch (...)
186 {
187 std::cerr << "Dune reported an unknown error." << std::endl;
188 exit(1);
189 }
190