[petsc-users] EPSSolve finds eigensolutions in deflation space

[Jose E. Roman]

Giacomo:
The custom function for stopping based on a user-defined criterion will be available in slepc-3.7, to be released in a few days. You can already try this functionality in the master branch of the development repo. Hope it fits your needs.
Regarding the deflation space issue, I will try your example and get back to you soon.
Jose


> El 8 may 2016, a las 13:20, Giacomo Mulas <giacomo.mulas84 at gmail.com> escribió:
> 
> I found a problem in the default iterative eigensystem solver of SLEPC,
> namely the Krilov-Schur one: it does not respect the deflation space that it
> is given.  This makes it very inefficient even in a case in which it should
> be extremely effective.
> 
> The problem I want to solve: I have a symmetric operator (an approximated
> vibrational Hamiltonian, but that's unimportant) and I want to selectively
> find a subset of eigenvectors so that some predefined vectors (which I call
> "target" vectors) are very nearly contained in the eigenvector subspace.
> 
> I therefore have both a well defined selection criterion (find eigenstates
> which have maximum projection on the subspace of "target" vectors) and a
> well defined condition for stopping (say, >99% of each of the "target"
> vectors is contained in the subspace of the eigenstates found).
> 
> The first one can be explicitly cast in a function and given to SLEPC via
> EPSSetArbitrarySelection and EPSSetWhichEigenpairs, so that the solver is
> guided to the solutions I am most interested in.  I also make use of
> EPSSetInitialSpace to have the solver start from a guess which is "rich" in
> the "target" vectors (just the "target" vector which has the smallest
> projection on the eigenstate space so far).
> 
> As to stopping when I found just enough eigenstates, unfortunately it is
> thus far impossible to give a custom function to slepc so that it will keep
> finding eigenvectors until all are found or the user-defined function tells
> it it's enough.  Therefore, my second best solution was to repeatedly run
> EPSSolve, aiming at a small number of eigensolutions in each call, putting
> previously found solutions in the deflation space and passing it to the
> solver via EPSSetDeflationSpace.  Unfortunately, this did not simply work
> out of the box.
> 
> The first problem is that after the first few, easiest to find
> eigensolutions are found, EPSSolve frequently does not find any solutions,
> if run with a very conservative EPSSetDimensions (e.g.  like "find me just
> one eigensolution").  I thus wrote my code so that the number of requested
> solutions in EPSSolve are increased if no new solutions are found by an
> EPSSolve call, and decreased if too many were found in just one call. Inefficient as it may be, due to many useless EPSSolve calls, it works.
> 
> However, when trying to use this on real cases, another problem arose, the
> big one: frequently, EPSSolve would find and return eigenvectors that were
> in the deflation space!  Some eigensolutions would therefore appear more
> than once in the set of solutions, of course breaking everything down the
> road of my physics calculation before I figured out why. When I found out, I set up a check in my code that explicitly computes the
> projection of newly found eigenstates on previously found ones, and discards
> them if they are duplicates.  But this is ineffective for several reasons:
> first, such a check is a waste of time and should not be necessary,
> otherwise what is the point in setting up a deflation space?  second, and
> most important, I have to run EPSSolve looking for lots of solutions to only
> find one or two really new ones, because all the first ones are duplicates. Even if the initial guess is rich in the solution I actually want, and poor
> in the duplicate ones. As a result, even if I would just need to find a
> (relatively) small subset of the eigenvectors, lapack turns out to be always
> much faster, which should not be the case.
> 
> I am attaching a test code that shows this behaviour, in the hope that it
> will help the SLEPC developers in finding why the deflation space is not
> respected and, if possible, correct it.  I know that being able to reproduce
> a problem is the first, sometimes actually the most important, step to
> attempt to solve it, so there it goes.  Of course I also have a tiny hope
> that they might actually include support for a user-defined check for
> "enough solutions" in EPSSolve...  There is in fact something not too
> different already implemented, for EPSSolve can find all solutions with
> eigenvalues in a given interval, without knowing in advance how many they
> will be, without the user having to loop over many tentative EPSSolve calls
> and finding lots of unneeded additional solutions in the process.
> 
> To run the test code, just compile the source, and run it in the same
> directory of the two data files.  This is a small test case, so that it runs
> fast enough, but I have much larger ones and can produce them more or less
> as large as wanted on request.  Upon running the code, you will see that
> lots of tentative EPSSolve calls are made, many times finding no solutions
> or only solutions already in the deflation space (which contains previously
> found eigenvectors).
> In this code I privileged clarity over efficiency, no attempt was made to
> overoptimise anything.  I put lots of comments in it, so that hopefully it
> is self-explanatory.  But of course I am more than willing to explain
> whatever is not clear in it.  Just ignore the part that reads the input
> files, that's irrelevant.
> 
> Of course, I would also very much appreciate it if I could get some clever
> suggestions, from SLEPC developers or anyone else on the list, as to how to
> solve my problem more efficiently.  And of course I am perfectly willing
> to thank them and/or to include them as authors in the scientific publication
> describing my molecular spectroscopy code, if some substantial contribution
> is offered.
> 
> Best regards
> Giacomo Mulas
> 
> -- 
> _________________________________________________________________
> 
> Giacomo Mulas <gmulas at oa-cagliari.inaf.it>
> _________________________________________________________________
> 
> INAF - Osservatorio Astronomico di Cagliari
> via della scienza 5 - 09047 Selargius (CA)
> 
> tel.   +39 070 71180244
> mob. : +39 329  6603810
> _________________________________________________________________
> 
> "When the storms are raging around you, stay right where you are"
>                         (Freddy Mercury)
> _________________________________________________________________<slepctesttargets.dat><slepctestham.dat><slepctester2.c>