[Trilinos-Users] PyTrilinos + Anasazi - Generalized eigenproblem

[Marko Budisic]

Hi Chris,

yes, after I rescaled the auxiliary vector to unit norm, everything 
seems to work just fine on the larger problems, even without changing 
the "Orthogonalization" setting.

I believe I am set right now. If I run into some other problems, I will 
e-mail to the list. Thanks a lot once again for all your help.

Best,
Marko

Chris Baker wrote:
> Hi Marko,
> 
> It looks like the problem where the solver would freeze with an 
> auxiliary vector is some combination of the orthogonalization routine 
> and the small problem size. I suggest you start running with a larger 
> sample problem, at least 10 x 10. For the case that the solver seems to 
> freeze, it should be possible to fix this via:
> myPL.set("Orthogonalization","DGKS")
> 
> I played around on some row-stochastic matrices to verify that deflating 
> against the (normalized) ones() vector worked.
> 
> Thanks,
> Chris
> 
> On Mon, Dec 22, 2008 at 13:28, Marko Budisic 
> <mbudisic at engineering.ucsb.edu <mailto:mbudisic at engineering.ucsb.edu>> 
> wrote:
> 
>     Hi Chris,
> 
>     thank you for all your responses, especially those over the weekend.
>     I wasn't trying to rush you in my last e-mail, I was just
>     summarizing all my questions because I wanted to make them clear at
>     the end of the e-mail so you didn't have to decode what I was
>     actually asking for in the end. The discussion with you was
>     extremely helpful with my understanding of the code and it's a
>     pleasure to have such responsive support on the other end.
> 
>     As for the problem, I resorted to solving the generalized
>     eigenproblem in the exact way you suggested and I'm glad to say that
>     it is working well. My original matrix A might end up being
>     symmetric in some cases too, so I will try the Davidson solver in
>     those cases.
> 
>     All of the things being said, it seems that it is only Monday and
>     that the prospects are good for my problems to be resolved, thanks
>     to you and Heidi in the past weeks, so I am glad that is the case.
> 
>     Best,
> 
>     Marko
> 
>     Chris Baker wrote:
> 
>         Hello Marko,
> 
>         I am currently looking into the reason that Anasazi+PyTrilinos
>         becomes unresponsive when using Auxilliary vectors.
> 
>            Am I right to assume that one use of AuxVecs would be to used for
>            removing "known" eigenvectors from the process? In my case, I
>         will
>            often have row-stochastic matrices. Most of the times, they will
>            have the largest eigenvalue at 1 with eigenspace being the vector
>            filled with ones (the diagonal of the vector space). Would
>         passing a
>            vector filled with ones ( via say PutScalar(1.0) ) to AuxVecs
>         (and
>            choosing "LM" as the sorting method) make the code return the
>         second
>            largest evalue/evector as the result (assuming NEV = 1)?
> 
> 
>         Yes.
> 
>          
>            In any case, I have tried working a bit more with the AuxVecs
>         to see
>            why the code was not responsive before. It seems that AuxVecs
>            vectors have to be of length 1
> 
> 
>         Also correct. In fact, for multiple auxiliary vectors, they must
>         form an orthonormal basis.
> 
>            In this case however, this is the last thing the code prints
>         to the
>            screen (even with all verbosity options turned on) but it still
>            continues to run at 100% CPU. I waited for some time but it never
>            exited or printed other messages so I ended up killing the
>         process.
> 
> 
>         Yeah, I'm currently investigating this. We don't see this
>         problem with Anasazi in general, so I am leaning toward some
>         Python issue.
>          
>            When I pass the unit length vector to AuxVecs, the code runs but
>            fails (and exits to command prompt) with
> 
>            Anasazi::SVQBOrthoManager::findBasis(): Orthogonalization
>            constraints not feasible
> 
> 
>         This is because the auxiliary vectors add another constraint on
>         the orthogonalization routine. Your 5x5 example just doesn't
>         have enough room to move around: 1 vector for the auxiliary
>         vector plus 3 vectors minimum for BKS plus 1 additional vector
>         because the problem is non-Hermitian (this gives us extra room
>         to keep from splitting a conjugate pair).
> 
>         Having said that, I don't believe this is the cause for the
>         problem with setAuxVecs(). I ran on a larger matrix and still
>         saw the freezing.
> 
>            To summarize, my main question is the one stated at the
>         beginning,
>            about the use ot AuxVecs, the rest of the e-mail was more of a
>            followup to your answer. Of course, I am still interested in an
>            answer to my previous question about why the generalized
>            eigenproblem acts as the regular eigenproblem in PyTrilinos
>            (explained in my previous e-mails).
> 
> 
>         Solving generalized eigenvalue problems with BlockKrylovSchur
>         requires a little more work. I neglected to tell you this in my
>         previous email, because my brain doesn't work properly on the
>         weekend. In short, the generalized eigenvalue problem must be
>         transformed into a standard eigenvalue problem to work with BKS.
>         (/This is not the case with the other Anasazi eigensolvers, but
>         they don't solve non-Hermitian problems, so they are not
>         relevant to you./)
> 
>         In your particular case (searching for the largest magnitude
>         eigenvalues of (W,Dsum)), construct a linear operator
>         Op = inv(Dsum) W
>         This transforms the eigenvalue problem to
>         inv(Dsum) W v = v lambda
>         The eigenvectors and eigenvalues are unchanged from the
>         generalized eigenvalue problem.
> 
>         Pass Op to the BasicEigenproblem class, either by setOp() or
>         setA(). If a matrix is passed via setM(), it is used for the
>         inner product defining orthogonality. The benefit of doing that
>         is to restore symmetry in Op and to handle certain special
>         cases. In your case, W was not symmetric, so that symmetry in
>         inv(Dsum) W cannot be easily restored. I would not pass anything
>         to setM().
> 
>         However, in your case, because Dsum is a diagonal matrix,
>         applying the operator inv(Dsum) W serves to implicitly scale the
>         entries of W on every application of the operator. In this case,
>         it is probably better to simply scale them once in situ, as you
>         were doing before, and consider simply the standard eigenvalue
>         problem.
> 
>         I hope this helps. Even more, I hope this is correct; after all,
>         it is only Monday.
> 
>         I'm still looking at the setAuxVecs() problem with
>         Anasazi+PyTrilinos. I will respond to you and the list when it
>         is diagnosed.
> 
>         Chris
> 
> 
>         ------------------------------------------------------------------------
> 
> 
>         _______________________________________________
>         Trilinos-Users mailing list
>         Trilinos-Users at software.sandia.gov
>         <mailto:Trilinos-Users at software.sandia.gov>
>         http://software.sandia.gov/mailman/listinfo/trilinos-users
> 
> 
> 
> 
>