Anasazi  Version of the Day
LOBPCGEpetra.cpp

Use LOBPCG with Epetra test problem from Galeri.This example computes the eigenvalues of largest magnitude of an eigenvalue problem $A x = x$, using Anasazi's implementation of the LOBPCG method, with Epetra linear algebra. It uses the Galeri package to construct the test problem.

// ***********************************************************************
//
// Anasazi: Block Eigensolvers Package
//
// Under terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
// the U.S. Government retains certain rights in this software.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are
// met:
//
// 1. Redistributions of source code must retain the above copyright
// notice, this list of conditions and the following disclaimer.
//
// 2. Redistributions in binary form must reproduce the above copyright
// notice, this list of conditions and the following disclaimer in the
// documentation and/or other materials provided with the distribution.
//
// 3. Neither the name of the Corporation nor the names of the
// contributors may be used to endorse or promote products derived from
// this software without specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
//
// ***********************************************************************
// Include header for LOBPCG eigensolver
// Include header to define eigenproblem Ax = \lambda*x
// plan to use Tpetra objects instead of Epetra objects, include
// Thyra objects instead of Epetra objects.
// Include header for Epetra sparse matrix and multivector.
#include "Epetra_CrsMatrix.h"
#include "Epetra_MultiVector.h"
// The Trilinos package Galeri has many example problems.
#include "Galeri_Maps.h"
#include "Galeri_CrsMatrices.h"
// Include selected communicator class required by Epetra objects
#ifdef EPETRA_MPI
# include "Epetra_MpiComm.h"
#else
# include "Epetra_SerialComm.h"
#endif // EPETRA_MPI
int
main (int argc, char *argv[])
{
using Teuchos::RCP;
using Teuchos::rcp;
using std::cerr;
using std::cout;
using std::endl;
// Anasazi solvers have the following template parameters:
//
// - Scalar: The type of dot product results.
// - MV: The type of (multi)vectors.
// - OP: The type of operators (functions from multivector to
// multivector). A matrix (like Epetra_CrsMatrix) is an example
// of an operator; an Ifpack preconditioner is another example.
//
// Here, Scalar is double, MV is Epetra_MultiVector, and OP is
// Epetra_Operator.
typedef Epetra_MultiVector MV;
typedef Epetra_Operator OP;
// Anasazi's interface to the linear algebra implementation's
// vectors and multivectors is MultiVecTraits. This is a C++ traits
// interface. It is templated on Scalar and MV (see above).
#ifdef EPETRA_MPI
MPI_Init (&argc, &argv);
Epetra_MpiComm Comm (MPI_COMM_WORLD);
#else
#endif // EPETRA_MPI
const int MyPID = Comm.MyPID ();
//
// Set up the test problem.
//
// We use Trilinos' Galeri package to construct a test problem.
// Here, we use a discretization of the 2-D Laplacian operator.
// The global mesh size is nx * nx.
//
const int nx = 30;
GaleriList.set ("n", nx * nx);
GaleriList.set ("nx", nx);
GaleriList.set ("ny", nx);
RCP<Epetra_Map> Map = rcp (Galeri::CreateMap ("Linear", Comm, GaleriList));
RCP<Epetra_RowMatrix> A =
rcp (Galeri::CreateCrsMatrix ("Laplace2D", &*Map, GaleriList));
// Set eigensolver parameters.
const double tol = 1.0e-8; // convergence tolerance
const int nev = 10; // number of eigenvalues for which to solve
const int blockSize = 5; // block size (number of eigenvectors processed at once)
const int maxIters = 500; // maximum number of iterations
// Create a set of initial vectors to start the eigensolver.
// This needs to have the same number of columns as the block size.
RCP<MV> ivec = rcp (new MV (*Map, blockSize));
ivec->Random ();
// Create the eigenproblem. This object holds all the stuff about
// your problem that Anasazi will see. In this case, it knows about
// the matrix A and the inital vectors.
RCP<Anasazi::BasicEigenproblem<double, MV, OP> > problem =
// Tell the eigenproblem that the operator A is symmetric.
problem->setHermitian (true);
// Set the number of eigenvalues requested
problem->setNEV (nev);
// Tell the eigenproblem that you are finishing passing it information.
const bool boolret = problem->setProblem();
if (boolret != true) {
if (MyPID == 0) {
cerr << "Anasazi::BasicEigenproblem::setProblem() returned an error." << endl;
}
#ifdef EPETRA_MPI
MPI_Finalize ();
#endif // EPETRA_MPI
return -1;
}
// Create a ParameterList, to pass parameters into the LOBPCG
// eigensolver.
anasaziPL.set ("Which", "SM");
anasaziPL.set ("Block Size", blockSize);
anasaziPL.set ("Maximum Iterations", maxIters);
anasaziPL.set ("Convergence Tolerance", tol);
anasaziPL.set ("Full Ortho", true);
anasaziPL.set ("Use Locking", true);
anasaziPL.set ("Verbosity", Anasazi::Errors + Anasazi::Warnings +
// Create the LOBPCG eigensolver.
Anasazi::LOBPCGSolMgr<double, MV, OP> anasaziSolver (problem, anasaziPL);
// Solve the eigenvalue problem.
//
// Note that creating the eigensolver is separate from solving it.
// After creating the eigensolver, you may call solve() multiple
// times with different parameters or initial vectors. This lets
// you reuse intermediate state, like allocated basis vectors.
Anasazi::ReturnType returnCode = anasaziSolver.solve ();
if (returnCode != Anasazi::Converged && MyPID == 0) {
cout << "Anasazi eigensolver did not converge." << endl;
}
// Get the eigenvalues and eigenvectors from the eigenproblem.
Anasazi::Eigensolution<double,MV> sol = problem->getSolution ();
// Anasazi returns eigenvalues as Anasazi::Value, so that if
// Anasazi's Scalar type is real-valued (as it is in this case), but
// some eigenvalues are complex, you can still access the
// eigenvalues correctly. In this case, there are no complex
// eigenvalues, since the matrix pencil is symmetric.
std::vector<Anasazi::Value<double> > evals = sol.Evals;
RCP<MV> evecs = sol.Evecs;
// Compute residuals.
std::vector<double> normR (sol.numVecs);
if (sol.numVecs > 0) {
MV tempAevec (*Map, sol.numVecs);
T.putScalar (0.0);
for (int i=0; i<sol.numVecs; ++i) {
T(i,i) = evals[i].realpart;
}
A->Apply (*evecs, tempAevec);
MVT::MvTimesMatAddMv (-1.0, *evecs, T, 1.0, tempAevec);
MVT::MvNorm (tempAevec, normR);
}
// Print the results on MPI process 0.
if (MyPID == 0) {
cout << endl << "Solver manager returned "
<< (returnCode == Anasazi::Converged ? "converged." : "unconverged.")
<< endl << endl
<< "------------------------------------------------------" << endl
<< std::setw(16) << "Eigenvalue"
<< std::setw(18) << "Direct Residual"
<< endl
<< "------------------------------------------------------" << endl;
for (int i=0; i<sol.numVecs; ++i) {
cout << std::setw(16) << evals[i].realpart
<< std::setw(18) << normR[i] / evals[i].realpart
<< endl;
}
cout << "------------------------------------------------------" << endl;
}
#ifdef EPETRA_MPI
MPI_Finalize () ;
#endif // EPETRA_MPI
return 0;
}