[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