Schur system + MatShell

[tribur at vision.ee.ethz.ch]

Dear all,

Sorry for switching from Schur to Hypre and back, but I'm trying two 
approaches at the same time to find the optimal solution for our 
convection-diffusion/Stokes problems: a) solving the global stiffness 
matrix directly and in parallel using Petsc and a suitable 
preconditioner (???) and b) applying first non-overlapping domain 
decomposition and than solving the Schur complement system.

Being concerned with b in the moment, I managed to set up and solve the 
global Schur system using MATDENSE. The solving works well with, e.g., 
gmres+jacobi, but the assembling of the global Schur matrix takes too 
long. Therefore, I'm trying to use the matrix in unassembled form using 
MatShell. Not very successfully, however:

1) When I use KSPGMRES, I got the error
[1]PETSC ERROR: MatMult() line 1632 in src/mat/interface/matrix.c
[1]PETSC ERROR: PCApplyBAorAB() line 584 in src/ksp/pc/interface/precon.c
[1]PETSC ERROR: GMREScycle() line 159 in src/ksp/ksp/impls/gmres/gmres.c
[1]PETSC ERROR: KSPSolve_GMRES() line 241 in src/ksp/ksp/impls/gmres/gmres.c
[1]PETSC ERROR: KSPSolve() line 379 in src/ksp/ksp/interface/itfunc.c

2) Using KSPBICG, it iterates without error message, but the result is 
wrong (norm of residual 1.42768 instead of something like 1.0e-10), 
although my Mat-functions PETSC_SchurMatMult and 
PETSC_SchurMatMultTranspose seem to be correct. I tested the latter 
comparing the vectors y1 and y2 computed by, e.g., MatMult(S,x,y1) and 
PETS_SchurMatMult(S,x,y2). Norm(y1-y2) was < e-15 for both functions.


Could you please have a look at my code snippet below?

Thank you very much!
Kathrin



PS: My Code:

Vec gtot, x;
...
Mat Stot;    IS is;
ISCreateGeneral(PETSC_COMM_SELF, NPb, &uBId_global[0], &is);
localData ctx;
ctx.NPb  = NPb;   //size of local Schur system S
ctx.Sloc = &S[0];
ctx.is = is;
MatCreateShell(PETSC_COMM_WORLD,m,n,NPb_tot,NPb_tot,&ctx,&Stot);
MatShellSetOperation(Stot,MATOP_MULT,(void(*)(void)) 
PETSC_SchurMatMult); 
MatShellSetOperation(Stot,MATOP_MULT_TRANSPOSE,(void(*)(void))PETSC_SchurMatMultTranspose);
KSP ksp;
KSPCreate(PETSC_COMM_WORLD,&ksp);
PC prec;
KSPSetOperators(ksp,Stot,Stot,DIFFERENT_NONZERO_PATTERN);
KSPGetPC(ksp,&prec);
PCSetType(prec, PCNONE);
KSPSetType(ksp, KSPBICG);
KSPSetTolerances(ksp, 1.e-10, 1.e-50,PETSC_DEFAULT,PETSC_DEFAULT);
KSPSolve(ksp,gtot,x);
...