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

assert(mpisize == 1);
load "PETSc"
{
    real[int, int] A(2, 2);
    A(0, 0) = 2;
    A(0, 1) = 4;
    A(1, 1) = -4;
    A(1, 0) = 0.3;
    matrix sp = A;
    Mat dA(sp);
    real[int] list(2);
    list = 0;
    list(1) = 1;
    real[int, int] S(1, 1);
    set(dA, sparams = "-ksp_type preonly -pc_type lu -pc_factor_mat_solver_type mumps", schurList = list, schurComplement = S);
    cout << S << endl;
    assert(abs(S(0, 0) - (A(1, 1) - 1.0 / A(0, 0) * A(1, 0) * A(0, 1))) < 1e-10);
}
include "ExtractDofsonBorder.idp"
macro grad(u) [dx(u), dy(u)]//
real h = 0.1;
mesh Th = square(2/h, 1/h, [2*x, y]);
Th = trunc(Th, (abs(x-1) > 3*h) || (abs(y-0.45) > h/2), label = 10);
plot(Th, wait = 1);
fespace Vh(Th,P1);
varf vA(u,v) = int2d(Th)(grad(u)' * grad(v)) + on(1, 3, 4, u = 0);
{
    matrix A = vA(Vh, Vh, tgv = -2);
    Mat dA(A);
    int[int] I;
    ExtractDofsonBorder(2, Vh, I, 1);
    real[int] list(Vh.ndof);
    for [j, ij : I] list[ij] = j + 1;
    real[int, int] S(1, 1);
    set(dA, sparams = "-ksp_type preonly -pc_type lu -pc_factor_mat_solver_type mumps", schurList = list, schurComplement = S);
    cout << S << endl;
}
