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

include "macro_ddm.idp"
load "PETSc"

mesh Th = square(40, 40);
func Pk = P1;
fespace Vh(Th, Pk);
varf vPb(u, v) = int2d(Th)(dx(u) * dx(v) + dy(u) * dy(v)) + on(1, 2, 3, 4, u = 0.0);
DmeshCreate(Th);
Mat A;
MatCreate(Th, A, Pk);
matrix Loc = vPb(Vh, Vh, tgv = -2);
A = Loc;
mesh ThTrunc;
matrix R;
real[int] bb(4);
boundingbox(Th, bb);
fespace VhTrunc(ThTrunc, Pk);
matrix prod;
if(bb[0] < 0.5) { // careful because some subdomains may not satisfy this condition!
    int[int] n2o;
    ThTrunc = trunc(Th, x < 0.5, new2old = n2o);
    int[int] rest = restrict(VhTrunc, Vh, n2o);
    real[int] V(VhTrunc.ndof);
    V = 1;
    int[int] I = 0:VhTrunc.ndof-1;
    R = [I, rest, V];
    R.resize(VhTrunc.ndof, Vh.ndof);
    prod = Loc * R';
}
Mat B(A, restriction = R); // if the above condition is not met,
Mat C(A, B, prod);         // these two distributed matrices won't have any local unknowns
real[int] values;
Vh[int] lvec(10);
VhTrunc[int] rvec(10);
int nsv = SVDSolve(C, sparams = "-svd_largest -svd_view_values -svd_type cyclic -svd_nsv 10", values = values, lvectors = lvec, rvectors = rvec);
func real[int] prodFunc(real[int]& up) {
    real[int] u;
    ChangeNumbering(B, u, up, inverse = true, exchange = true); // from PETSc to FreeFEM numbering + need to exchange ghost values
    real[int] v = R' * u;
    u.resize(Loc.n);
    u = A * v;
    ChangeNumbering(A, u, up); // from FreeFEM to PETSc numbering
    return up;
}
func real[int] prodFuncTranspose(real[int]& utp) {
    real[int] ut;
    ChangeNumbering(A, ut, utp, inverse = true, exchange = true); // from PETSc to FreeFEM numbering + need to exchange ghost values
    real[int] v = A' * ut;
    ut.resize(R.n);
    ut = R * v;
    ChangeNumbering(B, ut, utp); // from FreeFEM to PETSc numbering
    return utp;
}
Mat MF(A, B, prodFunc, transpose = prodFuncTranspose); // defines a "rectangular" operator (I = A.size and J = B.size)
ObjectView(MF, format = "info");
real[int] valuesMF;
real[int, int] lvecMF(A.n, 10);
real[int, int] rvecMF(B.n, 10);
int nsvMF = SVDSolve(MF, sparams = "-svd_largest -svd_type cyclic -svd_nsv 10", values = valuesMF, larray = lvecMF, rarray = rvecMF);
values -= valuesMF;
assert(values.linfty < 1.0e-4 && nsv == nsvMF);

Vh u, uMF, du;
VhTrunc v, vMF, dv;
for (int idSV = 0; idSV < nsv; ++idSV)
{
  // check same (or opposite) left vectors
  ChangeNumbering(A, uMF[], lvecMF(:,idSV), inverse = true, exchange = true);
  du[] = lvec[idSV][] - uMF[];
  if(du[].linfty >= 1.0e-6) { // sign flip?
      du[] = lvec[idSV][];
      du[] *= -1.0;
      du[] -= uMF[];
  }
  assert(du[].linfty < 1.0e-4);

  // check same (or opposite) right vectors
  ChangeNumbering(B, vMF[], rvecMF(:,idSV), inverse = true, exchange = true);
  dv[] = rvec[idSV][] - vMF[];
  if(dv[].linfty >= 1.0e-6) {
      dv[] = rvec[idSV][];
      dv[] *= -1.0;
      dv[] -= vMF[];
  }
  assert(dv[].linfty < 1.0e-4);

  // check C v = sigma u
  real[int] temp1(A.n), temp2(lvecMF.n);
  MatMult(MF, rvecMF(:,idSV), temp1); // rectangular MF, so use vectors un PETSc numbering + MatMult
  temp2 = valuesMF[idSV] * lvecMF(:,idSV);
  temp2 -= temp1;
  assert(temp2.linfty < 1.0e-4);
  MatMult(MF, rvecMF(:,idSV), temp2);
  ChangeNumbering(A, temp1, temp2, inverse = true, exchange = true);
  temp2.resize(uMF[].n);
  temp2 = MF * vMF[];
  temp2 -= temp1;
  assert(temp2.linfty < 1.e-4);
}
