Anasazi Version of the Day
Loading...
Searching...
No Matches
BlockKrylovSchur/BlockKrylovSchurEpetraExGenAztecOO.cpp

Use Anasazi::BlockKrylovSchurSolMgr to solve a generalized eigenvalue problem, using Epetra data stuctures and the AztecOO package of iterative linear solvers and preconditioners.

Use Anasazi::BlockKrylovSchurSolMgr to solve a generalized eigenvalue problem, using Epetra data stuctures and the AztecOO package of iterative linear solvers and preconditioners.

// @HEADER
// ***********************************************************************
//
// Anasazi: Block Eigensolvers Package
// Copyright 2004 Sandia Corporation
//
// 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)
//
// ***********************************************************************
// @HEADER
// This example computes the eigenvalues of smallest magnitude of a
// generalized eigenvalue problem $K x = \lambda M x$, using Anasazi's
// implementation of the block Krylov-Schur method. It implements
// inverse iteration using an AztecOO iterative linear solver with an
// Ifpack preconditioner.
//
// Anasazi computes the smallest-magnitude eigenvalues using a
// shift-and-invert strategy. For simplicity, the code below uses a
// shift of zero. It illustrates the general pattern for using
// Anasazi for this problem:
//
// 1. Construct an "operator" A such that $Az = K^{-1} M z$.
// 2. Use Anasazi to solve $Az = \sigma z$, which is a spectral
// transformation of the original problem $K x = \lambda M x$.
// 3. The eigenvalues $\lambda$ of the original problem are the
// inverses of the eigenvalues $\sigma$ of the transformed
// problem.
//
// The "operator A such that $A z = K^{-1} M z$" is a subclass of
// Epetra_Operator. The Apply method of that operator takes the
// vector b, and computes $x = K^{-1} M b$. It does so first by
// applying the matrix M, and then by solving the linear system $K x =
// M b$ for x. Trilinos implements many different ways to solve
// linear systems. The example below uses an iterative linear solver
// provided by the AztecOO package, with an Ifpack preconditioner, to
// solve linear systems.
// Include header for block Krylov-Schur solver
// Include header to define basic eigenproblem Ax = \lambda*Bx
// Include header to provide Anasazi with Epetra adapters. If you
// plan to use Tpetra objects instead of Epetra objects, include
// AnasaziTpetraAdapter.hpp instead; do analogously if you plan to use
// Thyra objects instead of Epetra objects.
// Include header for Epetra sparse matrix, Map (representation of
// parallel distributions), and linear problem. AztecOO uses the
// latter to encapsulate linear systems to solve.
#include "Epetra_Map.h"
#include "Epetra_CrsMatrix.h"
#include "Epetra_LinearProblem.h"
// Include header for AztecOO iterative linear solver, and
// AztecOO_Operator. The latter wraps an AztecOO solver in an
// Epetra_Operator.
#include "AztecOO.h"
#include "AztecOO_Operator.h"
// Include header for Ifpack preconditioner factory
#include "Ifpack.h"
// Include header for Teuchos serial dense matrix, which we use to
// compute the eigenvectors.
#include "Teuchos_SerialDenseMatrix.hpp"
// Include header for the problem definition
#include "ModeLaplace2DQ2.h"
// Include selected communicator class and map required by Epetra objects
#ifdef EPETRA_MPI
# include "Epetra_MpiComm.h"
#else
# include "Epetra_SerialComm.h"
#endif // EPETRA_MPI
// ****************************************************************************
// BEGIN MAIN ROUTINE
// ****************************************************************************
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;
#ifdef EPETRA_MPI
// Initialize MPI
MPI_Init (&argc, &argv);
Epetra_MpiComm Comm (MPI_COMM_WORLD);
#else
Epetra_SerialComm Comm;
#endif // EPETRA_MPI
const int MyPID = Comm.MyPID ();
//
// Set up the test problem
//
// Dimensionality of the spatial domain to discretize
const int space_dim = 2;
// Size of each of the dimensions of the (discrete) domain
std::vector<double> brick_dim (space_dim);
brick_dim[0] = 1.0;
brick_dim[1] = 1.0;
// Number of elements in each of the dimensions of the domain
std::vector<int> elements (space_dim);
elements[0] = 10;
elements[1] = 10;
// Create the test problem
RCP<ModalProblem> testCase =
rcp (new ModeLaplace2DQ2 (Comm, brick_dim[0], elements[0],
brick_dim[1], elements[1]));
// Get the stiffness and mass matrices.
//
// rcp (T*, false) returns a nonowning (doesn't deallocate) RCP.
RCP<Epetra_CrsMatrix> K =
rcp (const_cast<Epetra_CrsMatrix* > (testCase->getStiffness ()), false);
RCP<Epetra_CrsMatrix> M =
rcp (const_cast<Epetra_CrsMatrix*> (testCase->getMass ()), false);
//
// Create linear solver for linear systems with K
//
// Anasazi uses shift and invert, with a "shift" of zero, to find
// the eigenvalues of least magnitude. In this example, we
// implement the "invert" part of shift and invert by using an
// AztecOO iterative linear solver with an Ifpack preconditioner.
//
//
// Construct Ifpack preconditioner
//
// An Ifpack "factory" knows how to create Ifpack preconditioners.
Ifpack Factory;
// Use the Factory to create the preconditioner. Factory.Create()
// returns a raw Ifpack2_Preconditioner pointer. The caller (that's
// us!) is responsible for deallocation, so we use an owning RCP to
// deallocate it automatically.
//
// The Factory needs a string name of the preconditioner, and the
// overlap level. (Almost) all Ifpack preconditioners are local to
// each MPI process. Ifpack automatically does an additive Schwarz
// decomposition across processes. The default overlap level for
// additive Schwarz is zero, but you can use a different overlap
// level here if you want. The overlap level must be nonnegative,
// and only matters if running with more than one MPI process.
std::string PrecType = "ICT"; // incomplete Cholesky
const int OverlapLevel = 0;
// Create the preconditioner.
RCP<Ifpack_Preconditioner> Prec =
rcp (Factory.Create (PrecType, &*K, OverlapLevel));
if (Prec.is_null ()) {
throw std::logic_error ("Ifpack's factory returned a NULL preconditioner!");
}
//
// Set Ifpack preconditioner parameters.
//
// Set parameters after creating the preconditioner. Please refer
// to Ifpack's documentation for a list of valid parameters.
//
Teuchos::ParameterList ifpackList;
ifpackList.set ("fact: drop tolerance", 1e-4);
ifpackList.set ("fact: ict level-of-fill", 0.0);
// The "combine mode" describes how to combine contributions from
// other MPI processes. It only matters if the overlap level is
// nonzero. See Epetra_CombineMode.h for documentation of of the
// various combine modes.
ifpackList.set ("schwarz: combine mode", "Add");
// Set the parameters.
IFPACK_CHK_ERR(Prec->SetParameters(ifpackList));
// Initialize the preconditioner. This only looks at the structure
// of the matrix, not the values. Nevertheless, the matrix must
// generally be fill complete at this point. Compare Initialize()
// to a symbolic factorization, and Compute() (see below) to a
// numeric factorization.
IFPACK_CHK_ERR(Prec->Initialize());
// Compute the preconditioner. This looks at both the structure and
// the values of the matrix. Compare Initialize() (see above) to a
// symbolic factorization, and Compute() (see below) to a numeric
// factorization.
IFPACK_CHK_ERR(Prec->Compute());
//
// Set up AztecOO iterative solver for solving linear systems with K.
//
// Create Epetra linear problem class for solving linear systems
// with K. This implements the inverse operator for shift and
// invert.
Epetra_LinearProblem precProblem;
// Tell the linear problem about the matrix K. Epetra_LinearProblem
// doesn't know about RCP, so we have to give it a raw pointer.
precProblem.SetOperator (K.getRawPtr ());
// Create AztecOO iterative solver for solving linear systems with K.
AztecOO precSolver (precProblem);
// Tell the solver to use the Ifpack preconditioner we created above.
precSolver.SetPrecOperator (Prec.get ());
// Set AztecOO solver options.
precSolver.SetAztecOption (AZ_output, AZ_none); // Don't print output
precSolver.SetAztecOption (AZ_solver, AZ_cg); // Use CG
// Use the above AztecOO solver to create the AztecOO_Operator.
// This is the place where we tell the AztecOO solver the maximum
// number of iterations (here, we use the matrix dimension; in
// practice, you'll want a smaller number) and the convergence
// tolerance (here, 1e-12).
RCP<AztecOO_Operator> precOperator =
rcp (new AztecOO_Operator (&precSolver, K->NumGlobalRows (), 1e-12));
// Create an Operator that computes y = K^{-1} M x.
//
// This Operator object is the operator we give to Anasazi. Thus,
// Anasazi just sees an operator that computes y = K^{-1} M x. The
// matrix K got absorbed into precOperator via precProblem (the
// Epetra_LinearProblem object). Later, when we set up the Anasazi
// eigensolver, we will need to tell it about M, so that it can
// orthogonalize basis vectors with respect to the inner product
// defined by M (since it is symmetric positive definite).
RCP<Anasazi::EpetraGenOp> Aop = rcp (new Anasazi::EpetraGenOp (precOperator, M));
//
// Set parameters for the block Krylov-Schur eigensolver
//
double tol = 1.0e-8;
int nev = 10;
int blockSize = 3;
int numBlocks = 3 * nev / blockSize;
int maxRestarts = 10;
// We're looking for the largest-magnitude eigenvalues of the
// _inverse_ operator, thus, the smallest-magnitude eigenvalues of
// the original operator.
std::string which = "LM";
// Create ParameterList to pass into eigensolver
Teuchos::ParameterList MyPL;
MyPL.set ("Verbosity", verbosity);
MyPL.set ("Which", which);
MyPL.set ("Block Size", blockSize);
MyPL.set ("Num Blocks", numBlocks);
MyPL.set ("Maximum Restarts", maxRestarts);
MyPL.set ("Convergence Tolerance", tol);
// Create an initial set of vectors to start the eigensolver. Note:
// This needs to have the same number of columns as the block size.
RCP<MV> ivec = rcp (new MV (K->Map (), blockSize));
MVT::MvRandom (*ivec);
// This object holds all the stuff about your problem that Anasazi
// will see.
//
// Anasazi only needs M so that it can orthogonalize basis vectors
// with respect to the M inner product. Wouldn't it be nice if
// Anasazi didn't require M in two different places? Alas, this is
// not currently the case.
RCP<Anasazi::BasicEigenproblem<double,MV,OP> > MyProblem =
rcp (new Anasazi::BasicEigenproblem<double,MV,OP> (Aop, M, ivec));
// Tell the eigenproblem that the matrix pencil (K,M) is symmetric.
MyProblem->setHermitian (true);
// Set the number of eigenvalues requested
MyProblem->setNEV (nev);
// Tell the eigenproblem that you are finished passing it information.
bool boolret = MyProblem->setProblem();
if (boolret != true) {
if (MyPID == 0) {
cout << "Anasazi::BasicEigenproblem::setProblem() returned with error." << endl;
}
#ifdef EPETRA_MPI
MPI_Finalize ();
#endif // EPETRA_MPI
return -1;
}
// Create the Block Krylov-Schur eigensolver.
Anasazi::BlockKrylovSchurSolMgr<double, MV, OP> MySolverMgr (MyProblem, MyPL);
// 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 = MySolverMgr.solve ();
if (returnCode != Anasazi::Converged && MyPID == 0) {
cout << "Anasazi::EigensolverMgr::solve() returned unconverged." << endl;
}
// Get the eigenvalues and eigenvectors from the eigenproblem.
Anasazi::Eigensolution<double,MV> sol = MyProblem->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;
Teuchos::RCP<MV> evecs = sol.Evecs;
int numev = sol.numVecs;
if (numev > 0) {
// Reconstruct the eigenvalues. The ones that Anasazi gave back
// are the inverses of the original eigenvalues. Reconstruct the
// eigenvectors too.
Teuchos::SerialDenseMatrix<int,double> dmatr(numev,numev);
MV tempvec (K->Map (), MVT::GetNumberVecs (*evecs));
K->Apply (*evecs, tempvec);
MVT::MvTransMv (1.0, tempvec, *evecs, dmatr);
if (MyPID == 0) {
double compeval = 0.0;
cout.setf (std::ios_base::right, std::ios_base::adjustfield);
cout << "Actual Eigenvalues (obtained by Rayleigh quotient) : " << endl;
cout << "------------------------------------------------------" << endl;
cout << std::setw(16) << "Real Part"
<< std::setw(16) << "Rayleigh Error" << endl;
cout << "------------------------------------------------------" << endl;
for (int i = 0; i < numev; ++i) {
compeval = dmatr(i,i);
cout << std::setw(16) << compeval
<< std::setw(16)
<< std::fabs (compeval - 1.0/evals[i].realpart)
<< endl;
}
cout << "------------------------------------------------------" << endl;
}
}
#ifdef EPETRA_MPI
MPI_Finalize ();
#endif // EPETRA_MPI
return 0;
}
Basic implementation of the Anasazi::Eigenproblem class.
The Anasazi::BlockKrylovSchurSolMgr class provides a user interface for the block Krylov-Schur eigens...
Declarations of Anasazi multi-vector and operator classes using Epetra_MultiVector and Epetra_Operato...
This provides a basic implementation for defining standard or generalized eigenvalue problems.
The Anasazi::BlockKrylovSchurSolMgr provides a flexible solver manager over the BlockKrylovSchur eige...
Adapter class for creating an operators often used in solving generalized eigenproblems.
Traits class which defines basic operations on multivectors.
ReturnType
Enumerated type used to pass back information from a solver manager.
Struct for storing an eigenproblem solution.
std::vector< Value< ScalarType > > Evals
The computed eigenvalues.