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