Anasazi Version of the Day
Loading...
Searching...
No Matches
LOBPCGEpetraExGenShifted.cpp

Use LOBPCG with Epetra, with shifted eigenvalue problem.

Use LOBPCG with Epetra, with shifted eigenvalue problem

This example computes the eigenvalues of largest magnitude of the discretized 2-D Laplacian operator, using Anasazi's implementation of the LOBPCG method. This problem constructs a shifted eigenproblem that targets the smallest eigenvalues around a certain value (sigma). This operator is discretized using linear finite elements and constructed as an Epetra matrix, then passed shifted using EpetraExt utilities.

// @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
// Include autoconfigured header
// Include header for LOBPCG solver
// Include header to define basic eigenproblem Ax = \lambda*Bx
// Include header to provide Anasazi with Epetra adapters
// Include header to provide basic Anasazi output manager
// Include header for Epetra compressed-row storage matrix and linear problem
#include "Epetra_CrsMatrix.h"
#include "Epetra_LinearProblem.h"
// Include header for Teuchos serial dense matrix
#include "Teuchos_SerialDenseMatrix.hpp"
// Include header for the problem definition
#include "ModeLaplace2DQ2.h"
// Include EpetraExt MatrixMatrix helpers.
#include "EpetraExt_MatrixMatrix.h"
// Include selected communicator class and map required by Epetra objects
#ifdef EPETRA_MPI
#include "Epetra_MpiComm.h"
#else
#include "Epetra_SerialComm.h"
#endif
#include "Epetra_Map.h"
#include "Teuchos_StandardCatchMacros.hpp"
int
main (int argc, char *argv[])
{
#ifdef EPETRA_MPI
// Initialize MPI
MPI_Init(&argc,&argv);
Epetra_MpiComm Comm(MPI_COMM_WORLD);
#else
Epetra_SerialComm Comm;
#endif
bool success = false;
try {
int i;
int MyPID = Comm.MyPID();
// Create an Anasazi output manager
//
printer.stream(Anasazi::Errors) << Anasazi::Anasazi_Version() << std::endl << std::endl;
// Number of dimension of the domain
int space_dim = 2;
// Size of each of the dimensions of the 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 problem
Teuchos::RCP<ModalProblem> testCase = Teuchos::rcp( new ModeLaplace2DQ2(Comm, brick_dim[0], elements[0], brick_dim[1], elements[1]) );
// Get the stiffness and mass matrices
Teuchos::RCP<Epetra_CrsMatrix> K = Teuchos::rcp( const_cast<Epetra_CrsMatrix *>(testCase->getStiffness()), false );
Teuchos::RCP<Epetra_CrsMatrix> M = Teuchos::rcp( const_cast<Epetra_CrsMatrix *>(testCase->getMass()), false );
// Create the shifted system K - sigma * M.
double sigma = 1.0;
Teuchos::RCP<Epetra_CrsMatrix> Kshift = Teuchos::rcp( new Epetra_CrsMatrix( *K ) );
int addErr = EpetraExt::MatrixMatrix::Add( *M, false, -sigma, *Kshift, 1.0 );
if (addErr != 0) {
printer.print(Anasazi::Errors,"EpetraExt::MatrixMatrix::Add returned with error.\n");
throw -1;
}
//
// ************************************
// Start the LOBPCG iteration
// ************************************
//
// Variables used for the LOBPCG Method
//
const int nev = 10;
const int blockSize = 5;
const int maxIters = 500;
const double tol = 1.0e-8;
std::string which = "SM";
//
// Create parameter list to pass into solver
//
Teuchos::ParameterList MyPL;
MyPL.set( "Verbosity", verbosity );
MyPL.set( "Which", which );
MyPL.set( "Block Size", blockSize );
MyPL.set( "Maximum Iterations", maxIters );
MyPL.set( "Convergence Tolerance", tol );
MyPL.set( "Full Ortho", true );
MyPL.set( "Use Locking", true );
MyPL.set( "Locking Tolerance", tol/10 );
typedef Epetra_MultiVector MV;
typedef Epetra_Operator OP;
// typedef Anasazi::OperatorTraits<double, MV, OP> OPT;
// Create an Epetra_MultiVector for an initial vector to start the solver.
// Note: This needs to have the same number of columns as the blocksize.
Teuchos::RCP<Epetra_MultiVector> ivec = Teuchos::rcp( new Epetra_MultiVector(K->Map(), blockSize) );
MVT::MvRandom( *ivec );
Teuchos::RCP<Anasazi::BasicEigenproblem<double,MV,OP> > MyProblem =
Teuchos::rcp( new Anasazi::BasicEigenproblem<double,MV,OP>(Kshift, M, ivec) );
// Inform the eigenproblem that the matrix pencil (K,M) is symmetric
MyProblem->setHermitian(true);
// Set the number of eigenvalues requested
MyProblem->setNEV( nev );
// Inform the eigenproblem that you are finished passing it information
bool boolret = MyProblem->setProblem();
if (!boolret) {
printer.print(Anasazi::Errors,"Anasazi::BasicEigenproblem::setProblem() returned an error.\n");
throw -1;
}
// Initialize the LOBPCG solver
Anasazi::LOBPCGSolMgr<double, MV, OP> MySolverMgr(MyProblem, MyPL);
// Solve the problem to the specified tolerances or length
Anasazi::ReturnType returnCode = MySolverMgr.solve();
if (returnCode != Anasazi::Converged && MyPID==0) {
std::cout << "Anasazi::EigensolverMgr::solve() returned unconverged." << std::endl;
}
// Get the eigenvalues and eigenvectors from the eigenproblem
Anasazi::Eigensolution<double,MV> sol = MyProblem->getSolution();
std::vector<Anasazi::Value<double> > evals = sol.Evals;
Teuchos::RCP<MV> evecs = sol.Evecs;
int numev = sol.numVecs;
if (numev > 0) {
// Undo shift transformation; computed eigenvalues are real
std::vector<double> compEvals(numev);
for (i=0; i<numev; ++i) {
compEvals[i] = evals[i].realpart + sigma;
}
//************************************
// Compute residuals, just for funsies
//************************************
//
std::vector<double> normR(sol.numVecs);
Teuchos::SerialDenseMatrix<int,double> T(sol.numVecs, sol.numVecs);
Epetra_MultiVector Kvec( K->OperatorDomainMap(), evecs->NumVectors() );
Epetra_MultiVector Mvec( M->OperatorDomainMap(), evecs->NumVectors() );
T.putScalar(0.0);
for (i=0; i<sol.numVecs; i++) {
T(i,i) = compEvals[i];
}
K->Apply( *evecs, Kvec );
M->Apply( *evecs, Mvec );
MVT::MvTimesMatAddMv( -1.0, Mvec, T, 1.0, Kvec );
MVT::MvNorm( Kvec, normR );
//************************************
// Print the results
//************************************
//
std::ostringstream os;
os.setf(std::ios_base::right, std::ios_base::adjustfield);
os<<"Solver manager returned " << (returnCode == Anasazi::Converged ? "converged." : "unconverged.") << std::endl;
os<<std::endl;
os<<"------------------------------------------------------"<<std::endl;
os<<std::setw(16)<<"Eigenvalue"
<<std::setw(18)<<"Direct Residual"
<<std::endl;
os<<"------------------------------------------------------"<<std::endl;
for (i=0; i<sol.numVecs; i++) {
os<<std::setw(16)<<compEvals[i]
<<std::setw(18)<<normR[i]/compEvals[i]
<<std::endl;
}
os<<"------------------------------------------------------"<<std::endl;
printer.print(Anasazi::Errors,os.str());
}
success = true;
}
TEUCHOS_STANDARD_CATCH_STATEMENTS(true, std::cerr, success);
#ifdef EPETRA_MPI
MPI_Finalize();
#endif
return ( success ? EXIT_SUCCESS : EXIT_FAILURE );
}
Basic implementation of the Anasazi::Eigenproblem class.
Basic output manager for sending information of select verbosity levels to the appropriate output str...
Anasazi header file which uses auto-configuration information to include necessary C++ headers.
Declarations of Anasazi multi-vector and operator classes using Epetra_MultiVector and Epetra_Operato...
The Anasazi::LOBPCGSolMgr provides a powerful solver manager for the LOBPCG eigensolver.
This provides a basic implementation for defining standard or generalized eigenvalue problems.
Anasazi's basic output manager for sending information of select verbosity levels to the appropriate ...
User interface for the LOBPCG eigensolver.
ReturnType solve()
This method performs possibly repeated calls to the underlying eigensolver's iterate() routine until ...
Traits class which defines basic operations on multivectors.
virtual Teuchos::FancyOStream & stream(MsgType type)
Create a stream for outputting to.
virtual void print(MsgType type, const std::string output)
Send output to the output manager.
ReturnType
Enumerated type used to pass back information from a solver manager.
Struct for storing an eigenproblem solution.
Teuchos::RCP< MV > Evecs
The computed eigenvectors.
int numVecs
The number of computed eigenpairs.
std::vector< Value< ScalarType > > Evals
The computed eigenvalues.