[Trilinos-Users] PyTrilinos + Anasazi - Generalized eigenproblem

Marko Budisic mbudisic at engineering.ucsb.edu
Sat Dec 20 13:05:58 MST 2008

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 

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.

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,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] )


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)

storage = Epetra.MultiVector( W.RowMap() , blocksize)
storage.Random() # initial guess must be nonzero
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


More information about the Trilinos-Users mailing list