[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
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
*************************************************
More information about the Trilinos-Users
mailing list