Commits

Mark Adams committed 86f5270

shift coarse grids to avoid zero diags from low rank coarse grid space

Comments (4)

    1. Mark Adams author

      I tried things like 1/PETSC_MAX_REAL (FP_TRAPs in MatShift) and then 10000000/PETSC_MAX_REAL (FP_TRAPs in GMRES). So this is what I got to work. Open to suggestions.

Files changed (1)

src/ksp/pc/impls/gamg/gamg.c

   MPI_Comm        comm;
   PetscMPIInt     rank,size,new_size,nactive=*a_nactive_proc;
   PetscInt        ncrs_eq,ncrs,f_bs;
+  const PetscScalar shift = PETSC_MACHINE_EPSILON*PETSC_MACHINE_EPSILON*PETSC_MACHINE_EPSILON;
 
   PetscFunctionBegin;
   ierr = PetscObjectGetComm((PetscObject)Amat_fine,&comm);CHKERRQ(ierr);
   ierr = MatGetBlockSize(Amat_fine, &f_bs);CHKERRQ(ierr);
   /* RAP */
   ierr = MatPtAP(Amat_fine, Pold, MAT_INITIAL_MATRIX, 2.0, &Cmat);CHKERRQ(ierr);
+  ierr = MatShift(Cmat,shift);CHKERRQ(ierr);
 
   /* set 'ncrs' (nodes), 'ncrs_eq' (equations)*/
   ierr = MatGetLocalSize(Cmat, &ncrs_eq, NULL);CHKERRQ(ierr);
       ierr = PCReset_MG(pc);CHKERRQ(ierr);
       pc->setupcalled = 0;
     } else {
+      const PetscScalar shift = PETSC_MACHINE_EPSILON*PETSC_MACHINE_EPSILON*PETSC_MACHINE_EPSILON;
       PC_MG_Levels **mglevels = mg->levels;
       /* just do Galerkin grids */
       Mat          B,dA,dB;
           if (pc_gamg->setup_count==2) {
             ierr = MatPtAP(dB,mglevels[level+1]->interpolate,MAT_INITIAL_MATRIX,2.0,&B);CHKERRQ(ierr);
             ierr = MatDestroy(&mglevels[level]->A);CHKERRQ(ierr);
-
             mglevels[level]->A = B;
           } else {
             ierr = KSPGetOperators(mglevels[level]->smoothd,NULL,&B);CHKERRQ(ierr);
             ierr = MatPtAP(dB,mglevels[level+1]->interpolate,MAT_REUSE_MATRIX,2.0,&B);CHKERRQ(ierr);
           }
+          ierr = MatShift(B,shift);CHKERRQ(ierr);
           ierr = KSPSetOperators(mglevels[level]->smoothd,B,B);CHKERRQ(ierr);
           dB   = B;
         }