[Trilinos-Users] Smallest Eigenvalues
Alicia Klinvex
aklinvex at purdue.edu
Wed Mar 25 11:43:35 MDT 2015
Hello Pate,
Have you gotten it to work on any matrices? My advice is to try it with a
small matrix first to see if your driver works at all.
I'm afraid I can't tell you whether your driver is correct based on the
code snippet you provided. I don't see the part where you create the
Anasazi eigensolver. As far as I can tell, what you are doing is creating
a GMRES solver and solving a linear system. Is the linear solve failing?
When you say it fails, do you mean it says "Belos: # iterations, but did
not finish"?
This is an example which does what you seem to want:
http://trilinos.org/docs/dev/packages/anasazi/doc/html/BlockKrylovSchur_2BlockKrylovSchurEpetraExGenBelos_8cpp-example.html.
That example shows how to use the Anasazi eigensolver BKS with a Belos
linear solver and Ifpack preconditioner.
As always, let me know if any of this is unclear, and I'd be happy to
elaborate :-)
Best wishes,
Alicia
On Wed, Mar 25, 2015 at 12:58 PM, Pate Motter <pate.motter at colorado.edu>
wrote:
> Hi again Alicia,
>
> I am having some issues with getting my eigenvalue solver working. It
> compiles, but fails to converge.
>
> I am currently using the default parameters which could be affecting it
> somewhat. Since this will be run on a variety of matrices with no common
> properties, I am not sure how to best code it so that it works as
> generically as possible. This was pieced together from a few different
> Trilinos example codes since there was not one doing the exact mixture of
> Belos/Anasazi/Tpetra. So there may be some slight issues when I combined
> them.
>
> Thanks for the help,
> Pate
>
> RCP<MV> calcEigenValues(const RCP<const MAT> &A, std::string eigenType) {
>
> // Preconditioner
> const std::string precType = "ILUT";
> Teuchos::ParameterList plist; // leave empty to use default params
>
> // Setup preconditioner using Ifpack2
> typedef Ifpack2::Preconditioner<ST, LO, GO, NT> prec_type;
> RCP<prec_type> prec;
> Ifpack2::Factory factory;
> prec = factory.create(precType, A);
> prec->setParameters(plist);
> prec->initialize();
> prec->compute();
> // Create vectors
> RCP<OP> M = prec;
> RCP<MV> X = rcp (new MV (A->getDomainMap (), 1) );
> RCP<MV> B = rcp (new MV (A->getRangeMap(), 1) );
> B->randomize();
>
> // Solve using Belos
> RCP<Teuchos::ParameterList> solverParams;
> Belos::SolverFactory<ST, MV, OP> solverFactory;
> RCP<Belos::SolverManager<ST, MV, OP> > solver =
> solverFactory.create("GMRES", solverParams);
> typedef Belos::LinearProblem<ST, MV, OP> problem_type;
> RCP<problem_type> problem =
> rcp (new problem_type (A, X, B) );
> problem->setRightPrec(M);
> problem->setProblem();
> solver->setProblem(problem);
> Belos::ReturnType result = solver->solve();
> const int numIters = solver->getNumIters();
> if (result == Belos::Converged) {
> *fos << "Belos: " << numIters << " iteration(s) to finish" <<
> std::endl;
> } else {
> *fos << "Belos: " << numIters << " iteration(s), but did not finish"
> << std::endl;
> }
>
> On 3/17/2015 5:33:01 PM, Alicia Klinvex <aklinvex at purdue.edu> wrote:
> If your matrix has the potential to be indefinite, here are some
> additional options (although BKS will certainly work too):
>
> 1. Use TraceMin-Davidson, and ask it to compute the harmonic Ritz
> values. This is brand new, so if you try this and see any wonky behavior,
> let me know.
> 2. Use spectrum folding. This method has the potential to be very slow
> (since you are in essence squaring the condition number of your matrix),
> but it will allow you to use any eigensolver you wish. Example:
> http://trilinos.org/docs/dev/packages/anasazi/doc/html/LOBPCGCustomStatusTest_8cpp-example.html
>
> If you decide to stick with BKS and have trouble getting your driver to
> work, let us know :-)
>
> - Alicia
>
> On Tue, Mar 17, 2015 at 6:51 PM, Pate Motter <pate.motter at colorado.edu>
> wrote:
>
>> Hi Alicia,
>>
>> Thank you for the information. This is not guaranteed to be an SPD matrix
>> so the first option is probably the best. I remember seeing an example
>> going through that process, but could not get it to work either. I will
>> retry that though as it could have been an error elsewhere.
>>
>> -Pate
>>
>> On 3/17/2015 4:47:03 PM, Alicia Klinvex <aklinvex at purdue.edu> wrote:
>> 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/20150325/96802044/attachment.html>
More information about the Trilinos-Users
mailing list