[Trilinos-Users] PyTrilinos + Anasazi - Generalized eigenproblem

[Marko Budisic]

Dear Trilinos list,

I have written a few days ago concerning my venture into using Anasazi 
from Python to compute solutions of eigenproblems. I am still using the 
same example as before, only I have modified the code to compute a 
generalized eigenproblem instead of manually scaling the matrix. The 
eigensolver used is BlockKrylovSchur.

In essence, I have a real matrix A (5x5) and a matrix D which is formed 
as a diagonal matrix that has row sums of A on the diagonal. I would 
like to compute Nev largest (in module or in real part) eigenvalues 
(lambda) and corresponding eigenvectors (x) using the generalized 
eigenproblem setup
Ax = lambda * Dx.
This is equivalent to computing the eigenproblem
(inv(D)*A )x = lambda * x
for the matrices considered in my code.

In the example used, matrix A is nonsingular with 3 real eigenvalues and 
a complex conjugate eigenvalue pair (the code is at the bottom of this 
e-mail).

My questions regarding PyTrilinos + Anasazi are as follows:

1) It seems that the setM() operation is ignored by the code, i.e. 
setting setM(D) still computes just the regular eigenproblem, instead of 
  the generalized one.

2) Whenever I try to set Nev > 1, e.g. Nev = 2, the solver complains 
about orthogonality constraints producing the following error:

"Anasazi::BlockKrylovSchurSolMgr: Potentially impossible orthogonality 
requests. Reduce basis size."

I am guessing there are some issues with my setting the blocksize. I am 
passing a randomized storage multivector to setInitVec(). I have tried 
setting the size of multivector to both 1 and to Nev (the number of 
eigenvalues requested) and in both cases I get the above complaint.

Could anyone shed some light on how one would go about figuring the 
right blocksize for this problem and in general?

3) I am setting only "Block Size", "Which" and "Verbosity" options for 
the parameter list passed to BlockKrylovSchurSolMgr. Should other 
options relating to blocksize be set too? Should I set even the "Block 
Size" manually?

4) Should AuxVec also be set prior to solving the problem/is there a 
benefit for doing so? If yes, how would one go about choosing the size 
of this multivector and should it be initialized in a certain way?

In the end, I am including my python code for this problem. At the very 
end, numpy.linalg is used to verify the results obtained by the Anasazi.

Any help is appreciated.

Best,
Marko Budisic

***********************************************
"""
Eigensolver test.
"""

import PyTrilinos
from PyTrilinos import Epetra
from PyTrilinos import Anasazi
from PyTrilinos import Teuchos
import numpy
from numpy import linalg

nx = 5
numglobals = nx
mycomm = Epetra.SerialComm()
mymap = Epetra.Map( numglobals, 0, mycomm )

A = numpy.array( [ [1,2,3,0,2],
                    [0,1,2,0,0],
                    [0,0,1,2,0],
                    [3,0,0,1,0],
                    [0,0,0,2,1] ] )

D = numpy.sum( A, axis=1 )

# populate matrices
W = Epetra.CrsMatrix(Epetra.Copy, mymap, 4)
Dsum = Epetra.CrsMatrix(Epetra.Copy, mymap, 1)

for i in xrange(A.shape[0]):
     boolnz = numpy.logical_not(A[i,:] == 0) # nonzeros
     indnz = numpy.arange( A.shape[1] )[boolnz]
     W.InsertGlobalValues( i, (A[i, boolnz]).tolist(), indnz.tolist() )
     Dsum.InsertGlobalValues( i, [ D[i] ], [i] )

Dsum.FillComplete()
W.FillComplete()

print W
print Dsum

### Create the eigenproblem
Nev = 2
blocksize = Nev
eigs = Anasazi.BasicEigenproblem()
eigs.setNEV( Nev )
eigs.setA( W )
eigs.setM( Dsum ) # get evectors of a normalized matrix W (or A)
eigs.setHermitian(False)

storage = Epetra.MultiVector( W.RowMap() , blocksize)
storage.Random() # initial guess must be nonzero
eigs.setInitVec(storage)
print "Problem set:", eigs.setProblem()

### Create the eigensolver
myPL = Teuchos.ParameterList()
myPL.set("Which", "LM")
myPL.set("Verbosity", Anasazi.Errors | Anasazi.Warnings | \
              Anasazi.FinalSummary )
myPL.set("Block Size", blocksize)

mngr = Anasazi.BlockKrylovSchurSolMgr( eigs, myPL )

# solve the eigenproblem
info = mngr.solve()

print "Problem solved? ", info == Anasazi.Converged

sol = eigs.getSolution()
evecs = sol.Evecs().ExtractCopy().T # convert to column-vectors
evals = sol.Evals()

print "Anasazi computed evals:", evals

##### Sanity check
lhand = numpy.dot( A, evecs )
rhand = numpy.dot( numpy.diag(D), evecs )

# this should yield evalues in columns when setM() is executed
print "Generalized eigenproblem check"
print lhand/rhand

# this should yield evalues in columns when setM() is not executed
print "Eigenproblem check"
print lhand/evecs

# numpy solution
# matrix A is left-normalized by row-sums because
# numpy.linalg doesn't have a generalized eigenproblem solver
Asc = numpy.dot( numpy.diag(1.0/D), A )
w,v = linalg.eig(Asc)

print "Numpy generalized evectors:"
print v
print "Numpy generalized evalues:"
print w

*************************************************