[Trilinos-Users] Smallest Eigenvalues

Alicia Klinvex aklinvex at purdue.edu
Tue Mar 17 16:47:01 MDT 2015


Hello Pate,

If you want to find the smallest eigenvalues, I'd recommend one of the
following options:
1.  Using BKS in shift-and-invert mode.  There are several examples
demonstrating how to do this on the Anasazi page.  If you want to go this
route, you will have to create a linear solver (either direct or iterative)
and pass it to the BKS solver manager.
2.  Pick a different eigensolver.  Is your matrix symmetric positive
definite?  If so, you can use LOBPCG, RTR, or TraceMin-Davidson.  If your
matrix is not SPD, your choices are more limited, and we can discuss your
options further.

Let me know which route you want to go, and I'll try to help you with it :-)

Best wishes,
Alicia

On Tue, Mar 17, 2015 at 6:16 PM, Pate Motter <pate.motter at colorado.edu>
wrote:

> 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;
> }
>
>
> _______________________________________________
> Trilinos-Users mailing list
> Trilinos-Users at software.sandia.gov
> https://software.sandia.gov/mailman/listinfo/trilinos-users
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://software.sandia.gov/pipermail/trilinos-users/attachments/20150317/bd4f9516/attachment.html>


More information about the Trilinos-Users mailing list