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"xreal.h"
6
7 namespace std
8 {
9
10 HPA::xreal abs(const HPA::xreal& t)
11 {
12 return t>=0?t:-t;
13 }
14
15 };
16 */
17
18 #include "anisotropic.hh"
19 #include <dune/common/timer.hh>
20 #include <dune/common/parallel/indexset.hh>
21 #include <dune/common/parallel/communication.hh>
22 #include <dune/istl/paamg/amg.hh>
23 #include <dune/istl/paamg/fastamg.hh>
24 #include <dune/istl/paamg/pinfo.hh>
25 #include <dune/istl/solvers.hh>
26 #include <cstdlib>
27 #include <ctime>
28 #include <pthread.h>
29 #define NUM_THREADS 3
30
31 typedef double XREAL;
32
33 /*
34 typedef HPA::xreal XREAL;
35
36 namespace Dune
37 {
38 template<>
39 struct DoubleConverter<HPA::xreal>
40 {
41 static double toDouble(const HPA::xreal& t)
42 {
43 return t._2double();
44 }
45 };
46 }
47 */
48
49 template<class M, class V>
randomize(const M & mat,V & b)50 void randomize(const M& mat, V& b)
51 {
52 V x=b;
53
54 srand((unsigned)std::clock());
55
56 typedef typename V::iterator iterator;
57 for(iterator i=x.begin(); i != x.end(); ++i)
58 *i=(rand() / (RAND_MAX + 1.0));
59
60 mat.mv(static_cast<const V&>(x), b);
61 }
62
63 typedef Dune::ParallelIndexSet<int,LocalIndex,512> ParallelIndexSet;
64 typedef Dune::FieldMatrix<XREAL,1,1> MatrixBlock;
65 typedef Dune::BCRSMatrix<MatrixBlock> BCRSMat;
66 typedef Dune::FieldVector<XREAL,1> VectorBlock;
67 typedef Dune::BlockVector<VectorBlock> Vector;
68 typedef Dune::MatrixAdapter<BCRSMat,Vector,Vector> Operator;
69 typedef Dune::CollectiveCommunication<void*> Comm;
70 typedef Dune::SeqSSOR<BCRSMat,Vector,Vector> Smoother;
71 typedef Dune::Amg::SmootherTraits<Smoother>::Arguments SmootherArgs;
72
73 struct thread_arg
74 {
75 MYAMG *amg;
76 Vector *b;
77 Vector *x;
78 Operator *fop;
79 };
80
81
solve(void * arg)82 void *solve(void* arg)
83 {
84 thread_arg *amgarg=(thread_arg*) arg;
85
86 Dune::GeneralizedPCGSolver<Vector> amgCG(*amgarg->fop,*amgarg->amg,1e-6,80,2);
87 //Dune::LoopSolver<Vector> amgCG(fop, amg, 1e-4, 10000, 2);
88 Dune::Timer watch;
89 Dune::InverseOperatorResult r;
90 amgCG.apply(*amgarg->x,*amgarg->b,r);
91
92 double solvetime = watch.elapsed();
93
94 std::cout<<"AMG solving took "<<solvetime<<" seconds"<<std::endl;
95
96 pthread_exit(NULL);
97 }
98
solve2(void * arg)99 void *solve2(void* arg)
100 {
101 thread_arg *amgarg=(thread_arg*) arg;
102 *amgarg->x=0;
103 (*amgarg->amg).pre(*amgarg->x,*amgarg->b);
104 (*amgarg->amg).apply(*amgarg->x,*amgarg->b);
105 (*amgarg->amg).post(*amgarg->x);
106 return 0;
107 }
108
109 template <int BS, typename AMG>
testAMG(int N,int coarsenTarget,int ml)110 void testAMG(int N, int coarsenTarget, int ml)
111 {
112
113 std::cout<<"N="<<N<<" coarsenTarget="<<coarsenTarget<<" maxlevel="<<ml<<std::endl;
114
115
116
117 ParallelIndexSet indices;
118 int n;
119
120 Comm c;
121 BCRSMat mat = setupAnisotropic2d<MatrixBlock>(N, indices, c, &n, 1);
122
123 Vector b(mat.N()), x(mat.M());
124
125 b=0;
126 x=100;
127
128 setBoundary(x, b, N);
129
130 x=0;
131 randomize(mat, b);
132
133 std::vector<Vector> xs(NUM_THREADS, x);
134 std::vector<Vector> bs(NUM_THREADS, b);
135
136 if(N<6) {
137 Dune::printmatrix(std::cout, mat, "A", "row");
138 Dune::printvector(std::cout, x, "x", "row");
139 }
140
141 Dune::Timer watch;
142
143 watch.reset();
144 Operator fop(mat);
145
146 typedef Dune::Amg::CoarsenCriterion<Dune::Amg::UnSymmetricCriterion<BCRSMat,Dune::Amg::FirstDiagonal> >
147 Criterion;
148 //typedef Dune::SeqSOR<BCRSMat,Vector,Vector> Smoother;
149 //typedef Dune::SeqJac<BCRSMat,Vector,Vector> Smoother;
150 //typedef Dune::SeqOverlappingSchwarz<BCRSMat,Vector,Dune::MultiplicativeSchwarzMode> Smoother;
151 //typedef Dune::SeqOverlappingSchwarz<BCRSMat,Vector,Dune::SymmetricMultiplicativeSchwarzMode> Smoother;
152 //typedef Dune::SeqOverlappingSchwarz<BCRSMat,Vector> Smoother;
153
154 SmootherArgs smootherArgs;
155
156 smootherArgs.iterations = 1;
157
158 //smootherArgs.overlap=SmootherArgs::vertex;
159 //smootherArgs.overlap=SmootherArgs::none;
160 //smootherArgs.overlap=SmootherArgs::aggregate;
161
162 smootherArgs.relaxationFactor = 1;
163
164 Criterion criterion(15,coarsenTarget);
165 criterion.setDefaultValuesIsotropic(2);
166 criterion.setAlpha(.67);
167 criterion.setBeta(1.0e-4);
168 criterion.setMaxLevel(ml);
169 criterion.setSkipIsolated(false);
170
171 Dune::SeqScalarProduct<Vector> sp;
172
173 Smoother smoother(mat,1,1);
174 AMG amg(fop, criterion);
175
176
177 double buildtime = watch.elapsed();
178
179 std::cout<<"Building hierarchy took "<<buildtime<<" seconds"<<std::endl;
180 std::vector<AMG> amgs(NUM_THREADS, amg);
181 std::vector<thread_arg> args(NUM_THREADS);
182 std::vector<pthread_t> threads(NUM_THREADS);
183 for(int i=0; i < NUM_THREADS; ++i)
184 {
185 args[i].amg=&amgs[i];
186 args[i].b=&bs[i];
187 args[i].x=&xs[i];
188 args[i].fop=&fop;
189 pthread_create(&threads[i], NULL, solve, (void*) &args[i]);
190 }
191 void* retval;
192
193 for(int i=0; i < NUM_THREADS; ++i)
194 pthread_join(threads[i], &retval);
195
196 amgs.clear();
197 args.clear();
198 amg.pre(x, b);
199 amgs.resize(NUM_THREADS, amg);
200 for(int i=0; i < NUM_THREADS; ++i)
201 {
202 args[i].amg=&amgs[i];
203 args[i].b=&bs[i];
204 args[i].x=&xs[i];
205 args[i].fop=&fop;
206 pthread_create(&threads[i], NULL, solve2, (void*) &args[i]);
207 }
208 for(int i=0; i < NUM_THREADS; ++i)
209 pthread_join(threads[i], &retval);
210
211 }
212
213
main(int argc,char ** argv)214 int main(int argc, char** argv)
215 try
216 {
217 int N=100;
218 int coarsenTarget=1200;
219 int ml=10;
220
221 if(argc>1)
222 N = atoi(argv[1]);
223
224 if(argc>2)
225 coarsenTarget = atoi(argv[2]);
226
227 if(argc>3)
228 ml = atoi(argv[3]);
229
230 testAMG<1,MYAMG>(N, coarsenTarget, ml);
231
232 //testAMG<2>(N, coarsenTarget, ml);
233
234 return 0;
235 }
236 catch (std::exception &e)
237 {
238 std::cout << "ERROR: " << e.what() << std::endl;
239 return 1;
240 }
241 catch (...)
242 {
243 std::cerr << "Dune reported an unknown error." << std::endl;
244 exit(1);
245 }
246