Anasazi Version of the Day
Loading...
Searching...
No Matches
BlockDavidson/BlockDavidsonEpetraEx.cpp

This is an example of how to use the Anasazi::BlockDavidsonSolMgr solver manager to solve a standard eigenvalue problem, using Epetra data structures.

// @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 "Epetra_CrsMatrix.h"
#include "Teuchos_CommandLineProcessor.hpp"
#include "Teuchos_StandardCatchMacros.hpp"
#include "Teuchos_Assert.hpp"
#ifdef HAVE_MPI
#include "Epetra_MpiComm.h"
#include <mpi.h>
#else
#include "Epetra_SerialComm.h"
#endif
#include "Epetra_Map.h"
int main(int argc, char *argv[]) {
#ifdef HAVE_MPI
// Initialize MPI
//
MPI_Init(&argc,&argv);
#endif
bool success = false;
try {
// Create an Epetra communicator
//
#ifdef HAVE_MPI
Epetra_MpiComm Comm(MPI_COMM_WORLD);
#else
Epetra_SerialComm Comm;
#endif
// Create an Anasazi output manager
//
printer.stream(Anasazi::Errors) << Anasazi::Anasazi_Version() << std::endl << std::endl;
// Get the sorting string from the command line
//
std::string which("SM");
Teuchos::CommandLineProcessor cmdp(false,true);
cmdp.setOption("sort",&which,"Targetted eigenvalues (SM or LM).");
if (cmdp.parse(argc,argv) != Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL) {
throw -1;
}
// Dimension of the matrix
//
// Discretization points in any one direction.
//
const int nx = 10;
// Size of matrix nx*nx
//
const int NumGlobalElements = nx*nx;
// Construct a Map that puts approximately the same number of
// equations on each processor.
//
Epetra_Map Map(NumGlobalElements, 0, Comm);
// Get update list and number of local equations from newly created Map.
//
int NumMyElements = Map.NumMyElements();
std::vector<int> MyGlobalElements(NumMyElements);
Map.MyGlobalElements(&MyGlobalElements[0]);
// Create an integer vector NumNz that is used to build the Petra Matrix.
// NumNz[i] is the Number of OFF-DIAGONAL term for the ith global equation
// on this processor
//
std::vector<int> NumNz(NumMyElements);
/* We are building a matrix of block structure:
| T -I |
|-I T -I |
| -I T |
| ... -I|
| -I T|
where each block is dimension nx by nx and the matrix is on the order of
nx*nx. The block T is a tridiagonal matrix.
*/
for (int i=0; i<NumMyElements; i++) {
if (MyGlobalElements[i] == 0 || MyGlobalElements[i] == NumGlobalElements-1 ||
MyGlobalElements[i] == nx-1 || MyGlobalElements[i] == nx*(nx-1) ) {
NumNz[i] = 3;
}
else if (MyGlobalElements[i] < nx || MyGlobalElements[i] > nx*(nx-1) ||
MyGlobalElements[i]%nx == 0 || (MyGlobalElements[i]+1)%nx == 0) {
NumNz[i] = 4;
}
else {
NumNz[i] = 5;
}
}
// Create an Epetra_Matrix
//
Teuchos::RCP<Epetra_CrsMatrix> A = Teuchos::rcp( new Epetra_CrsMatrix(Epetra_DataAccess::Copy, Map, &NumNz[0]) );
// Compute coefficients for discrete convection-diffution operator
//
const double one = 1.0;
std::vector<double> Values(4);
std::vector<int> Indices(4);
double rho = 0.0;
double h = one /(nx+1);
double h2 = h*h;
double c = 5.0e-01*rho/ h;
Values[0] = -one/h2 - c; Values[1] = -one/h2 + c; Values[2] = -one/h2; Values[3]= -one/h2;
double diag = 4.0 / h2;
int NumEntries;
for (int i=0; i<NumMyElements; i++)
{
if (MyGlobalElements[i]==0)
{
Indices[0] = 1;
Indices[1] = nx;
NumEntries = 2;
int info = A->InsertGlobalValues(MyGlobalElements[i], NumEntries, &Values[1], &Indices[0]);
TEUCHOS_TEST_FOR_EXCEPTION( info != 0, std::runtime_error, "Failure in InsertGlobalValues()" );
}
else if (MyGlobalElements[i] == nx*(nx-1))
{
Indices[0] = nx*(nx-1)+1;
Indices[1] = nx*(nx-2);
NumEntries = 2;
int info = A->InsertGlobalValues(MyGlobalElements[i], NumEntries, &Values[1], &Indices[0]);
TEUCHOS_TEST_FOR_EXCEPTION( info != 0, std::runtime_error, "Failure in InsertGlobalValues()" );
}
else if (MyGlobalElements[i] == nx-1)
{
Indices[0] = nx-2;
NumEntries = 1;
int info = A->InsertGlobalValues(MyGlobalElements[i], NumEntries, &Values[0], &Indices[0]);
assert( info==0 );
Indices[0] = 2*nx-1;
info = A->InsertGlobalValues(MyGlobalElements[i], NumEntries, &Values[2], &Indices[0]);
TEUCHOS_TEST_FOR_EXCEPTION( info != 0, std::runtime_error, "Failure in InsertGlobalValues()" );
}
else if (MyGlobalElements[i] == NumGlobalElements-1)
{
Indices[0] = NumGlobalElements-2;
NumEntries = 1;
int info = A->InsertGlobalValues(MyGlobalElements[i], NumEntries, &Values[0], &Indices[0]);
assert( info==0 );
Indices[0] = nx*(nx-1)-1;
info = A->InsertGlobalValues(MyGlobalElements[i], NumEntries, &Values[2], &Indices[0]);
TEUCHOS_TEST_FOR_EXCEPTION( info != 0, std::runtime_error, "Failure in InsertGlobalValues()" );
}
else if (MyGlobalElements[i] < nx)
{
Indices[0] = MyGlobalElements[i]-1;
Indices[1] = MyGlobalElements[i]+1;
Indices[2] = MyGlobalElements[i]+nx;
NumEntries = 3;
int info = A->InsertGlobalValues(MyGlobalElements[i], NumEntries, &Values[0], &Indices[0]);
TEUCHOS_TEST_FOR_EXCEPTION( info != 0, std::runtime_error, "Failure in InsertGlobalValues()" );
}
else if (MyGlobalElements[i] > nx*(nx-1))
{
Indices[0] = MyGlobalElements[i]-1;
Indices[1] = MyGlobalElements[i]+1;
Indices[2] = MyGlobalElements[i]-nx;
NumEntries = 3;
int info = A->InsertGlobalValues(MyGlobalElements[i], NumEntries, &Values[0], &Indices[0]);
TEUCHOS_TEST_FOR_EXCEPTION( info != 0, std::runtime_error, "Failure in InsertGlobalValues()" );
}
else if (MyGlobalElements[i]%nx == 0)
{
Indices[0] = MyGlobalElements[i]+1;
Indices[1] = MyGlobalElements[i]-nx;
Indices[2] = MyGlobalElements[i]+nx;
NumEntries = 3;
int info = A->InsertGlobalValues(MyGlobalElements[i], NumEntries, &Values[1], &Indices[0]);
TEUCHOS_TEST_FOR_EXCEPTION( info != 0, std::runtime_error, "Failure in InsertGlobalValues()" );
}
else if ((MyGlobalElements[i]+1)%nx == 0)
{
Indices[0] = MyGlobalElements[i]-nx;
Indices[1] = MyGlobalElements[i]+nx;
NumEntries = 2;
int info = A->InsertGlobalValues(MyGlobalElements[i], NumEntries, &Values[2], &Indices[0]);
assert( info==0 );
Indices[0] = MyGlobalElements[i]-1;
NumEntries = 1;
info = A->InsertGlobalValues(MyGlobalElements[i], NumEntries, &Values[0], &Indices[0]);
TEUCHOS_TEST_FOR_EXCEPTION( info != 0, std::runtime_error, "Failure in InsertGlobalValues()" );
}
else
{
Indices[0] = MyGlobalElements[i]-1;
Indices[1] = MyGlobalElements[i]+1;
Indices[2] = MyGlobalElements[i]-nx;
Indices[3] = MyGlobalElements[i]+nx;
NumEntries = 4;
int info = A->InsertGlobalValues(MyGlobalElements[i], NumEntries, &Values[0], &Indices[0]);
TEUCHOS_TEST_FOR_EXCEPTION( info != 0, std::runtime_error, "Failure in InsertGlobalValues()" );
}
// Put in the diagonal entry
int info = A->InsertGlobalValues(MyGlobalElements[i], 1, &diag, &MyGlobalElements[i]);
TEUCHOS_TEST_FOR_EXCEPTION( info != 0, std::runtime_error, "Failure in InsertGlobalValues()" );
}
// Finish up
int info = A->FillComplete();
TEUCHOS_ASSERT( info==0 );
A->SetTracebackMode(1); // Shutdown Epetra Warning tracebacks
// Create a identity matrix for the temporary mass matrix
Teuchos::RCP<Epetra_CrsMatrix> M = Teuchos::rcp( new Epetra_CrsMatrix(Epetra_DataAccess::Copy, Map, 1) );
for (int i=0; i<NumMyElements; i++)
{
Values[0] = one;
Indices[0] = i;
NumEntries = 1;
info = M->InsertGlobalValues(MyGlobalElements[i], NumEntries, &Values[0], &Indices[0]);
TEUCHOS_ASSERT( info==0 );
}
// Finish up
info = M->FillComplete();
TEUCHOS_ASSERT( info==0 );
M->SetTracebackMode(1); // Shutdown Epetra Warning tracebacks
//************************************
// Call the Block Davidson solver manager
//***********************************
//
// Variables used for the Block Davidson Method
//
const int nev = 4;
const int blockSize = 5;
const int numBlocks = 8;
const int maxRestarts = 100;
const double tol = 1.0e-8;
typedef Epetra_MultiVector MV;
typedef Epetra_Operator OP;
// 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(Map, blockSize) );
ivec->Random();
// Create the eigenproblem.
//
Teuchos::RCP<Anasazi::BasicEigenproblem<double, MV, OP> > MyProblem =
Teuchos::rcp( new Anasazi::BasicEigenproblem<double, MV, OP>(A, ivec) );
// Inform the eigenproblem that the operator A is symmetric
//
MyProblem->setHermitian(true);
// Set the number of eigenvalues requested
//
MyProblem->setNEV( nev );
// Inform the eigenproblem that you are finishing passing it information
//
bool boolret = MyProblem->setProblem();
if (boolret != true) {
printer.print(Anasazi::Errors,"Anasazi::BasicEigenproblem::setProblem() returned an error.\n");
throw -1;
}
// Create parameter list to pass into the solver manager
//
Teuchos::ParameterList MyPL;
MyPL.set( "Which", which );
MyPL.set( "Block Size", blockSize );
MyPL.set( "Num Blocks", numBlocks );
MyPL.set( "Maximum Restarts", maxRestarts );
MyPL.set( "Convergence Tolerance", tol );
MyPL.set( "Verbosity", verbosity );
//
// Create the solver manager
Anasazi::BlockDavidsonSolMgr<double, MV, OP> MySolverMan(MyProblem, MyPL);
// Solve the problem
//
Anasazi::ReturnType returnCode = MySolverMan.solve();
// 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;
// Compute residuals.
//
std::vector<double> normR(sol.numVecs);
if (sol.numVecs > 0) {
Teuchos::SerialDenseMatrix<int,double> T(sol.numVecs, sol.numVecs);
Epetra_MultiVector 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
//
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 (int i=0; i<sol.numVecs; i++) {
os<<std::setw(16)<<evals[i].realpart
<<std::setw(18)<<normR[i]/evals[i].realpart
<<std::endl;
}
os<<"------------------------------------------------------"<<std::endl;
printer.print(Anasazi::Errors,os.str());
success = true;
}
TEUCHOS_STANDARD_CATCH_STATEMENTS(true, std::cerr, success);
#ifdef HAVE_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...
The Anasazi::BlockDavidsonSolMgr provides a solver manager for the BlockDavidson eigensolver.
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...
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 ...
The BlockDavidsonSolMgr provides a powerful solver manager over the BlockDavidson 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.