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