[Trilinos-Users] Smallest Eigenvalues

[Pate Motter]

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 [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 [mailto: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 [mailto: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 [mailto: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 [mailto:Trilinos-Users at software.sandia.gov]
https://software.sandia.gov/mailman/listinfo/trilinos-users [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/eff3afca/attachment.html>