// NBPROC 4
load "MPICG"

// the grad of the bilinear part of J (the RHS in remove)
   int n=3;
  real[int]  b(n),u(n);
  b=1;
   u=0; // set  right hand side and initial gest

  func real[int] matId(real[int] &u) { return u;};
  func real[int] DJ0(real[int] &u)
    { 
      for (int i=0;i<u.n;i++)
	u[i]=10*mpirank+(i+1)*u[i];
      return u;  // return of global variable ok 
   };
  func real error(real[int] & u,real[int] & b)
   {
   real s=0;
     for (int i=0;i<u.n;i++)
	s += abs(10*mpirank+(i+1)*u[i] - b[i]);
   return s;    
   }

  b=1; u=0; // set  right hand side and initial gest
  MPILinearCG(DJ0,u,b,eps=1.e-6,nbiter=20,precon=matId);
  cout << "LinearCG (Linear)  err=" << error(u,b) << endl;
  assert(error(u,b) < 1e-5);


{
  int[int] procs=[0,1];
  mpiGroup gpr(procs);
  cout << " before  ...  " << mpirank << endl;
  mpiComm comm(gpr);
  if(comm)  // ..  comm 
    {
      int rank=mpiRank(comm);
      cout << " rank = " << rank << endl; 
      mesh[int]  aTh(2);
      int[int]  ll=[1,2-rank,1,1+rank];
      aTh[rank]=square(4,4,[rank+x,y], label=ll);	  
      broadcast(processor(0,comm),aTh[0]);
      broadcast(processor(1,comm),aTh[1]);
      mesh Th=aTh[rank];
      mesh Thi=aTh[1-rank]; 
      mesh Thh = Th+Thi;
      
      
      fespace Vh(Th,P1);	
      fespace Vhi(Thi,P1);
      
      matrix Ihi=interpolate(Vh,Vhi, inside=1);// Vh <-Vhi
      cout << rank << " Ihi = "<< Ihi << " ---- " << endl;
      
      // compresse row numbering to min. th comm.
      int [int] Ii(Vhi.ndof),II(Vh.ndof);
      real[int] dhi(Vh.ndof);
      int nn,NN;
      { 
	Vh uh=1.; Vhi ui=1;
	real[int] vi(0:ui[].n-1);
	ui[] =Ihi'*uh[];//';
	ui[] = ui[] ? 1 : 0;
	nn = rint(ui[].sum);
	vi = ui[] ? vi : 1e30;
	vi.sort;
	//	   cout << vi[nn] << " " << vi[nn-1] << endl;
	assert(vi[nn]>1e10 && vi[nn-1] < 1e10);
	vi.resize(nn);
	Ii.resize(nn);
	for(int i=0;i<nn;++i)
	  Ii[i]=vi[i];
	processor(1-rank,comm) << nn;
	processor(1-rank,comm) >> NN;
	assert( nn == NN) ;
	//   II.resize(NN);
	processor(1-rank,comm) <<Ii;
	processor(1-rank,comm) >>II;	     
	//cout << Ii << endl;
	real Ihioo=0;
	try {
	  Ihioo= Ihi(Vh.ndof-1,Ii[nn-1]);
	}
	   catch(...) {;};
	int[int] I(0:Vh.ndof-1);
	Ihi = Ihi(I,Ii);
	//assert(Ihi(Vh.ndof-1,nn-1) == Ihioo);
	// cout << rank << " Ihi = "<< Ihi << " ---- " << endl;
	dhi=1;
	vi=1;
	dhi += Ihi*vi;
	real[int] one(dhi);one=1.;
	dhi = one ./ dhi;
	cout << dhi << endl;
      }
      
      real[int] tosend(NN),torecv(nn);
      real tgv=1e3;
      macro  Grad(u) [dx(u),dy(u)] //
	varf va(u,v)=int2d(Th)(u*v+Grad(u)'*Grad(v)) - int2d(Th)( x*v ) -int1d(Th,1)( N.x*v)  ;	      ;//');
      matrix A=va(Vh,Vh,tgv=tgv);
      real[int] b=va(0,Vh,tgv=tgv);   
      int kiter=0;
      func real[int] projC(real[int] &u)
      {
	verbosity=0;
	processor(1-rank,comm) << (tosend = u(II));
	processor(1-rank,comm) >> torecv;
	//cout << "to send " << tosend << endl;
	u += Ihi*torecv;
	//  u = u .* dhi;
	Vh uu; uu[]=u;
	//     cout << kiter++ << " " << mpirank <<  " ---  " <<  uu(1,0.5) << endl; 	      
	return   u;	      
      }	  
      
      func real[int] projD(real[int] &u)
      {
	verbosity=0;
	processor(1-rank,comm) << (tosend = u(II));
	processor(1-rank,comm) >> torecv;
	u += Ihi*torecv;
	u = u .* dhi;
	Vh uu; uu[]=u;
	//   cout << kiter++ << " " << mpirank <<  " ---  " <<  uu(1,0.5) << endl; 	      
	return   u;	      
      }	
  
      func bool  Plot(real[int] & uu)
      {
	Vh u; u[]=uu;
	Vhi ui;
	if(rank==0)
	  {     
	    processor(1,comm) >> ui[];
	    plot(u,ui,wait=1);
	  }
	else
	  processor(0,comm) << u[];
	return true;
      }
      
      func real[int] DJJ(real[int] &v)
      { 
	verbosity=0;
	// Plot(v);
	v = projD(v);
	real[int] u=A*v;
	u += b; 
	u = projC(u);
	verbosity=1;
	return u;  // return 
      };		  
      
      func real[int] matId(real[int] &u) { 
	verbosity=0;
	u =projD(u);
	verbosity=1;
	return u;};
      
      
      
      
      Vh u=0;
      
      u=0;
      verbosity=1;
      MPIAffineCG(DJJ,u[],eps=1.e-4,nbiter=200,precon=matId,comm=comm);
      projD(u[]);
      Plot(u[]);
      if(rank==0)
	{
	  fespace Vhh(Thh,P1);
	  Vhh uuu,vvv;
	  // u = x , dn(x) = N.x 
	  solve PPPP(uuu,vvv) = int2d(Thh)(uuu*vvv+Grad(uuu)'*Grad(vvv)) - int2d(Thh)( x*vvv ) -int1d(Thh,1)( N.x*vvv)  ; //');
	  cout << " err = " << int2d(Th)( square(uuu-u)) << " " << int2d(Thh)( square(uuu-x)) << endl;
	 // savemesh(Thh,"/tmp/Thh.msh");
	}
    }
}

