[petsc-users] Preconditioning with SLEPc
Matthew Knepley
knepley at gmail.com
Wed Aug 8 11:02:18 CDT 2018
On Wed, Aug 8, 2018 at 11:51 AM Kenneth Hall, SC.D. <kenneth.c.hall at duke.edu>
wrote:
> Roman,
>
> When I tried what you suggested, I get the following error:
>
> call
> PCSetOperators(pc,PETSC_NULL_MAT,my_matmul_k,ierr);CHKERRA(ierr)
> 1
> Error: Invalid procedure argument at (1)
> make: *** [/Users/hall/git/Mustang/Mustang/build_dir/slepc_eigen2.o] Error
> 1
>
>
> I assume this is because PCSetOperators is expecting a MAT type rather
> than a subroutine. I modified your code snippet using a shell matrix, with
> the result:
>
> 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)
> call EPSSetUp(eps,ierr);CHKERRA(ierr)
>
> ! begin modified code
>
>
> call MatCreateShell(PETSC_COMM_WORLD,n,n,n,n,PETSC_NULL_INTEGER,P,ierr) ; CHKERRA(ierr)
> call MatShellSetOperation(P,MATOP_MULT,my_matmul_k,ierr); CHKERRA(ierr)
> call PCSetOperators(pc,PETSC_NULL_MAT,P,ierr);CHKERRA(ierr)
>
> ! end modified code
>
> call PCSetType(pc,PCMAT,ierr);CHKERRA(ierr)
>
>
>
> Is this the correct thing to do here?
>
Yes.
Matt
> Again, thank you for your help.
> Kenneth
>
>
> On Aug 7, 2018, at 11:13 AM, Jose E. Roman <jroman at dsic.upv.es> wrote:
>
> Currently there is no straightforward way to do this in SLEPc. For a
> preconditioned eigensolver such as JD you could use this:
> https://urldefense.proofpoint.com/v2/url?u=http-
> 3A__slepc.upv.es_documentation_current_docs_manualpages_ST_STPrecondSetMatForPC.html&d=DwIFaQ&c=imBPVzF25OnBgGmVOlcsiEgHoG1i6YHLR0Sj_gZ4adc&r=HM_IVChvU7XXVyLZB9xFCbw5_4LZNuzGEpVw5dUCyyo&m=L-kD1fTipDH_5BEC-T3FTfZKFQBUO6_ymvkpAiFfe38&s=kZTO6rmsUygc5qVZm5IkYMfjEMPVXloxoE0xi16ccUQ&e=
> 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)
>
>
>
>
--
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener
https://www.cse.buffalo.edu/~knepley/ <http://www.caam.rice.edu/~mk51/>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20180808/dfd56000/attachment.html>
More information about the petsc-users
mailing list