1 //ff-c++-LIBRARY-dep: cxx11 [mkl|blas] mpi pthread htool
2 //ff-c++-cpp-dep:
3 // for def M_PI under windows in <cmath>
4 #define _USE_MATH_DEFINES
5 #include <ff++.hpp>
6 #include <AFunction_ext.hpp>
7
8 #include <htool/lrmat/partialACA.hpp>
9 #include <htool/types/matrix.hpp>
10 #include <htool/types/hmatrix.hpp>
11
12 #include "PlotStream.hpp"
13
14 extern FILE *ThePlotStream;
15
16 //#include "solve.hpp"
17
18 using namespace std;
19 using namespace htool;
20
21 template<class K>
22 class MyMatrix: public IMatrix<K>{
23 const Fem2D::Mesh & ThU; // line
24 const Fem2D::Mesh & ThV; // colunm
25
26 public:
MyMatrix(const FESpace * Uh,const FESpace * Vh)27 MyMatrix(const FESpace * Uh , const FESpace * Vh ):IMatrix<K>(Uh->Th.nv,Vh->Th.nv),ThU(Uh->Th), ThV(Vh->Th) {}
28
get_coef(const int & i,const int & j) const29 K get_coef(const int& i, const int& j)const {return 1./(0.01+(ThU.vertices[i]-ThV.vertices[j]).norme2());}
30
31 };
32
33 template<template<class> class LR, class K>
34 class init_Op : public E_F0mps {
35 public:
36 Expression A;
37 static const int n_name_param = 1;
38 static basicAC_F0::name_and_type name_param[];
39 Expression nargs[n_name_param];
init_Op(const basicAC_F0 & args,Expression param1)40 init_Op(const basicAC_F0& args, Expression param1) : A(param1) {
41 args.SetNameParam(n_name_param, name_param, nargs);
42 }
43
44 AnyType operator()(Stack stack) const;
45 };
46
47 template<template<class> class LR, class K>
48 basicAC_F0::name_and_type init_Op<LR, K>::name_param[] = {
49 };
50
51 template<template<class> class LR, class K>
52 class init : public OneOperator {
53 public:
init()54 init() : OneOperator(atype<HMatrix<LR,K>**>(), atype<HMatrix<LR,K>**>(), atype<Matrice_Creuse<K>*>()) { }
55
code(const basicAC_F0 & args) const56 E_F0* code(const basicAC_F0& args) const {
57 return new init_Op<LR, K>(args, t[0]->CastTo(args[0]));
58 }
59 };
60
61 template<template<class> class LR, class K>
62 class assembleHMatrix : public OneOperator { public:
63
64 class Op : public E_F0info {
65 public:
66 Expression a,b;
67
68 static const int n_name_param = 7;
69 static basicAC_F0::name_and_type name_param[] ;
70 Expression nargs[n_name_param];
arg(int i,Stack stack,bool a) const71 bool arg(int i,Stack stack,bool a) const{ return nargs[i] ? GetAny<bool>( (*nargs[i])(stack) ): a;}
argl(int i,Stack stack,long a) const72 long argl(int i,Stack stack,long a) const{ return nargs[i] ? GetAny<long>( (*nargs[i])(stack) ): a;}
arg(int i,Stack stack,double a) const73 double arg(int i,Stack stack,double a) const{ return nargs[i] ? GetAny<double>( (*nargs[i])(stack) ): a;}
arg(int i,Stack stack,KN_<long> a) const74 KN_<long> arg(int i,Stack stack,KN_<long> a ) const{ return nargs[i] ? GetAny<KN_<long> >( (*nargs[i])(stack) ): a;}
arg(int i,Stack stack,pcommworld a) const75 pcommworld arg(int i,Stack stack,pcommworld a ) const{ return nargs[i] ? GetAny<pcommworld>( (*nargs[i])(stack) ): a;}
76 public:
Op(const basicAC_F0 & args,Expression aa,Expression bb)77 Op(const basicAC_F0 & args,Expression aa,Expression bb) : a(aa),b(bb) {
78 args.SetNameParam(n_name_param,name_param,nargs); }
79 };
80
assembleHMatrix()81 assembleHMatrix() : OneOperator(atype<const typename assembleHMatrix<LR,K>::Op *>(),
82 atype<pfes *>(),
83 atype<pfes *>()) {}
84
code(const basicAC_F0 & args) const85 E_F0 * code(const basicAC_F0 & args) const
86 {
87 return new Op(args,t[0]->CastTo(args[0]),
88 t[1]->CastTo(args[1]));
89 }
90 };
91
92 template<template<class> class LR, class K>
93 basicAC_F0::name_and_type assembleHMatrix<LR, K>::Op::name_param[]= {
94 { "epsilon", &typeid(double)},
95 { "eta", &typeid(double)},
96 { "minclustersize", &typeid(long)},
97 { "maxblocksize", &typeid(long)},
98 { "mintargetdepth", &typeid(long)},
99 { "minsourcedepth", &typeid(long)},
100 { "comm", &typeid(pcommworld)}
101 };
102
103 /*
104 template <template<class> class LR, class K, EquationEnum Eq, BIOpKernelEnum Op, BIOpKernelEnum local_Op, int dim, typename Discretization>
105 HMatrix<LR,K>* generateBIO() {
106 double kappa = 1;
107 Geometry node("Th.msh");
108 bemtool::Mesh<dim> mesh; mesh.Load(node,0);
109 Orienting(mesh);
110 int nb_elt = NbElt(mesh);
111 std::vector<bemtool::R3> normals= NormalTo(mesh);
112
113 // Dof
114 Dof<Discretization> dof(mesh);
115 int nb_dof = NbDof(dof);
116 std::vector<htool::R3> x(nb_dof);
117 for (int i=0;i<nb_dof;i++){
118 x[i][0]=dof(((dof.ToElt(i))[0])[0])[((dof.ToElt(i))[0])[1]][0];
119 x[i][1]=dof(((dof.ToElt(i))[0])[0])[((dof.ToElt(i))[0])[1]][1];
120 x[i][2]=dof(((dof.ToElt(i))[0])[0])[((dof.ToElt(i))[0])[1]][2];
121 }
122
123 BIO_Generator<BIOpKernel<Eq,Op,dim+1,Discretization,Discretization>,Discretization> generator_BIO(dof,kappa);
124 HMatrix<LR,K>* H = new HMatrix<LR,K>(generator_BIO,x);
125
126
127 // if (apply_boundary(Op)){
128 // std::vector<int> boundary=is_boundary_nodes(dof);
129 // if (*max_element(boundary.begin(),boundary.end())!=0){
130 // H->apply_dirichlet(boundary);
131 // }
132 // }
133
134 return H;
135 }
136 */
137
138 template<template<class> class LR, class K>
SetHMatrix(Stack stack,Expression emat,Expression einter,int init)139 AnyType SetHMatrix(Stack stack,Expression emat,Expression einter,int init)
140 {
141 using namespace Fem2D;
142
143 HMatrix<LR,K>** Hmat =GetAny<HMatrix<LR,K>** >((*emat)(stack));
144 const typename assembleHMatrix<LR, K>::Op * mi(dynamic_cast<const typename assembleHMatrix<LR, K>::Op *>(einter));
145
146 double epsilon=mi->arg(0,stack,htool::Parametres::epsilon);
147 double eta=mi->arg(1,stack,htool::Parametres::eta);
148 int minclustersize=mi->argl(2,stack,htool::Parametres::minclustersize);
149 int maxblocksize=mi->argl(3,stack,htool::Parametres::eta);
150 int mintargetdepth=mi->argl(4,stack,htool::Parametres::mintargetdepth);
151 int minsourcedepth=mi->argl(5,stack,htool::Parametres::minsourcedepth);
152 pcommworld pcomm=mi->arg(6,stack,nullptr);
153
154 MPI_Comm comm = pcomm ? *(MPI_Comm*)pcomm : MPI_COMM_WORLD;
155
156 ffassert(einter);
157 pfes * pUh = GetAny< pfes * >((* mi->a)(stack));
158 FESpace * Uh = **pUh;
159 int NUh =Uh->N;
160
161 pfes * pVh = GetAny< pfes * >((* mi->b)(stack));
162 FESpace * Vh = **pVh;
163 int NVh =Vh->N;
164
165 ffassert(Vh);
166 ffassert(Uh);
167
168 int n=Uh->NbOfDF;
169 int m=Vh->NbOfDF;
170
171 const Fem2D::Mesh & ThU =Uh->Th; // line
172 const Fem2D::Mesh & ThV =Vh->Th; // colunm
173 bool samemesh = &Uh->Th == &Vh->Th; // same Fem2D::Mesh
174
175 //cout << samemesh << " " << NUh << " " << n << " " << m << endl;
176
177 vector<htool::R3> p1(Uh->Th.nv);
178 vector<htool::R3> p2(Vh->Th.nv);
179 for (int i=0; i<ThU.nv; i++)
180 p1[i] = {ThU.vertices[i].x, ThU.vertices[i].y, 0};
181
182 for (int i=0; i<ThV.nv; i++)
183 p2[i] = {ThV.vertices[i].x, ThV.vertices[i].y, 0};
184
185
186 MyMatrix<K> A(Uh,Vh);
187
188 SetMaxBlockSize(maxblocksize);
189 SetMinClusterSize(minclustersize);
190 SetEpsilon(epsilon);
191 SetEta(eta);
192 SetMinTargetDepth(mintargetdepth);
193 SetMinSourceDepth(minsourcedepth);
194
195 if (init)
196 delete *Hmat;
197 *Hmat = new HMatrix<LR,K>(A,p1,p2,-1,comm);
198
199 //*Hmat = generateBIO<LR,K,EquationEnum::YU,BIOpKernelEnum::HS_OP,BIOpKernelEnum::SL_OP,2,P1_2D>();
200
201 return Hmat;
202 }
203
204 template<template<class> class LR, class K, int init>
SetHMatrix(Stack stack,Expression emat,Expression einter)205 AnyType SetHMatrix(Stack stack,Expression emat,Expression einter)
206 { return SetHMatrix<LR,K>(stack,emat,einter,init);}
207
208 template<class R, class A, class B> R Build(A a, B b) {
209 return R(a, b);
210 }
211
212 template<template<class> class LR, class K>
ToDense(Stack stack,Expression emat,Expression einter,int init)213 AnyType ToDense(Stack stack,Expression emat,Expression einter,int init)
214 {
215 ffassert(einter);
216 HMatrix<LR,K>** Hmat =GetAny<HMatrix<LR,K>** >((*einter)(stack));
217 ffassert(Hmat && *Hmat);
218 HMatrix<LR,K>& H = **Hmat;
219 Matrix<K> mdense = H.to_dense_perm();
220 const std::vector<K>& vdense = mdense.get_mat();
221
222 KNM<K>* M =GetAny<KNM<K>*>((*emat)(stack));
223
224 for (int i=0; i< mdense.nb_rows(); i++)
225 for (int j=0; j< mdense.nb_cols(); j++)
226 (*M)(i,j) = mdense(i,j);
227
228 return M;
229 }
230
231 template<template<class> class LR, class K, int init>
ToDense(Stack stack,Expression emat,Expression einter)232 AnyType ToDense(Stack stack,Expression emat,Expression einter)
233 { return ToDense<LR,K>(stack,emat,einter,init);}
234
235 template<class V, template<class> class LR, class K>
236 class Prod {
237 public:
238 const HMatrix<LR ,K>* h;
239 const V u;
Prod(HMatrix<LR,K> ** v,V w)240 Prod(HMatrix<LR ,K>** v, V w) : h(*v), u(w) {}
241
prod(V x) const242 void prod(V x) const {h->mvprod_global(*(this->u), *x);};
243
mv(V Ax,Prod<V,LR,K> A)244 static V mv(V Ax, Prod<V, LR, K> A) {
245 *Ax = K();
246 A.prod(Ax);
247 return Ax;
248 }
init(V Ax,Prod<V,LR,K> A)249 static V init(V Ax, Prod<V, LR, K> A) {
250 Ax->init(A.u->n);
251 return mv(Ax, A);
252 }
253
254 };
255
256 template<template<class> class LR, class K>
get_infos(HMatrix<LR,K> ** const & H)257 std::map<std::string, std::string>* get_infos(HMatrix<LR ,K>** const& H) {
258 return new std::map<std::string, std::string>((*H)->get_infos());
259 }
260
get_info(std::map<std::string,std::string> * const & infos,string * const & s)261 string* get_info(std::map<std::string, std::string>* const& infos, string* const& s){
262 return new string((*infos)[*s]);
263 }
264
operator <<(ostream & out,const std::map<std::string,std::string> & infos)265 ostream & operator << (ostream &out, const std::map<std::string, std::string> & infos)
266 {
267 for (std::map<std::string,std::string>::const_iterator it = infos.begin() ; it != infos.end() ; ++it){
268 out<<it->first<<"\t"<<it->second<<std::endl;
269 }
270 out << std::endl;
271 return out;
272 }
273
274 template<class A>
275 struct PrintPinfos: public binary_function<ostream*,A,ostream*> {
fPrintPinfos276 static ostream* f(ostream* const & a,const A & b) { *a << *b;
277 return a;
278 }
279 };
280
281 template<template<class> class LR, class K>
282 class plotHMatrix : public OneOperator {
283 public:
284
285 class Op : public E_F0info {
286 public:
287 Expression a;
288
289 static const int n_name_param = 1;
290 static basicAC_F0::name_and_type name_param[] ;
291 Expression nargs[n_name_param];
arg(int i,Stack stack,bool a) const292 bool arg(int i,Stack stack,bool a) const{ return nargs[i] ? GetAny<bool>( (*nargs[i])(stack) ): a;}
arg(int i,Stack stack,long a) const293 long arg(int i,Stack stack,long a) const{ return nargs[i] ? GetAny<long>( (*nargs[i])(stack) ): a;}
294
295 public:
Op(const basicAC_F0 & args,Expression aa)296 Op(const basicAC_F0 & args,Expression aa) : a(aa) {
297 args.SetNameParam(n_name_param,name_param,nargs);
298 }
299
operator ()(Stack stack) const300 AnyType operator()(Stack stack) const{
301
302 bool wait = arg(0,stack,false);
303
304 HMatrix<LR,K>** H =GetAny<HMatrix<LR,K>** >((*a)(stack));
305
306 PlotStream theplot(ThePlotStream);
307
308 if (mpirank == 0 && ThePlotStream) {
309 theplot.SendNewPlot();
310 theplot << 3L;
311 theplot <= wait;
312 theplot.SendEndArgPlot();
313 theplot.SendMeshes();
314 theplot << 0L;
315 theplot.SendPlots();
316 theplot << 1L;
317 theplot << 31L;
318 }
319
320 if (!H || !(*H)) {
321 if (mpirank == 0&& ThePlotStream) {
322 theplot << 0;
323 theplot << 0;
324 theplot << 0L;
325 theplot << 0L;
326 }
327 }
328 else {
329 const std::vector<LR<K>*>& lrmats = (*H)->get_MyFarFieldMats();
330 const std::vector<SubMatrix<K>*>& dmats = (*H)->get_MyNearFieldMats();
331
332 int nbdense = dmats.size();
333 int nblr = lrmats.size();
334
335 int sizeworld = (*H)->get_sizeworld();
336 int rankworld = (*H)->get_rankworld();
337
338 int nbdenseworld[sizeworld];
339 int nblrworld[sizeworld];
340 MPI_Allgather(&nbdense, 1, MPI_INT, nbdenseworld, 1, MPI_INT, (*H)->get_comm());
341 MPI_Allgather(&nblr, 1, MPI_INT, nblrworld, 1, MPI_INT, (*H)->get_comm());
342 int nbdenseg = 0;
343 int nblrg = 0;
344 for (int i=0; i<sizeworld; i++) {
345 nbdenseg += nbdenseworld[i];
346 nblrg += nblrworld[i];
347 }
348
349 int* buf = new int[4*(mpirank==0?nbdenseg:nbdense) + 5*(mpirank==0?nblrg:nblr)];
350
351 for (int i=0;i<dmats.size();i++) {
352 const SubMatrix<K>& l = *(dmats[i]);
353 buf[4*i] = l.get_offset_i();
354 buf[4*i+1] = l.get_offset_j();
355 buf[4*i+2] = l.nb_rows();
356 buf[4*i+3] = l.nb_cols();
357 }
358
359 int displs[sizeworld];
360 int recvcounts[sizeworld];
361 displs[0] = 0;
362
363 for (int i=0; i<sizeworld; i++) {
364 recvcounts[i] = 4*nbdenseworld[i];
365 if (i > 0) displs[i] = displs[i-1] + recvcounts[i-1];
366 }
367 MPI_Gatherv(rankworld==0?MPI_IN_PLACE:buf, recvcounts[rankworld], MPI_INT, buf, recvcounts, displs, MPI_INT, 0, (*H)->get_comm());
368
369 int* buflr = buf + 4*(mpirank==0?nbdenseg:nbdense);
370 double* bufcomp = new double[mpirank==0?nblrg:nblr];
371
372 for (int i=0;i<lrmats.size();i++) {
373 const LR<K>& l = *(lrmats[i]);
374 buflr[5*i] = l.get_offset_i();
375 buflr[5*i+1] = l.get_offset_j();
376 buflr[5*i+2] = l.nb_rows();
377 buflr[5*i+3] = l.nb_cols();
378 buflr[5*i+4] = l.rank_of();
379 bufcomp[i] = l.compression();
380 }
381
382 for (int i=0; i<sizeworld; i++) {
383 recvcounts[i] = 5*nblrworld[i];
384 if (i > 0) displs[i] = displs[i-1] + recvcounts[i-1];
385 }
386
387 MPI_Gatherv(rankworld==0?MPI_IN_PLACE:buflr, recvcounts[rankworld], MPI_INT, buflr, recvcounts, displs, MPI_INT, 0, (*H)->get_comm());
388
389 for (int i=0; i<sizeworld; i++) {
390 recvcounts[i] = nblrworld[i];
391 if (i > 0) displs[i] = displs[i-1] + recvcounts[i-1];
392 }
393
394 MPI_Gatherv(rankworld==0?MPI_IN_PLACE:bufcomp, recvcounts[rankworld], MPI_DOUBLE, bufcomp, recvcounts, displs, MPI_DOUBLE, 0, (*H)->get_comm());
395
396 if (mpirank == 0 && ThePlotStream ) {
397
398 int si = (*H)->nb_rows();
399 int sj = (*H)->nb_cols();
400
401 theplot << si;
402 theplot << sj;
403 theplot << (long)nbdenseg;
404 theplot << (long)nblrg;
405
406 for (int i=0;i<nbdenseg;i++) {
407 theplot << buf[4*i];
408 theplot << buf[4*i+1];
409 theplot << buf[4*i+2];
410 theplot << buf[4*i+3];
411 }
412
413 for (int i=0;i<nblrg;i++) {
414 theplot << buflr[5*i];
415 theplot << buflr[5*i+1];
416 theplot << buflr[5*i+2];
417 theplot << buflr[5*i+3];
418 theplot << buflr[5*i+4];
419 theplot << bufcomp[i];
420 }
421
422 theplot.SendEndPlot();
423
424 }
425 delete [] buf;
426 delete [] bufcomp;
427
428 }
429
430 return 0L;
431 }
432 };
433
plotHMatrix()434 plotHMatrix() : OneOperator(atype<long>(),atype<HMatrix<LR ,K> **>()) {}
435
code(const basicAC_F0 & args) const436 E_F0 * code(const basicAC_F0 & args) const
437 {
438 return new Op(args,t[0]->CastTo(args[0]));
439 }
440 };
441
442 template<template<class> class LR, class K>
443 basicAC_F0::name_and_type plotHMatrix<LR, K>::Op::name_param[]= {
444 { "wait", &typeid(bool)}
445 };
446
447 template<template<class> class LR, class K>
add(const char * namec)448 void add(const char* namec) {
449 Dcl_Type<HMatrix<LR ,K>**>(Initialize<HMatrix<LR,K>*>, Delete<HMatrix<LR,K>*>);
450 Dcl_TypeandPtr<HMatrix<LR ,K>*>(0,0,::InitializePtr<HMatrix<LR ,K>*>,::DeletePtr<HMatrix<LR ,K>*>);
451 //TheOperators->Add("<-", new init<LR,K>);
452 Dcl_Type<const typename assembleHMatrix<LR, K>::Op *>();
453 Add<const typename assembleHMatrix<LR, K>::Op *>("<-","(", new assembleHMatrix<LR, K>);
454
455 TheOperators->Add("=",
456 new OneOperator2_<HMatrix<LR ,K>**,HMatrix<LR ,K>**,const typename assembleHMatrix<LR, K>::Op*,E_F_StackF0F0>(SetHMatrix<LR, K, 1>));
457 TheOperators->Add("<-",
458 new OneOperator2_<HMatrix<LR ,K>**,HMatrix<LR ,K>**,const typename assembleHMatrix<LR, K>::Op*,E_F_StackF0F0>(SetHMatrix<LR, K, 0>));
459 Global.Add(namec,"(",new assembleHMatrix<LR, K>);
460
461 //atype<HMatrix<LR ,K>**>()->Add("(","",new OneOperator2_<string*, HMatrix<LR ,K>**, string*>(get_infos<LR,K>));
462 Add<HMatrix<LR ,K>**>("infos",".",new OneOperator1_<std::map<std::string, std::string>*, HMatrix<LR ,K>**>(get_infos));
463
464 Dcl_Type<Prod<KN<K>*, LR, K>>();
465 TheOperators->Add("*", new OneOperator2<Prod<KN<K>*, LR, K>, HMatrix<LR ,K>**, KN<K>*>(Build));
466 TheOperators->Add("=", new OneOperator2<KN<K>*, KN<K>*, Prod<KN<K>*, LR, K>>(Prod<KN<K>*, LR, K>::mv));
467 TheOperators->Add("<-", new OneOperator2<KN<K>*, KN<K>*, Prod<KN<K>*, LR, K>>(Prod<KN<K>*, LR, K>::init));
468
469 Global.Add("display","(",new plotHMatrix<LR, K>);
470
471 // to dense:
472 TheOperators->Add("=",
473 new OneOperator2_<KNM<K>*, KNM<K>*, HMatrix<LR ,K>**,E_F_StackF0F0>(ToDense<LR, K, 1>));
474 TheOperators->Add("<-",
475 new OneOperator2_<KNM<K>*, KNM<K>*, HMatrix<LR ,K>**,E_F_StackF0F0>(ToDense<LR, K, 0>));
476 }
477
Init_Schwarz()478 static void Init_Schwarz() {
479 Dcl_Type<std::map<std::string, std::string>*>( );
480 TheOperators->Add("<<",new OneBinaryOperator<PrintPinfos<std::map<std::string, std::string>*>>);
481 Add<std::map<std::string, std::string>*>("[","",new OneOperator2_<string*, std::map<std::string, std::string>*, string*>(get_info));
482
483 add<partialACA,double>("assemble");
484 add<partialACA,std::complex<double>>("assemblecomplex");
485
486 zzzfff->Add("HMatrix", atype<HMatrix<partialACA ,double>**>());
487 map_type_of_map[make_pair(atype<HMatrix<partialACA ,double>**>(), atype<double*>())] = atype<HMatrix<partialACA ,double>**>();
488 map_type_of_map[make_pair(atype<HMatrix<partialACA ,double>**>(), atype<Complex*>())] = atype<HMatrix<partialACA ,std::complex<double>>**>();
489 }
490
491 LOADFUNC(Init_Schwarz)
492