[Trilinos-Users] PyTrilinos + Anasazi - Generalized eigenproblem

[Marko Budisic]

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
> http://software.sandia.gov/mailman/listinfo/trilinos-users