1NewMacro DMMDeffuncAndGlobals(prefix,comm,jpart,Whi,Vhc,Pknbcomp,Ai,vPbC,onG10,Pii,Usend,Vrecv,U2V)
2/*
3
4*/
5
6int prefix#rgmres=0;
7int prefix#kiter=-1;
8/*  the global name user... */
9/*
10  jpart,njpart : partition
11  Usend, Vrecv : buffer
12  Ai , Bi
13  rMj, sMj : matrices
14
15  onG10:  tgv only on DDM s Gamma and not on Gamma
16*/
17/*----  for coarse  solver ... */
18
19
20matrix prefix#AC,prefix#Rci,prefix#Pci;/**/
21{
22int[int] Sigma11=U2V;
23prefix#Pci=   interpolate(Whi,VhC,U2Vc=Sigma11);
24prefix#Rci =  prefix#Pci'*Pii;
25/*Rci=   interpolate(Whi,VhC,t=1,U2Vc=Sigma11(0:Pknbcomp-1));
26Pci =  Pii*Rci'; */
27}
28
29/*----End of Global Def -------------*/
30
31/*-----------------*/
32/*-----------------*/
33
34
35func bool prefix#Update(real[int] &ui, real[int] &vi)
36{
37  for(int j=0;j<jpart.n;++j)
38    Usend[j][]=sMj[j]*ui;
39   SendRecvUV(comm,jpart,Usend,Vrecv)
40     vi = Pii*ui;
41   for(int j=0;j<jpart.n;++j)
42     vi += rMj[j]*Vrecv[j][];
43   return true;
44}
45/*-----------------*/
46/*-----------------*/
47/*-----------------*/
48
49func bool  prefix#CoarseSolve(real[int]& V,real[int]& U,mpiComm& comm)
50{
51    if(prefix#AC.n==0 && mpiRank(comm)==0)
52      prefix#AC = vPbC(VhC,VhC,solver=sparsesolver);
53   /*  solvibg the coarse probleme */
54   real[int] Uc(prefix#Rci.n),Bc(Uc.n);
55   Uc= prefix#Rci*U;
56   mpiReduce(Uc,Bc,processor(0,comm),mpiSUM);
57   if(mpiRank(comm)==0)
58      Uc = prefix#AC^-1*Bc;
59    broadcast(processor(0,comm),Uc);
60   V = prefix#Pci*Uc;
61}
62
63/*-----------------*/
64/*-----------------*/
65/*-----------------*/
66
67func real[int] prefix#DJ(real[int]& U)
68{
69  ++prefix#kiter;
70  real[int] V(U.n);
71  V =  Ai*U;
72  V = onG10 ? 0.: V;
73  return V;
74}
75
76func real[int] prefix#PDJ(real[int]& U) /* C1*/
77{
78  real[int] V(U.n);
79
80  real[int] b= onG10 ? 0. :  U;
81  V =  Ai^-1*b;
82  prefix#Update(V,U);
83  return U;
84}
85
86func real[int] prefix#PDJC(real[int]& U) /**/
87{ /* Precon  C1= Precon , C2  precon Coarse
88   Idea : F. Nataf.
89     0 ~  (I-C1A)(I-C2A) => I ~  - C1AC2A +C1A +C2A
90     New Prec P= C1+C2 - C1AC2   = C1(I- A C2) +C2
91   (  C1(I- A C2) +C2 ) Uo
92     V =  - C2*Uo
93  .... */
94  real[int] V(U.n);
95  prefix#CoarseSolve(V,U,comm);
96  V = -V; /*  -C2*Uo  */
97  U  += Ai*V; /* U =  (I-A C2) Uo */
98  real[int] b= onG10 ? 0. :  U;
99  U =  Ai^-1*b;	/*  ( C1( I -A C2) Uo */
100  V = U -V; /**/
101  prefix#Update(V,U);
102  return U;
103}
104
105/*-----------------*/
106/*-----------------*/
107/*-----------------*/
108
109func real[int] prefix#PDJC2(real[int]& U) /*  bogus ???? */
110{ /* Precon  C1= precon Coarse C2  precon Precon
111   Idea : F. Nataf.
112      0 ~  (I C1A)(I-C2A) => I ~  - C1AC2A +C1A +C2A
113      New Prec P= C1+C2 - C1AC2   = C1(I- A C2) +C2
114      (  C1(I- A C2) +C2 ) Uo
115      V =  - C2*Uo
116      ....
117      V = - C2 Uo
118      W = Uo + A V
119      V + C1 W
120      */
121  real[int] V(U.n);
122  real[int] b= onG10 ? 0. :  U;
123  V =  Ai^-1*b;
124  b=U;
125
126  V = -V;
127  prefix#Update(V,U);
128  U  += Ai*V;
129
130  prefix#CoarseSolve(U,b,comm);
131  V = U -V;
132  prefix#Update(V,U);
133  return U;
134}
135/*-----------------*/
136/*-----------------*/
137/*-----------------*/
138
139 func real[int] prefix#DJ0(real[int]& U)
140{
141  ++prefix#kiter;
142  real[int] V(U.n);
143  real[int] b= onG10 .* U;
144  b  = onG10 ? b : Bi ;
145  V = Ai^-1*b;
146  prefix#Update(V,U);
147  V -= U;
148   return V;
149}
150
151/*-----------------*/
152/*-----------------*/
153/*-----------------*/
154
155func bool  prefix#CheckUpdate()
156{ /* verification.....*/
157
158  Whi defPk#Pknbcomp(u,) =Times#Pknbcomp(1),defPk#Pknbcomp(v,);
159  prefix#Update(u[],v[]);
160  u[]-=v[];
161  if(mpirank==0 || (u[].linfty>1e-6))
162     cout << "CheckUpdate  " << u[].linfty << " rank: " << mpirank <<endl;
163  if(u[].linfty>1e-6) plot(u,cmm="bug ????");
164  assert( u[].linfty<1e-6);
165
166  return 1; }
167
168func bool prefix#Saveff(string sff,real epss,int ksplit,int nloc,int sizeoverlaps)
169{  if(sff != "")
170  {
171    ofstream ff(sff+".txt",append);
172    cout << " ++++  " ;
173    cout  << mpirank <<"/" <<  mpisize << " k=" <<  ksplit << " n= " << nloc << " "
174           << sizeoverlaps << " it=  " << prefix#kiter  ;
175    for (int i=1; i<ittt;++i)
176      cout << " " << ttt[i]-ttt[i-1] << " ";
177    cout << epss << " " << Ai.nbcoef << " " << Ai.n << endl;
178
179    /*
180      1 mpirank
181      2 mpisize
182      3 ksplit
183      4 nloc
184      5 sizeoverlaps
185      6 kiter
186      7 mesh & part build
187      8 build the partion
188      9 build mesh, transfere , and the fine mesh ..
189      10 build the matrix,  the trans matrix, factorizatioon
190      11 GMRES
191    */
192    ff   << mpirank << " " << mpisize << " " << sPk << " " ;
193    ff <<  ksplit << " " << nloc << " " << sizeoverlaps << " " << prefix#kiter  ;
194    for (int i=1; i<ittt;++i)
195      ff << " " << ttt[i]-ttt[i-1] << " ";
196    ff << epss << " " << Ai.nbcoef << " " << Ai.n << " " << gmres << endl;
197
198  }
199  return 1;
200}
201//
202/*-----------------*/
203/*-----------------*/
204/*-----------------*/
205/*-----------------*/
206/*-----------------*/
207
208
209NewMacro prefix#DDMSolver(Bi,u,v,gmres,epss,vdebug)
210{
211 settt
212 int ipart=mpiRank(comm);
213 if(gmres==1)
214  {
215   prefix#rgmres=MPIAffineGMRES(prefix#DJ0,u[],veps=epss,nbiter=300,comm=comm,dimKrylov=100,verbosity=ipart ? 0: 50);
216   real[int] b= onG10 .* u[];
217   b  = onG10 ? b : Bi ;
218   v[] = Ai^-1*b;
219   prefix#Update(v[],u[]);
220  }
221else if(gmres==2)
222  prefix#rgmres= MPILinearGMRES(prefix#DJ,precon=prefix#PDJ,u[],Bi,veps=epss,nbiter=1000,comm=comm,
223                         dimKrylov=100,verbosity=ipart ? 0: 50);
224else if(gmres==3)
225  prefix#rgmres= MPILinearGMRES(prefix#DJ,precon=prefix#PDJC,u[],Bi,veps=epss,nbiter=1000,comm=comm,
226                          dimKrylov=100,verbosity=ipart ? 0: 50);
227else if(gmres==4)
228   prefix#rgmres= MPILinearGMRES(prefix#DJ,precon=prefix#PDJC2,u[],Bi,veps=epss,nbiter=1000,comm=comm,
229                          dimKrylov=100,verbosity=ipart ? 0: 50);
230else /*algo Shwarz for demo */
231   for(int iter=0;iter <100; ++iter)
232     {
233       prefix#kiter=iter;
234       real[int] b= onG10 .* u[];
235       b  = onG10 ? b : Bi ;
236       v[] = Ai^-1*b;
237       b=v[];
238       prefix#Update(v[],u[]);
239       if(vdebug) plotMPIall(Thi,u[],"u-"+iter);
240        b -= u[];
241
242       real err = b.l1;
243       real umax = u[].max;
244       real[int] aa=[err,umax], bb(2);
245       mpiAllReduce(aa,bb,comm,mpiMAX);
246       real errg = bb[0];
247       real umaxg = bb[1];
248
249       if(ipart==0)
250	     cout << ipart << " " << iter << " err = " << errg << " u. max  " << umaxg << endl;
251       if(errg< 1e-10 ) break;
252     }
253     prefix#Update(u[],v[]);
254
255if(vdebug) plotMPIall(Thi,v[],"u-final");
256 u[]=v[];
257 settt
258}
259EndMacro  /*prefix#DDMSolver*/
260
261EndMacro
262
263
264