<div dir="ltr"><div class="gmail_quote"><div dir="ltr">On Wed, Aug 8, 2018 at 11:51 AM Kenneth Hall, SC.D. <<a href="mailto:kenneth.c.hall@duke.edu">kenneth.c.hall@duke.edu</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<div style="word-wrap:break-word;line-break:after-white-space">
Roman,<br>
<br>
When I tried what you suggested, I get the following error:<br>
<br>
<font face="Courier"> call PCSetOperators(pc,PETSC_NULL_MAT,my_matmul_k,ierr);CHKERRA(ierr) <br>
1<br>
Error: Invalid procedure argument at (1)<br>
make: *** [/Users/hall/git/Mustang/Mustang/build_dir/slepc_eigen2.o] Error 1<br>
<br>
</font><br>
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:<br>
<br>
call EPSSetWhichEigenpairs(eps,EPS_TARGET_REAL,ierr);CHKERRA(ierr)<br>
call EPSSetTarget(eps,sigma,ierr);CHKERRA(ierr)<br>
call EPSGetST(eps,st,ierr);CHKERRA(ierr)<br>
call STSetType(st,STSINVERT,ierr);CHKERRA(ierr) <br>
call STGetKSP(st,ksp,ierr);CHKERRA(ierr) <br>
call KSPGetPC(ksp,pc,ierr);CHKERRA(ierr) <br>
call KSPSetType(ksp,KSPBCGS,ierr);CHKERRA(ierr) <br>
call PCSetType(pc,PCNONE,ierr);CHKERRA(ierr) <br>
call EPSSetUp(eps,ierr);CHKERRA(ierr) <br>
<br>
<div>! begin modified code</div>
<div><br>
call MatCreateShell(PETSC_COMM_WORLD,n,n,n,n,PETSC_NULL_INTEGER,P,ierr) ; CHKERRA(ierr)<br>
call MatShellSetOperation(P,MATOP_MULT,my_matmul_k,ierr); CHKERRA(ierr)<br>
call PCSetOperators(pc,PETSC_NULL_MAT,P,ierr);CHKERRA(ierr)</div>
<div><br>
! end modified code<br>
<br>
</div>
<div>call PCSetType(pc,PCMAT,ierr);CHKERRA(ierr)<br>
<br>
<br>
<br>
</div>
<div>Is this the correct thing to do here?</div></div></blockquote><div><br></div><div>Yes.</div><div><br></div><div> Matt</div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div style="word-wrap:break-word;line-break:after-white-space">
<div>Again, thank you for your help.<br>
<div>Kenneth<br>
<br>
<br>
<blockquote type="cite">On Aug 7, 2018, at 11:13 AM, Jose E. Roman <<a href="mailto:jroman@dsic.upv.es" target="_blank">jroman@dsic.upv.es</a>> wrote:<br>
<br>
Currently there is no straightforward way to do this in SLEPc. For a preconditioned eigensolver such as JD you could use this:<br>
<a href="https://urldefense.proofpoint.com/v2/url?u=http-" target="_blank">https://urldefense.proofpoint.com/v2/url?u=http-</a>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=<br>
in combination with pc=PCMAT.<br>
<br>
For a Krylov solver you will have to do something like this:<br>
<br>
call EPSSetWhichEigenpairs(eps,EPS_TARGET_REAL,ierr);CHKERRA(ierr)<br>
call EPSSetTarget(eps,sigma,ierr);CHKERRA(ierr)<br>
call EPSGetST(eps,st,ierr);CHKERRA(ierr)<br>
call STSetType(st,STSINVERT,ierr);CHKERRA(ierr)<br>
call STGetKSP(st,ksp,ierr);CHKERRA(ierr)<br>
call KSPGetPC(ksp,pc,ierr);CHKERRA(ierr)<br>
call KSPSetType(ksp,KSPBCGS,ierr);CHKERRA(ierr)<br>
call PCSetType(pc,PCNONE,ierr);CHKERRA(ierr) ! dummy preconditioner<br>
<br>
! provide the preconditioner matrix after setup<br>
call EPSSetUp(eps,ierr);CHKERRA(ierr)<br>
call PCSetOperators(pc,PETSC_NULL_MAT,my_matmul_K,ierr);CHKERRA(ierr)<br>
call PCSetType(pc,PCMAT,ierr);CHKERRA(ierr)<br>
<br>
We may change the interface in future versions so that this can be done in a more natural way.<br>
<br>
Jose<br>
<br>
<br>
<blockquote type="cite">El 7 ago 2018, a las 16:42, Kenneth Hall, SC.D. <<a href="mailto:kenneth.c.hall@duke.edu" target="_blank">kenneth.c.hall@duke.edu</a>> escribió:<br>
<br>
Hello,<br>
<br>
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<br>
<br>
A.x = lambda.B.x<br>
<br>
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<br>
<br>
K^-1 \approx (A - sigma.B)^-1<br>
<br>
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.<br>
<br>
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. <br>
<br>
Thanks in advance.<br>
Kenneth Hall<br>
<br>
<br>
===============================================================<br>
<br>
!<br>
!.... initialize slepc<br>
call SlepcInitialize(PETSC_NULL_CHARACTER,ierr) ; if (ierr .ne. 0) then ; stop 'SlepcInitialize failed' ; end if<br>
!<br>
!.... get size and rank of MPI process<br>
call MPI_Comm_size(PETSC_COMM_WORLD,sz ,ierr) ; CHKERRA(ierr)<br>
call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr) ; CHKERRA(ierr)<br>
!<br>
if (sz .ne. 1) then ; SETERRA(PETSC_COMM_SELF,1,'This is a uniprocessor example only!'); endif<br>
!<br>
!.... get user options<br>
call PetscOptionsGetInt (PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-n' ,n ,flg,ierr) ; CHKERRA(ierr)<br>
call PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-terse',terse,ierr) ; CHKERRA(ierr)<br>
!<br>
!.... write out option list<br>
if (rank .eq. 0 ) then<br>
write(*,*) 'n = ',n<br>
end if<br>
!<br>
!.... Register the matrix-vector subroutine for the operator that defines the eigensystem, Ax=kBx<br>
call MatCreateShell(PETSC_COMM_WORLD,n,n,n,n,PETSC_NULL_INTEGER,A,ierr) ; CHKERRA(ierr)<br>
call MatShellSetOperation(A,MATOP_MULT,my_matmul_A,ierr) ; CHKERRA(ierr)<br>
!<br>
!.... set up B matrix<br>
call MatCreate(PETSC_COMM_WORLD,B,ierr) ; CHKERRA(ierr)<br>
call MatSetSizes(B,n,n,n,n,ierr) ; CHKERRA(ierr)<br>
call MatSetFromOptions(B,ierr) ; CHKERRA(ierr)<br>
call MatSetUp(B,ierr) ; CHKERRA(ierr)<br>
call My_Matrix_B(B,n)<br>
!<br>
!.... set up eigenvectors<br>
call VecCreate(PETSC_COMM_WORLD,xr,ierr) ; CHKERRA(ierr)<br>
call VecSetSizes(xr,n,n,ierr) ; CHKERRA(ierr)<br>
call VecSetFromOptions(xr,ierr) ; CHKERRA(ierr)<br>
call VecSetUp(xr,ierr) ; CHKERRA(ierr)<br>
!<br>
call VecCreate(PETSC_COMM_WORLD,xi,ierr) ; CHKERRA(ierr)<br>
call VecSetSizes(xi,n,n,ierr) ; CHKERRA(ierr)<br>
call VecSetFromOptions(xi,ierr) ; CHKERRA(ierr)<br>
call VecSetUp(xi,ierr) ; CHKERRA(ierr)<br>
!<br>
!.... Create eigensolver context<br>
call EPSCreate(PETSC_COMM_WORLD,eps,ierr) ; CHKERRA(ierr)<br>
!<br>
!.... Set operators. In this case, it is a generalized eigenvalue problem<br>
call EPSSetOperators(eps,A,B,ierr) ; CHKERRA(ierr)<br>
call EPSSetProblemType(eps,EPS_GNHEP,ierr) ; CHKERRA(ierr)<br>
!<br>
!.... set default options<br>
call EPSSetType(eps,EPSKRYLOVSCHUR,ierr) ; CHKERRA(ierr)<br>
!<br>
!.... set default eigenvalues to find<br>
call EPSSetWhichEigenpairs(eps,EPS_LARGEST_REAL,ierr) ; CHKERRA(ierr)<br>
!<br>
!.... number of requested eigenvalues<br>
call EPSSetDimensions(eps,nev,PETSC_DEFAULT_INTEGER,PETSC_DEFAULT_INTEGER,ierr) ; CHKERRA(ierr)<br>
!<br>
!.... set stopping criteria<br>
call EPSSetTolerances(eps,tol,maxit,ierr) ; CHKERRA(ierr)<br>
!<br>
!.... set up monitors<br>
call PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_DEFAULT,vf,ierr) ; CHKERRA(ierr)<br>
call EPSMonitorSet(eps,EPSMONITORFIRST,vf,PetscViewerAndFormatDestroy,ierr) ; CHKERRA(ierr)<br>
call SlepcConvMonitorCreate(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_DEFAULT,ctx,ierr) ; CHKERRA(ierr)<br>
call EPSMonitorSet(eps,EPSMONITORCONVERGED,ctx,SlepcConvMonitorDestroy,ierr) ; CHKERRA(ierr)<br>
!<br>
!.... Set solver parameters at runtime<br>
call EPSSetFromOptions(eps,ierr) ; CHKERRA(ierr)<br>
!<br>
!.... Solve the eigensystem<br>
call EPSSolve(eps,ierr) ; CHKERRA(ierr)<br>
<br>
</blockquote>
<br>
</blockquote>
<br>
</div>
</div>
</div>
</blockquote></div><br clear="all"><div><br></div>-- <br><div dir="ltr" class="gmail_signature" data-smartmail="gmail_signature"><div dir="ltr"><div><div dir="ltr"><div>What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>-- Norbert Wiener</div><div><br></div><div><a href="http://www.caam.rice.edu/~mk51/" target="_blank">https://www.cse.buffalo.edu/~knepley/</a><br></div></div></div></div></div></div>