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

load "PETSc"

int[int] I(1), J(1);
I(0) = 0;
J(0) = 0;

if(mpisize != 4)
    exit(1);

int[int] Ji;
matrix[int] Ai(0);
if(mpirank == 0) {
    Ji.resize(2);
    Ai.resize(2);
    Ji[0] = 0; Ji[1] = 2;
    real[int] C(1);
    C(0) = 10.0;
    Ai[0] = [I, J, C];
    C(0) = 2.0;
    Ai[1] = [I, J, C];
}
else if(mpirank == 1) {
    Ji.resize(2);
    Ai.resize(2);
    Ji[0] = 1; Ji[1] = 2;
    real[int] C(1);
    C(0) = -10.0;
    Ai[0] = [I, J, C];
    C(0) = 3.0;
    Ai[1] = [I, J, C];
}
else if(mpirank == 2) {
    Ji.resize(2);
    Ai.resize(2);
    Ji[0] = 0; Ji[1] = 2;
    real[int] C(1);
    C(0) = -1.0;
    Ai[0] = [I, J, C];
    C(0) = 5.0;
    Ai[1] = [I, J, C];
}
else if(mpirank == 3) {
    Ji.resize(1);
    Ai.resize(1);
    Ji[0] = 3;
    real[int] C(1);
    C(0) = 7.0;
    Ai[0] = [I, J, C];
}
Mat A(Ai, columns = Ji);
real[int] x(1), b(1);
b = (mpirank + 1) * 10;
x = A^-1 * b;
for(int i = 0; i < mpisize; ++i) { mpiBarrier(mpiCommWorld); if(i == mpirank) cout << "rhs on process #" << mpirank << ": " << b[0] << ", solution: " << x[0] << endl; mpiBarrier(mpiCommWorld); }
real[int] y(1);
y = A * x;
for(int i = 0; i < mpisize; ++i) { mpiBarrier(mpiCommWorld); if(i == mpirank) cout << "A * x on process #" << mpirank << ": " << y[0] << ", rhs: " << b[0] << endl; mpiBarrier(mpiCommWorld); }
