[petsc-users] Preconditioning with SLEPc

[Jose E. Roman]

Currently there is no straightforward way to do this in SLEPc. For a preconditioned eigensolver such as JD you could use this:
http://slepc.upv.es/documentation/current/docs/manualpages/ST/STPrecondSetMatForPC.html
in combination with pc=PCMAT.

For a Krylov solver you will have to do something like this:

      call EPSSetWhichEigenpairs(eps,EPS_TARGET_REAL,ierr);CHKERRA(ierr)
      call EPSSetTarget(eps,sigma,ierr);CHKERRA(ierr)
      call EPSGetST(eps,st,ierr);CHKERRA(ierr)
      call STSetType(st,STSINVERT,ierr);CHKERRA(ierr)
      call STGetKSP(st,ksp,ierr);CHKERRA(ierr)
      call KSPGetPC(ksp,pc,ierr);CHKERRA(ierr)
      call KSPSetType(ksp,KSPBCGS,ierr);CHKERRA(ierr)
      call PCSetType(pc,PCNONE,ierr);CHKERRA(ierr)  ! dummy preconditioner

      ! provide the preconditioner matrix after setup
      call EPSSetUp(eps,ierr);CHKERRA(ierr)
      call PCSetOperators(pc,PETSC_NULL_MAT,my_matmul_K,ierr);CHKERRA(ierr)
      call PCSetType(pc,PCMAT,ierr);CHKERRA(ierr)

We may change the interface in future versions so that this can be done in a more natural way.

Jose


> El 7 ago 2018, a las 16:42, Kenneth Hall, SC.D. <kenneth.c.hall at duke.edu> escribió:
> 
> Hello,
> 
> I am new to SLEPc and Petsc, so I would appreciate any help I can get with my question. I have a large eigenvalue problem I am solving, based on the linearization of a nonlinear code.  The eigenvalue problem takes the form
> 
> A.x = lambda.B.x
> 
> I can generate the matrix B explicitly, but I can only form products of A.x, and not A itself, so I use a shell matrix.   I also have an iterative solver I could potentially use as a preconditioner. That is, I have an operator K^-1 where
> 
> K^-1 \approx  (A - sigma.B)^-1
> 
> But again, I don’t have K^-1 explicitly, rather I can form products of K^-1 with a vector, so I would need to use a shell matrix.
> 
> I have had success with the un-preconditioned case, but can’t figure out how to do the preconditioned case.  Any suggestions would be greatly appreciated.  I work with Fortran.  A code snippet with the un-preconditioned set-up follows. 
> 
> Thanks in advance.
> Kenneth Hall
> 
> 
> ===============================================================
> 
> !
> !.... initialize slepc
>       call SlepcInitialize(PETSC_NULL_CHARACTER,ierr) ; if (ierr .ne. 0) then ; stop 'SlepcInitialize failed' ; end if
> !
> !.... get size and rank of MPI process
>       call MPI_Comm_size(PETSC_COMM_WORLD,sz  ,ierr)                                               ; CHKERRA(ierr)
>       call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)                                               ; CHKERRA(ierr)
> !
>       if (sz .ne. 1) then ; SETERRA(PETSC_COMM_SELF,1,'This is a uniprocessor example only!'); endif
> !
> !.... get user options
>       call PetscOptionsGetInt (PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-n'   ,n ,flg,ierr)        ; CHKERRA(ierr)
>       call PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-terse',terse,ierr)        ; CHKERRA(ierr)
> !
> !.... write out option list
>       if (rank .eq. 0 ) then
>          write(*,*) 'n        = ',n
>       end if
> !
> !.... Register the matrix-vector subroutine for the operator that defines the eigensystem, Ax=kBx
>       call MatCreateShell(PETSC_COMM_WORLD,n,n,n,n,PETSC_NULL_INTEGER,A,ierr)                      ; CHKERRA(ierr)
>       call MatShellSetOperation(A,MATOP_MULT,my_matmul_A,ierr)                                     ; CHKERRA(ierr)
> !
> !.... set up B matrix
>       call MatCreate(PETSC_COMM_WORLD,B,ierr)                                                      ; CHKERRA(ierr)
>       call MatSetSizes(B,n,n,n,n,ierr)                                                             ; CHKERRA(ierr)
>       call MatSetFromOptions(B,ierr)                                                               ; CHKERRA(ierr)
>       call MatSetUp(B,ierr)                                                                        ; CHKERRA(ierr)
>       call My_Matrix_B(B,n)
> !
> !.... set up eigenvectors
>       call VecCreate(PETSC_COMM_WORLD,xr,ierr)                                                      ; CHKERRA(ierr)
>       call VecSetSizes(xr,n,n,ierr)                                                                 ; CHKERRA(ierr)
>       call VecSetFromOptions(xr,ierr)                                                               ; CHKERRA(ierr)
>       call VecSetUp(xr,ierr)                                                                        ; CHKERRA(ierr)
> !
>       call VecCreate(PETSC_COMM_WORLD,xi,ierr)                                                      ; CHKERRA(ierr)
>       call VecSetSizes(xi,n,n,ierr)                                                                 ; CHKERRA(ierr)
>       call VecSetFromOptions(xi,ierr)                                                               ; CHKERRA(ierr)
>       call VecSetUp(xi,ierr)                                                                        ; CHKERRA(ierr)
> !
> !.... Create eigensolver context
>       call EPSCreate(PETSC_COMM_WORLD,eps,ierr)                                                    ; CHKERRA(ierr)
> !
> !.... Set operators. In this case, it is a generalized eigenvalue problem
>       call EPSSetOperators(eps,A,B,ierr)                                                           ; CHKERRA(ierr)
>       call EPSSetProblemType(eps,EPS_GNHEP,ierr)                                                   ; CHKERRA(ierr)
> !
> !.... set default options
>       call EPSSetType(eps,EPSKRYLOVSCHUR,ierr)                                                     ; CHKERRA(ierr)
> !
> !.... set default eigenvalues to find
>       call EPSSetWhichEigenpairs(eps,EPS_LARGEST_REAL,ierr)                                        ; CHKERRA(ierr)
> !
> !.... number of requested eigenvalues
>       call EPSSetDimensions(eps,nev,PETSC_DEFAULT_INTEGER,PETSC_DEFAULT_INTEGER,ierr)              ; CHKERRA(ierr)
> !
> !.... set stopping criteria
>       call EPSSetTolerances(eps,tol,maxit,ierr)                                                    ; CHKERRA(ierr)
> !
> !.... set up monitors
>       call PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_DEFAULT,vf,ierr)      ; CHKERRA(ierr)
>       call EPSMonitorSet(eps,EPSMONITORFIRST,vf,PetscViewerAndFormatDestroy,ierr)                  ; CHKERRA(ierr)
>       call SlepcConvMonitorCreate(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_DEFAULT,ctx,ierr)         ; CHKERRA(ierr)
>       call EPSMonitorSet(eps,EPSMONITORCONVERGED,ctx,SlepcConvMonitorDestroy,ierr)                 ; CHKERRA(ierr)
> !
> !.... Set solver parameters at runtime
>       call EPSSetFromOptions(eps,ierr)                                                             ; CHKERRA(ierr)
> !
> !.... Solve the eigensystem
>       call EPSSolve(eps,ierr)                                                                      ; CHKERRA(ierr)
>