[Trilinos-Users] Smallest Eigenvalues
Pate Motter
pate.motter at colorado.edu
Tue Mar 17 16:16:21 MDT 2015
Hi,
I am having difficulties converging when finding the smallest magnitude
eigenvalues in a sparse, square, real matrix. I am trying to use a
combination of Tpetra and Anasazi to accomplish this based on a few of the
examples I was able to find online. It currently works for finding largest
magnitude.
Thanks,
Pate Motter
RCP<MV> calcEigenValues(const RCP<MAT> &A, std::string eigenType) {
// Get norm
ST mat_norm = A->getFrobeniusNorm();
// Eigensolver parameters
int nev = 1;
int blockSize = 1;
int numBlocks = 10*nev / blockSize;
ST tol = 1e-6;
// Create parameters to pass to the solver
Teuchos::ParameterList MyPL;
MyPL.set("Block Size", blockSize );
MyPL.set("Convergence Tolerance", tol);
MyPL.set("Num Blocks", numBlocks);
// Default to largest magnitude
if (eigenType.compare("SM") == 0) {
MyPL.set("Which", "SM");
} else if (eigenType.compare("SR") == 0) {
MyPL.set("Which", "SR");
} else if (eigenType.compare("LR") == 0) {
MyPL.set("Which", "LR");
} else {
MyPL.set("Which", "LM");
}
// Create multivector for a initial vector to start the solver
RCP<MV> ivec = rcp (new MV(A->getRowMap(), blockSize));
MVT::MvRandom(*ivec);
// Create eigenproblem
RCP<Anasazi::BasicEigenproblem<ST, MV, OP> > MyProblem =
Teuchos::rcp(new Anasazi::BasicEigenproblem<ST, MV, OP>(A, ivec));
MyProblem->setHermitian(false);
MyProblem->setNEV(nev);
MyProblem->setProblem();
Anasazi::BlockKrylovSchurSolMgr<ST, MV, OP> MySolverMgr(MyProblem, MyPL);
// Solve the problem
Anasazi::ReturnType returnCode = MySolverMgr.solve();
if (returnCode != Anasazi::Converged) {
*fos << "unconverged" << ", ";
}
// Get the results
Anasazi::Eigensolution<ST, MV> sol = MyProblem->getSolution();
std::vector<Anasazi::Value<ST> > evals = sol.Evals;
RCP<MV> evecs = sol.Evecs;
int numev = sol.numVecs;
// Compute residual
if (numev > 0) {
Teuchos::SerialDenseMatrix<int, ST> T(numev,numev);
for (int i = 0; i < numev; i++) {
T(i,i) = evals[i].realpart;
}
std::vector<ST> normR(sol.numVecs);
MV Kvec(A->getRowMap(), MVT::GetNumberVecs(*evecs));
OPT::Apply(*A, *evecs, Kvec);
MVT::MvTimesMatAddMv(-1.0, *evecs, T, 1.0, Kvec);
MVT::MvNorm(Kvec, normR);
for (int i=0; i<numev; i++) {
*fos << evals[i].realpart << ", " << normR[i]/mat_norm << ", ";
}
}
return evecs;
}
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://software.sandia.gov/pipermail/trilinos-users/attachments/20150317/fe4b4b62/attachment.html>
More information about the Trilinos-Users
mailing list