
// ff-mpirun -np 4 MPIGMRES2D.edp -glut ffglut  -n 11 -k 1  -d 1 -ns -gmres 1
/*
  a first true parallele example fisrt freefem++ 
  Ok up to 200 proc for a Poisson equation.. 
  See the Doc for full explaiantion

  F Hecht Dec. 2010. 
  -------------------
usage :
ff-mpirun [mpi parameter] MPIGMRES2d.edp  [-glut ffglut]  [-n N] [-k K]  [-d D] [-ns] [-gmres [0|1]
 argument: 
   -glut ffglut : to see graphicaly the process
   -n N:  set the mesh cube split NxNxN
   -d D:  set debug flag D must be one for mpiplot 
   -k K:  to refined by K all  elemnt
   -ns: reomove script dump
   -gmres 0   : use iterative schwarz algo.  
          1   :  Algo GMRES on residu of schwarz algo.
          2   :  DDM GMRES 
          3   :  DDM GMRES with coarse grid preconditionner (Good one)  
*/
load "MPICG"  load "medit"  load "metis"
include "getARGV.idp"
include "MPIplot.idp"
include "MPIGMRESmacro.idp"
//include "AddLayer2d.idp" 

searchMethod=0; // more safe seach algo (warning can be very expensive in case lot of ouside point) 
assert(version >=3.11);
real[int] ttt(10);int ittt=0;
macro settt {ttt[ittt++]=mpiWtime();}//


verbosity=getARGV("-vv",0);
int vdebug=getARGV("-d",1);
int ksplit=getARGV("-k",3);
int nloc = getARGV("-n",11);
string sff=getARGV("-p,","");
int gmres=getARGV("-gmres",2); 
bool dplot=getARGV("-dp",0);
int nC = getARGV("-N" ,max(nloc/10,4)); 

if(mpirank==0 && verbosity)
{
  cout << "ARGV : ";
  for(int i=0;i<ARGV.n;++i)
    cout << ARGV[i] <<" ";
  cout << endl;
}


if(mpirank==0 && verbosity)
  cout << " vdebug: " << vdebug << " kspilt "<< ksplit << " nloc "<< nloc << " sff "<< sff <<"."<< endl;

int withplot=0;
bool withmetis=1;
bool RAS=1; 
string sPk="P2-2gd";     
func Pk=P2;

func bool  plotMPIall(mesh &Th,real[int] & u,string  cm)
{if(vdebug) PLOTMPIALL(mesh,Pk, Th, u,{ cmm=cm,nbiso=20,fill=1,dim=3,value=1}); return 1;}

int sizeoverlaps=1; // size of overlap

mpiComm comm(mpiCommWorld,0,0);// trick : make a no split mpiWorld 
 
int npart=mpiSize(comm); // total number of partion
int ipart= mpiRank(comm); // current partition number 

int njpart=0; // nb of part with intersection (a jpart) with ipart without ipart
int[int] jpart(npart); //  list of jpart ..
if(ipart==0)  cout << " Final N=" << ksplit*nloc << " nloc =" << nloc << " split =" << ksplit <<  endl;
int[int] l111=[1,1,1,1]; 
settt 

mesh Thg=square(nloc,nloc,label=l111);
mesh ThC=square(nC,nC,label=l111);//   Caarse Mesh

mesh Thi,Thin;//  with overlap, without olverlap  
fespace Phg(Thg,P0);
fespace Vhg(Thg,P1);
fespace VhC(ThC,P1); // of the coarse problem.. 


Phg  part;

// build the partitioning ...  
{    
 int[int] nupart(Thg.nt);
 nupart=0; 
 if(npart>1 && ipart==0)
   metisdual(nupart,Thg,npart); 
   
 broadcast(processor(0,comm),nupart);
 for(int i=0;i<nupart.n;++i)
    part[][i]=nupart[i];
  
} // build ... 


if(withplot>1)
  plot(part,fill=1,cmm="dual",wait=1);

