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

load "MUMPS_mpi"  // load the library dynamic correspond to MUMPS interface
verbosity = 0;
int[int] ICNTL(40); // declaration of ICNTL parameter for MUMPS

// get value of ICNTL from file
if(mpirank == 0)
{
	ifstream ff("ffmumps_fileparam.txt");
	string line;
	getline(ff,line);
	getline(ff,line);
	for(int iii=0; iii<40;iii++){
	ff >> ICNTL[iii];
	getline(ff,line);
	}
}

broadcast(processor(0),ICNTL);

// Given data of MUMPS solver in array lparams(SYM,PAR,ICNTL)
// There is no symmetric storage for a matrix associated with a sparse solver.
// Therefore, the matrix will be considered unsymmetric for parallel sparse solver even if symmetric.
// This fact will be change in new version of FreeFEM
{
    int SYM = 0;
    int PAR = 1;
    matrix A =
      [[ 40,  0,     45,  0, 0],
       [ 0,    12,     0,  0 , 0],
       [ 0,     0,  40,  0, 0],
       [ 12,     0,    0,  22, 0],
       [ 0,     0,    20,  0., 22]];

    // construction of integer parameter for MUMPS
    int[int] MumpsLParams(42);
    MumpsLParams[0] = SYM;
    MumpsLParams[1] = PAR;
    for(int ii=0; ii< 40; ii++) MumpsLParams[ii+2] = ICNTL[ii]; // ICNTL begin with index 0 here

    real[int] xx = [ 1,32,45,7,2], x(5), b(5), di(5);
    b=A*xx;
    if(mpirank ==0) cout << "xx=" << xx << endl;

    set(A,solver=sparsesolver,lparams=MumpsLParams); // We take the default value for CNTL MUMPS parameter

    if(mpirank ==0)  cout << "solving solution" << endl;
    x = A^-1*b;
    if(mpirank ==0) cout << "b=" << b << endl;
    if(mpirank ==0) cout << "x=" << endl; cout << x << endl;
    di = xx-x;
    if(mpirank==0){
    cout << "x-xx="<< endl; cout << "Linf "<< di.linfty << " L2 " << di.l2 << endl;
    }
}

// FFCS - reference value for regression tests
real regtest;

// Read parameter of MUMPS solver in file ffmumps_fileparam.txt

{
matrix A =
      [[ 40,  0,     45,  0, 0],
       [ 0,    12,     0,  0 , 0],
       [ 0,     0,  40,  0, 0],
       [ 12,     0,    0,  22, 0],
       [ 0,     0,    20,  0., 22]];


    real[int] xx = [ 1,32,45,7000,2], x(5), b(5), di(5); // xb(4),bbb(4);
    b=A*xx;
    cout << "b="  << b  << endl;
    cout << "xx=" << xx << endl;

    set(A,solver=sparsesolver);

    cout << "solving solution" << endl;
    x = A^-1*b;

    cout << "b=" << b << endl;
    cout << "x=" << endl; cout << x << endl;
    di = xx-x;
    if(mpirank==0){
    cout << "x-xx="<< endl; cout << "Linf "<< di.linfty << " L2 " << di.l2 << endl;
    regtest=di.l2;
    }
}
