Schur system + MatShell

[Barry Smith]

On Apr 21, 2008, at 5:54 AM, 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.

    Even if GMRES+Jacobi works reasonably well, GMRES without Jacobi  
can be much worse, this is a danger of matrix free without
some kind of preconditioner.

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

Run it with -ksp_monitor_true_residual (-ksp_truemonitor) for PETSc  
pre 2.3.3) and -ksp_converged_reason to see
what is happening. Note that KSPSolve() does NOT generate an error if  
it fails to converge, you need to check with
KSPGetConvergedReason() or -ksp_converged_reason after the solve to  
see if KSP thinks it has converged or
why it did not converge.

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