// overlapping partition

 Phg suppi= abs(part-ipart)<0.1; 
 Vhg unssd;                       // boolean function 1 in the subdomain 0 elswhere 
 Thin=trunc(Thg,suppi>0,label=10); // non-overlapping mesh, interfaces have label 10
 int nnn = sizeoverlaps*2;// to be sure 
 AddLayers(Thg,suppi[],nnn,unssd[]);    // see above ! suppi and unssd are modified  
 unssd[] *= nnn;  //  to put value nnn a 0  
 real nnn0 = nnn - sizeoverlaps +  0.001   ;
 Thi=trunc(Thg,unssd>nnn0 ,label=10); // overlapping mesh, interfaces have label 10
 
 settt 

 fespace Vhi(Thi,P1);
 int npij=npart;
 Vhi[int] pij(npij);// local partition of unit + pii 
 Vhi pii;  
 
 real nnn1=  + 0.001  ;
 { /*
   construction of the partition of the unit,
    let phi_i P1 FE function 1 on Thin and zero ouside of Thi and positive
    the partition is build with  
  $$  p_i = phi_i/ \sum phi_i 
    
    to build the partition of one domain i
    we nned to find all j such that supp(phi_j) \cap supp(phi_j) is not empty
    <=> int phi_j
 */
 //   build a local mesh of thii such that all compuation of the unit partition are
 //   exact in thii
 mesh Thii=trunc(Thg,unssd>nnn1 ,label=10); // overlapping mesh, interfaces have label 10


 {  
   // find all j  mes (supp(p_j) \cap supp(p_i)) >0  
   // compute all phi_j on Thii
   //  remark supp p_i include in Thi 
   //  
   fespace Phii(Thii,P0);
   fespace Vhii(Thii,P1);
   Vhi sumphi=0;
   jpart=0;
   njpart=0;
   int nlayer=RAS?1:sizeoverlaps;
   if(ipart==0) 
     cout <<" nlayer=" << nlayer << endl;
   pii= max(unssd-nnn+nlayer,0.)/nlayer;
   if(dplot) plot(pii,wait=1,cmm=" 0000");
   sumphi[] +=  pii[];
   if(dplot) plot(sumphi,wait=1,cmm=" summ 0000");
   Vhii phii=0;
   real epsmes=1e-10*Thii.area;
   for (int i=0;i<npart;++i)
     if(i != ipart ) 
       { 
	    Phii suppii=abs(i-part)<0.2; 
	    if(suppii[].max > 0.5) 
	    {
	     AddLayers(Thii,suppii[],nlayer,phii[]);
	     assert(phii[].min >=0);
	     real interij = int2d(Thi)(  phii); 
	     if(interij>epsmes) 
	       {  
		     pij[njpart]=abs(phii);	 
		     if(vdebug>2) cout << " ***** " << int2d(Thi)(real(pij[njpart])<0) << " " <<pij[njpart][].min << " " << phii[].min << endl;
		     assert(int2d(Thi)(real(pij[njpart])<0) ==0);
		     if(dplot)  plot(pij[njpart],wait=1,cmm=" j = "+ i + " " + njpart);
		     sumphi[] += pij[njpart][]; 
		     if(dplot)  plot(sumphi,wait=1,cmm=" sum j = "+ i + " " + njpart);
		     jpart[njpart++]=i;
	       }}}
   
   if(dplot) plot(sumphi,wait=1,dim=3,cmm="sum ",fill=1 );
   pii[]=pii[] ./ sumphi[];
   for (int j=0;j<njpart;++j)
     pij[j][] = pij[j][] ./ sumphi[];
   jpart.resize(njpart);	
   for (int j=0;j<njpart;++j)
     assert(pij[j][].max<=1);
   {
     cout << ipart << " number of jpart " << njpart << " : ";
     for (int j=0;j<njpart;++j)
       cout << jpart[j] << " ";
  	cout << endl;
   }
   sumphi[]=pii[];
   for (int j=0;j<njpart;++j)
     sumphi[]+= pij[j][];
   if(vdebug>2)  
     cout << " sum min " <<sumphi[].min << " " << sumphi[].max << endl;
   assert(sumphi[].min> 1.-1e-6 && sumphi[].max< 1.+1e-6);  
   //  verification 
 }}// (Thii is remove here) 
  // end of the construction of the local partition of the unity ...
  // on Thi ...  
  // -----------------------------------------------------------------
if(ipart==0) cout << " *** end build partition " << endl;

//  computation of  number of intersection ..
// ------------------------------------------

// here  pii and the pij is the locate partition of the unite on 
// Thi ( mesh with overlap )....
//Thi=splitmesh(Thi,2);
if(dplot ) 
  { plot(Thi,wait=1); 
    for(int j=0;j<njpart;++j)
      plot(pij[j],cmm=" j="+j ,wait=1); } 
      
//  Partition of the unity on Thi ..
// computation of message.
// all j> we have to recive 
// data on intersection of the support of pij[0] and pij[j]
settt 

 if(vdebug) plotMPIall(Thi,pii[],"pi_i");

mesh[int] aThij(njpart); 
matrix Pii;  
matrix[int] sMj(njpart); // M of send to j
matrix[int] rMj(njpart); // M to recv from j
fespace Whi(Thi,Pk);
mesh Thij=Thi;
fespace Whij(Thij,Pk);// 

