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

load "PETSc"                        // PETSc plugin
macro trueRestrict()true//
macro removeZeros()true//
include "macro_ddm.idp"             // additional DDM functions

macro def(i)i// EOM                 // scalar field definition
macro init(i)i// EOM                // scalar field initialization
macro grad(u)[dx(u), dy(u)]// EOM   // two-dimensional gradient
func Pk = P1;                       // finite element space


mesh Th = square(getARGV("-global", 40), getARGV("-global", 40)); // global mesh
mesh ThBackup = Th;
int[int] n2o;
macro ThN2O()n2o//
DmeshCreate(Th);
fespace Wh(Th, Pk);           // local finite element space

varf vPb(u, v) = int2d(Th)(u * v + grad(u)' * grad(v) - 10 * dy(v) * u) + int2d(Th)(v * x) + on(1, 2, 3, 4, u = 0.0);
varf vPbT(u, v) = int2d(Th)(u * v + grad(u)' * grad(v) - 10 * dy(u) * v) + int2d(Th)(v * x) + on(1, 2, 3, 4, u = 0.0);
matrix<real> Loc  = vPb(Wh, Wh);
matrix<real> LocT = vPbT(Wh, Wh);
real[int] rhs = vPb(0, Wh);

Mat A;
MatCreate(Th, A, Pk);
A = Loc;
Mat AT(A, LocT);
Wh<real> def(u);  // local solution
Wh<real> def(v);  // local solution

u[] = A^-1 * rhs;
plotMPI(Th, u, Pk, def, real, cmm = "Global solution");
u[] = A'^-1 * rhs;
plotMPI(Th, u, Pk, def, real, cmm = "Global solution transpose");
v[] = AT^-1 * rhs;
u[] -= v[];
real local = u[].linfty;
real norm;
mpiAllReduce(local, norm, mpiCommWorld, mpiMAX);
plotMPI(Th, v, Pk, def, real, cmm = "Global solution transposed problem (linfty of error = " + norm + ")");
if(!NoGraphicWindow) {
    fespace WhBackup(ThBackup, Pk);
    int[int] rest = restrict(Wh, WhBackup, n2o);
    u[] = v[];
    u[] .*= A.D;
    WhBackup sol, solReduced;
    for[i, v : rest] sol[][v] = u[][i];
    mpiReduce(sol[], solReduced[], processor(0), mpiSUM);
    plot(solReduced, wait = 1);
}
