Schur system + MatShell

[Satish Balay]

On Mon, 21 Apr 2008, tribur at vision.ee.ethz.ch wrote:

> 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.

Hmm - with dense - if you have some other efficient way of assembling
the matrix - you can specify this directly to MatCreateMPIDense() - [or
use MatGetArray() - and set the values directly into this array]

> 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

Which version of PETSc is this? I can't place the line numbers
correctly with latest petsc-2.3.3. [Can you send the complete error
trace?]

> 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.

Not sure what the problem could be. Can you confirm that the code is
valgrind clean? It could explain the issue 1 aswell.

With mpich2 you can do the following on linux:

 mpiexec -np 2 valgrind --tool=memcheck ./executable

Satish

> 
> 
> 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);
> ...
> 
> 
>