Anasazi  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
TraceMinDavidsonSpecTransEx.cpp

This is an example of how to use the TraceMinDavidsonSolMgr solver manager to compute the largest eigenpairs of a matrix via an invisible spectral transformation, using Tpetra data stuctures.

// @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 demonstrates how to use TraceMin-Davidson to find the
// largest eigenvalues of a matrix via spectral transformation
// Include autoconfigured header
// Include header for TraceMin-Davidson solver
// Include header to define basic eigenproblem Ax = \lambda*Bx
// Include header to provide Anasazi with Tpetra adapters
// Include header for Tpetra compressed-row storage matrix
#include "Tpetra_CrsMatrix.hpp"
#include "Tpetra_Core.hpp"
#include "Tpetra_Version.hpp"
#include "Tpetra_Map.hpp"
#include "Tpetra_MultiVector.hpp"
#include "Tpetra_Operator.hpp"
#include "Tpetra_Vector.hpp"
// Include headers for reading and writing matrix-market files
#include <MatrixMarket_Tpetra.hpp>
// Include header for sparse matrix operations
#include <TpetraExt_MatrixMatrix_def.hpp>
// Include header for Teuchos serial dense matrix
//#include "Teuchos_ArrayViewDecl.hpp"
int
main (int argc, char *argv[])
{
using Teuchos::RCP;
using std::cout;
using std::cin;
//
// Specify types used in this example
//
typedef Tpetra::CrsMatrix<>::scalar_type Scalar;
typedef Tpetra::CrsMatrix<Scalar> CrsMatrix;
typedef Tpetra::MultiVector<Scalar> MV;
typedef Tpetra::Operator<Scalar> OP;
//
// Initialize the MPI session
//
Tpetra::ScopeGuard tpetraScope(&argc,&argv);
//
// Get the default communicator
//
RCP<const Teuchos::Comm<int> > comm = Tpetra::getDefaultComm();
const int myRank = comm->getRank();
//
// Get parameters from command-line processor
//
std::string filenameA("/home/amklinv/matrices/bcsstk06.mtx");
Scalar tol = 1e-6;
int nev = 4;
int blockSize = 1;
bool verbose = true;
std::string whenToShift = "Always";
cmdp.setOption("fileA",&filenameA, "Filename for the Matrix-Market stiffness matrix.");
cmdp.setOption("tolerance",&tol, "Relative residual used for solver.");
cmdp.setOption("nev",&nev, "Number of desired eigenpairs.");
cmdp.setOption("blocksize",&blockSize, "Number of vectors to add to the subspace at each iteration.");
cmdp.setOption("verbose","quiet",&verbose, "Whether to print a lot of info or a little bit.");
cmdp.setOption("whenToShift",&whenToShift, "When to perform Ritz shifts. Options: Never, After Trace Levels, Always.");
return -1;
}
//
// Read the matrices from a file
//
RCP<const CrsMatrix> K = Tpetra::MatrixMarket::Reader<CrsMatrix>::readSparseFile(filenameA, comm);
//
// Compute the norm of the matrix
//
Scalar mat_norm = K->getFrobeniusNorm();
//
// ************************************
// Start the block Arnoldi iteration
// ************************************
//
// Variables used for the Block Arnoldi Method
//
int verbosity;
int numRestartBlocks = 2*nev/blockSize;
int numBlocks = 10*nev/blockSize;
if(verbose)
else
//
// Create parameter list to pass into solver
//
MyPL.set( "Verbosity", verbosity ); // How much information should the solver print?
MyPL.set( "Saddle Solver Type", "Projected Krylov"); // Use projected minres/gmres to solve the saddle point problem
MyPL.set( "Block Size", blockSize ); // Add blockSize vectors to the basis per iteration
MyPL.set( "Convergence Tolerance", tol*mat_norm ); // How small do the residuals have to be
MyPL.set( "Relative Convergence Tolerance", false); // Don't scale residuals by eigenvalues (when checking for convergence)
MyPL.set( "Use Locking", true); // Use deflation
MyPL.set( "Relative Locking Tolerance", false); // Don't scale residuals by eigenvalues (when checking whether to lock a vector)
MyPL.set("Num Restart Blocks", numRestartBlocks); // When we restart, we start back up with 2*nev blocks
MyPL.set("Num Blocks", numBlocks); // Maximum number of blocks in the subspace
MyPL.set("When To Shift", whenToShift);
MyPL.set("Which", "LM");
//
// 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.
//
RCP<MV> ivec = Teuchos::rcp( new MV(K->getRowMap(), numRestartBlocks*blockSize) );
MVT::MvRandom( *ivec );
//
// Create the eigenproblem
//
RCP<Anasazi::BasicEigenproblem<Scalar,MV,OP> > MyProblem =
//
// 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 != true) {
if (myRank == 0) {
cout << "Anasazi::BasicEigenproblem::setProblem() returned with error." << std::endl;
}
return -1;
}
//
// Initialize the TraceMin-Davidson solver
//
//
// Solve the problem to the specified tolerances
//
Anasazi::ReturnType returnCode = MySolverMgr.solve();
if (returnCode != Anasazi::Converged && myRank == 0) {
cout << "Anasazi::EigensolverMgr::solve() returned unconverged." << std::endl;
}
else if (myRank == 0)
cout << "Anasazi::EigensolverMgr::solve() returned converged." << std::endl;
//
// Get the eigenvalues and eigenvectors from the eigenproblem
//
Anasazi::Eigensolution<Scalar,MV> sol = MyProblem->getSolution();
std::vector<Anasazi::Value<Scalar> > evals = sol.Evals;
RCP<MV> evecs = sol.Evecs;
int numev = sol.numVecs;
//
// Compute the residual, just as a precaution
//
if (numev > 0) {
for(int i=0; i < numev; i++)
T(i,i) = evals[i].realpart;
std::vector<Scalar> normR(sol.numVecs);
MV Kvec( K->getRowMap(), MVT::GetNumberVecs( *evecs ) );
OPT::Apply( *K, *evecs, Kvec );
MVT::MvTimesMatAddMv( -1.0, *evecs, T, 1.0, Kvec );
MVT::MvNorm( Kvec, normR );
if (myRank == 0) {
cout.setf(std::ios_base::right, std::ios_base::adjustfield);
cout<<"Actual Eigenvalues (obtained by Rayleigh quotient) : "<<std::endl;
cout<<"------------------------------------------------------"<<std::endl;
cout<<std::setw(16)<<"Real Part"
<<std::setw(16)<<"Error"<<std::endl;
cout<<"------------------------------------------------------"<<std::endl;
for (int i=0; i<numev; i++) {
cout<<std::setw(16)<<evals[i].realpart
<<std::setw(16)<<normR[i]/mat_norm
<<std::endl;
}
cout<<"------------------------------------------------------"<<std::endl;
}
}
return 0;
}