//  run with MPI:  ff-mpirun -np 4 script.edp
// NBPROC 4

load "MUMPS_mpi" // add _mpi for windows do day
// load "symmetrizeCSR"
int[int] l = [1, 1, 2, 2];
mesh Th = square(150, 150, label = l);
int master = 0; //  no choose see MUMPS doc sect 4.3 
if(mpirank != master) // no need to store the matrix on ranks other than 0
    Th = square(1, 1, label = l);
fespace Vh(Th, P2);
real[int] rinfo(40);
 int[int] info(40); 

varf lap(u,v) = int2d(Th)(dx(u)*dx(v) + dy(u)*dy(v) - 100*u*v) + int2d(Th)(v) + on(1, u = 1);
real[int] b = lap(0, Vh);
matrix A,B;
//if( mpirank==master) // => trap pb of size  in  A^-1 * b;
{ // def only on master .. 
	A = lap(Vh, Vh); B = A;
	//symmetrizeCSR(B);
	set(B,sym=1);
	
	cout<< " A nnz "<< A.nnz << endl;
	cout<< " B nnz "<< B.nnz << endl;
} 
verbosity = 2;
// If you want, you can pass an array of integer (int[int] t)
// for specifying the ordering you want to use
// 1 : AMD
// 3 : AMF
// 4 : SCOTCH
// 5 : PORD
// 6 : METIS
// 7 : QAMD
// 10 : PT-SCOTCH
// 11 : ParMETIS
//
// You can also set on which MPI rank the matrix is defined
// as the current interface only supports centralized matrices
//
// Example :
// If one wants to use SCOTCH with a matrix stored on mpirank == 0
// int[int] t(2); t(0) = 4; t(1) = 1;
// set(A, solver = sparsesolver, lparams = t);
set(A, solver = sparsesolver,master=master);
set(B, solver = sparsesolver,master=master,rinfo=rinfo,info=info);
verbosity = 0;
Vh u;
u[] = A^-1 * b;
if(mpirank==master) cout << " A :  u " << u[].l1 << " " << u[].linfty << endl; 
if(mpirank==0-master)
plot(Th, u, wait = 1, dim = 3, fill = 1, cmm = "sparsesolver unsymmetric", value = 1);
u[] = B^-1 * b;
if(mpirank==master) cout << " B :  u " << u[].l1 << " " << u[].linfty << endl; 
if(mpirank==master)
for(int i=0; i<40;++i)
   {
    cout << "mumps: info B" << i+1 << " = " << info[i] << " , / r " << rinfo[i] <<endl;
   }

if(mpirank==master)
plot(Th, u, wait = 1, dim = 3, fill = 1, cmm = "sparsesolver symmetric", value = 1);

//  try distributed  version ...
 Th = square(150, 150, label = l);
u=0; // resize u..
b.resize(Vh.ndof);
Th=change(Th,fregion= nuTriangle%mpisize);
varf lapd(u,v) = int2d(Th,mpirank)(dx(u)*dx(v) + dy(u)*dy(v)) + int2d(Th,mpirank)(v) + on(1, u = 1);
A = lapd(Vh, Vh,solver=CG);
set(A, solver = sparsesolver,master=-1);// distruteb
b = lapd(0, Vh);
u[] = A^-1 * b;
cout << " D u " << u[].l1 << " " << u[].linfty << " " << mpirank << endl; 
cout << " On Fini" << endl; 