// construction of the mesh intersect i,j part 
for(int jp=0;jp<njpart;++jp)
  aThij[jp]  = trunc(Thi,pij[jp]>1e-6,label=10); // mesh of the supp of pij

for(int jp=0;jp<njpart;++jp)
  aThij[jp]  = trunc(aThij[jp],1,split=ksplit);

Thi =   trunc(Thi,1,split=ksplit);

settt 

if(ipart==0) cout << " *** end build mesh  intersection  " << endl;
// construction of transfert  matrix
{
  Whi wpii=pii;
  Pii = wpii[];
  for(int jp=0;jp<njpart;++jp)
    {
      int j=jpart[jp];
      Thij = aThij[jp];
      matrix I = interpolate(Whij,Whi); // Whji <- Whi
      sMj[jp] = I*Pii;  // Whi -> s Whij  
      rMj[jp] = interpolate(Whij,Whi,t=1);   // Whji -> Whi
      if(vdebug>10) {
      {Whi uuu=1;Whij vvv=-1; vvv[]+=I*uuu[]; cout << jp << " %%% " << vvv[].linfty << endl; assert(vvv[].linfty < 1e-6);}
      {Whi uuu=1;Whij vvv=-1; vvv[]+=rMj[jp]'*uuu[]; cout << jp << " ### " << vvv[].linfty << endl; assert(vvv[].linfty < 1e-6);}}
    }}
if(ipart==0) cout << " *** end build transfert matrix " << endl;
// alloc array of send and recv data ..

InitU(njpart,Whij,Thij,aThij,Usend)  // initU(n,Vh,Th,aTh,U)
InitU(njpart,Whij,Thij,aThij,Vrecv) // ...
if(ipart==0) cout << " *** end init  data for send/revc  " << endl;

Whi ui,vi; 

func bool Update(real[int] &ui, real[int] &vi)
{
  for(int j=0;j<njpart;++j)
    Usend[j][]=sMj[j]*ui; 
   SendRecvUV(comm,jpart,Usend,Vrecv)
     vi = Pii*ui;
   for(int j=0;j<njpart;++j)
     vi += rMj[j]*Vrecv[j][]; 
   return true;
}


// the definition of the Problem ....
func G=x*0.1; // ok 
func F=1.; // ok 
macro grad(u) [dx(u),dy(u)] // 
varf vBC(U,V)=  on(1,U=G); 
varf vPb(U,V)= int2d(Thi)(grad(U)'*grad(V)) + int2d(Thi)(F*V) + on(10,U=0)+on(1,U=G) ; //');// for emacs
varf vPbC(U,V)= int2d(ThC)(grad(U)'*grad(V))  +on(1,U=0) ; //');// for emacs
varf vPbon(U,V)=on(10,U=1)+on(1,U=1);
varf vPbon10only(U,V)=on(10,U=1)+on(1,U=0); 
// remark the order is important we want 0 part on 10 and 1


//----  

matrix Ai = vPb(Whi,Whi,solver=sparsesolver);
matrix AC,Rci,Pci;// 


if(mpiRank(comm)==0)
  AC = vPbC(VhC,VhC,solver=sparsesolver);

Pci=   interpolate(Whi,VhC);
Rci =  Pci'*Pii; 

real[int] onG10 = vPbon10only(0,Whi);
real[int] onG = vPbon(0,Whi);
real[int] Bi=vPb(0,Whi);





int kiter=-1;


func bool  CoarseSolve(real[int]& V,real[int]& U,mpiComm& comm)
{
   //  solving the coarse probleme 
   real[int] Uc(Rci.n),Bc(Uc.n); 
   Uc= Rci*U;
   mpiReduce(Uc,Bc,processor(0,comm),mpiSUM);
   if(mpiRank(comm)==0) 
      Uc = AC^-1*Bc;
    broadcast(processor(0,comm),Uc);
   V = Pci*Uc;
}//EOF ...
func real[int] DJ(real[int]& U)
{ 
  ++kiter;
  real[int] V(U.n); 
   V =  Ai*U;
  V = onG10 ? 0.: V;  // remove internal boundary  
  return V; 
}

func real[int] PDJ(real[int]& U) // C1
{ 
  real[int] V(U.n); 
  
  real[int] b= onG10 ? 0. :  U; 
  V =  Ai^-1*b;	
  Update(V,U);
  return U; 
}

func real[int] PDJC(real[int]& U) // 
{ // Precon  C1= Precon //, C2  precon Coarse
// Idea : F. Nataf. 
  //  0 ~  (I C1A)(I-C2A) => I ~  - C1AC2A +C1A +C2A 
  //  New Prec P= C1+C2 - C1AC2   = C1(I- A C2) +C2
  // (  C1(I- A C2) +C2 ) Uo 
  //   V =  - C2*Uo
  // .... 
  real[int] V(U.n); 
  CoarseSolve(V,U,comm);
  V = -V; //  -C2*Uo 
  U  += Ai*V; // U =  (I-A C2) Uo 
  real[int] b= onG10 ? 0. :  U; 
  U =  Ai^-1*b;	//  ( C1( I -A C2) Uo 
  V = U -V; //  
  Update(V,U);
  return U; 
}


 func real[int] DJ0(real[int]& U)
{ 
  ++kiter;
  real[int] V(U.n); 
  real[int] b= onG .* U;
  b  = onG ? b : Bi ;  
  V = Ai^-1*b;	
  Update(V,U);
  V -= U;
   return V; 
}


Whi u=0,v;
{ // verification.....
  Whi u=1,v;
  Update(u[],v[]);
  u[]-=v[];
  assert( u[].linfty<1e-6); }
  


settt
u[]=vBC(0,Whi,tgv=1); // set u with tge BC value ...

real epss=1e-6;
int rgmres=0;
if(gmres==1)
  {
   rgmres=MPIAffineGMRES(DJ0,u[],veps=epss,nbiter=300,comm=comm,dimKrylov=100,verbosity=ipart ? 0: 50);
   real[int] b= onG .* u[];
   b  = onG ? b : Bi ;
   v[] = Ai^-1*b;	
   Update(v[],u[]);
  }
else if(gmres==2)
  rgmres= MPILinearGMRES(DJ,precon=PDJ,u[],Bi,veps=epss,nbiter=300,comm=comm,dimKrylov=100,verbosity=ipart ? 0: 50);
else if(gmres==3)
   rgmres= MPILinearGMRES(DJ,precon=PDJC,u[],Bi,veps=epss,nbiter=300,comm=comm,dimKrylov=100,verbosity=ipart ? 0: 50);
else // algo Shwarz for demo ...
   for(int iter=0;iter <10; ++iter)
     {
       real[int] b= onG .* u[];
       b  = onG ? b : Bi ;
       v[] = Ai^-1*b;	
       
       Update(v[],u[]);
       if(vdebug) plotMPIall(Thi,u[],"u-"+iter);
        v[] -= u[];
       
       real err = v[].linfty;
       real umax = u[].max;
       real[int] aa=[err,umax], bb(2);
       mpiAllReduce(aa,bb,comm,mpiMAX);
       real errg = bb[0];
       real umaxg = bb[1];
       
       if(ipart==0)
	     cout << ipart << " err = " << errg << " u. max  " << umaxg << endl;
       if(errg< 1e-5) break;
     }
if(vdebug) plotMPIall(Thi,u[],"u-final");

settt

real errg =1,umaxg;
{ 
  real umax = u[].max,umaxg;
  real[int] aa=[umax], bb(1);
  mpiAllReduce(aa,bb,comm,mpiMAX);
  errg=bb[0];
  if(ipart==0)
    cout << " umax global  = " << bb[0] << " Wtime = " << (ttt[ittt-1]-ttt[ittt-2])  << " s " <<  " " << kiter <<  endl;
}

if(sff != "")
  {
    ofstream ff(sff+".txt",append);
    cout << " ++++  " ;
    cout  << mpirank <<"/" <<  mpisize << " k=" <<  ksplit << " n= " << nloc << " " << sizeoverlaps << " it=  " << kiter  ;  
    for (int i=1; i<ittt;++i)	
      cout << " " << ttt[i]-ttt[i-1] << " ";
    cout << epss << " " << Ai.nbcoef << " " << Ai.n << endl;
    
    /*
      1 mpirank
      2 mpisize 
      3 ksplit 
      4 nloc
      5 sizeoverlaps
      6 kiter
      7 mesh & part build  
      8 build the partion 
      9 build mesh, transfere , and the fine mesh .. 
      10 build the matrix,  the trans matrix, factorizatioon 
      11 GMRES 
    */
    ff   << mpirank << " " << mpisize << " " << sPk << " " ;
    ff <<  ksplit << " " << nloc << " " << sizeoverlaps << " " << kiter  ;  
    for (int i=1; i<ittt;++i)	
      ff << " " << ttt[i]-ttt[i-1] << " ";
    ff << epss << " " << Ai.nbcoef << " " << Ai.n << " " << gmres << endl;
    
  }
