doc/html/BlockKrylovSchur_2BlockKrylovSchurEpetraExGenAmesos_8cpp-example.html

// @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 the KLU sparse direct linear solver, and

// accesses KLU through Trilinos' Amesos package.

//

// 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 the sparse direct solver

// KLU to do so.  Trilinos' Amesos package has an interface to KLU.


// Include header for block Krylov-Schur solver

#include "AnasaziBlockKrylovSchurSolMgr.hpp"

// Include header to define basic eigenproblem Ax = \lambda*Bx

#include "AnasaziBasicEigenproblem.hpp"

// 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 "AnasaziEpetraAdapter.hpp"

// Include header for Epetra sparse matrix, Map (representation of

// parallel distributions), and linear problem.  Amesos uses the

// latter to encapsulate linear systems to solve.

#include "Epetra_Map.h"

#include "Epetra_CrsMatrix.h"

#include "Epetra_LinearProblem.h"

// Include header for Amesos solver.

#include "Amesos.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 required by Epetra objects

#ifdef EPETRA_MPI

#  include "Epetra_MpiComm.h"

#else

#  include "Epetra_SerialComm.h"

#endif


// \class AmesosGenOp

// \brief Operator that computes \f$Y = K^{-1} M X\f$ or \f$Y = M

//   K^{-T} X\f$, where solves with the Epetra_CrsMatrix K use an

//   Amesos solver, and M is an Epetra_Operator.

// \author Heidi K. Thornquist

// \date Jan 9, 2006

//

// This operator calls an inverse method in Apply().  It may use a

// Belos solve, Amesos solve, or AztecOO solve.  This is used to

// implement shift and invert, with a shift to zero.  Shift before you

// invert (i.e., solve).  This operator goes into Block Krylov-Schur.

//

// Anasazi will ask for the largest eigenvalues of the inverse, thus,

// the smallest eigenvalues of the original (shifted) matrix.  If you

// wanted a shift, just apply the shift first, and send that to this

// operator.

class AmesosGenOp : public virtual Epetra_Operator {

public:

  // \brief Constructor

  //

  // \param solver [in] The Amesos solver to use for solving linear

  //   systems with K.  It must already have its Epetra_LinearProblem

  //   set, and that Epetra_LinearProblem must have K.

  // \param massMtx [in] The "mass matrix" M.

  // \param useTranspose [in] If false, apply \f$K^{-1} M\f$; else,

  //   apply \f$M K^{-T}\f$.

  AmesosGenOp (const Teuchos::RCP<Amesos_BaseSolver>& solver,

               const Teuchos::RCP<Epetra_Operator>& massMtx,

               const bool useTranspose = false );

  // Destructor

  virtual ~AmesosGenOp () {}


  // \brief Compute Y such that \f$K^{-1} M X = Y\f$ or \f$M K^{-T} X

  //   = Y\f$, where solves with the Epetra_CrsMatrix K use an Amesos

  //   solver.

  int Apply (const Epetra_MultiVector& X, Epetra_MultiVector& Y) const;


  // The Operator's (human-readable) label.

  const char* Label() const {

    return "Operator that applies K^{-1} M or M K^{-T}";

  }


  // Whether to apply \f$K^{-1} M\f$ (false) or \f$K^{-T} M\f$ (true).

  bool UseTranspose() const { return useTranspose_; };


  // \brief Set whether to apply \f$K^{-1} M\f$ (false) or

  //   \f$K^{-T} M\f$ (true).

  //

  // \param useTranspose [in] If true, tell this Operator to apply

  //   \f$K^{-T} M\f$, from now until the next call to

  //   SetUseTranspose().  If false, tell this Operator to apply

  //   \f$K^{-1} M\f$, until the next call to SetUseTranspose().

  int SetUseTranspose (bool useTranspose);


  // The Operator's communicator.

  const Epetra_Comm& Comm () const {

    return solver_->Comm ();

  }


  // The Operator's domain Map.

  const Epetra_Map& OperatorDomainMap () const {

    return massMtx_->OperatorDomainMap ();

  }


  // The Operator's range Map.

  const Epetra_Map& OperatorRangeMap () const {

    return massMtx_->OperatorRangeMap ();

  }


  // \brief NOT IMPLEMENTED: Apply \f$( K^{-1} M )^{-1}\f$ or

  //   \f$( K^{-1} M )^{-T}\f$.

  //

  // Returning a nonzero value means that this Operator doesn't know

  // how to apply its inverse.

  int ApplyInverse (const Epetra_MultiVector& X, Epetra_MultiVector& Y) const {

    return -1;

  };


  // NOT IMPLEMENTED: Whether this Operator can compute its infinity norm.

  bool HasNormInf() const { return false; }


  // \brief NOT IMPLEMENTED: Infinity norm of \f$( K^{-1} M )^{-1}\f$

  //   or \f$( K^{-1} M )^{-T}\f$.

  //

  // Returning -1.0 means that the Operator does not know how to

  // compute its infinity norm.

  double NormInf () const { return -1.0; }


private:

  // Default constructor: You may not call this.

  AmesosGenOp () {};


  // Copy constructor: You may not call this.

  AmesosGenOp (const AmesosGenOp& genOp);


  Teuchos::RCP<Amesos_BaseSolver> solver_;

  Teuchos::RCP<Epetra_Operator> massMtx_;

  Epetra_LinearProblem* problem_;

  bool useTranspose_;

};


// ****************************************************************************

// 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;

  typedef Anasazi::MultiVecTraits<double, MV> MVT;


#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

  // Amesos direct solver.

  //


  // Create Epetra linear problem class for solving linear systems

  // with K.  This implements the inverse operator for shift and

  // invert.

  Epetra_LinearProblem AmesosProblem;

  // Tell the linear problem about the matrix K.  Epetra_LinearProblem

  // doesn't know about RCP, so we have to give it a raw pointer.

  AmesosProblem.SetOperator (K.get ());


  // Create Amesos factory and solver for solving linear systems with

  // K.  The solver uses the KLU library to do a sparse LU

  // factorization.

  //

  // Note that the AmesosProblem object "absorbs" K.  Anasazi doesn't

  // see K, just the operator that implements K^{-1} M.

  Amesos amesosFactory;

  RCP<Amesos_BaseSolver> AmesosSolver;


  // Amesos can interface to many different solvers.  The following

  // loop picks a solver that Amesos supports.  The loop order

  // reflects solver preference, only in the sense that using LAPACK

  // here is a suboptimal fall-back.  (With the LAPACK option, Amesos

  // makes a dense version of the sparse matrix and uses LAPACK to

  // compute the factorization.  The other options are true sparse

  // direct factorizations.)

  const int numSolverNames = 9;

  const char* solverNames[9] = {

    "Klu", "Umfpack", "Superlu", "Superludist", "Mumps",

    "Paradiso", "Taucs", "CSparse", "Lapack"

  };

  for (int k = 0; k < numSolverNames; ++k) {

    const char* const solverName = solverNames[k];

    if (amesosFactory.Query (solverName)) {

      AmesosSolver = rcp (amesosFactory.Create (solverName, AmesosProblem));

      if (MyPID == 0) {

        cout << "Amesos solver: \"" << solverName << "\"" << endl;

      }

    }

  }

  if (AmesosSolver.is_null ()) {

    throw std::runtime_error ("Amesos appears not to have any solvers enabled.");

  }


  // The AmesosGenOp class assumes that the symbolic and numeric

  // factorizations have already been performed on the linear problem.

  AmesosSolver->SymbolicFactorization ();

  AmesosSolver->NumericFactorization ();


  //

  // Set parameters for the block Krylov-Schur eigensolver

  //


  double tol = 1.0e-8; // convergence tolerance

  int nev = 10; // number of eigenvalues for which to solve

  int blockSize = 3; // block size (number of eigenvectors processed at once)

  int numBlocks = 3 * nev / blockSize; // restart length

  int maxRestarts = 5; // maximum number of restart cycles


  // 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";

  int verbosity = Anasazi::Errors + Anasazi::Warnings + Anasazi::FinalSummary;


  // 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));

  ivec->Random ();


  // Create the Epetra_Operator for the spectral transformation using

  // the Amesos direct solver.

  //

  // The AmesosGenOp 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 AmesosProblem (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<AmesosGenOp> Aop = rcp (new AmesosGenOp (AmesosSolver, M));


  // Create the eigenproblem.  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.

  const bool boolret = MyProblem->setProblem ();

  if (boolret != true) {

    if (MyPID == 0) {

      cerr << "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 eigensolver did not converge." << 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;

  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.

    MV tempvec (K->Map (), MVT::GetNumberVecs (*evecs));

    K->Apply (*evecs, tempvec);

    Teuchos::SerialDenseMatrix<int,double> dmatr (numev, numev);

    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;

}


AmesosGenOp::

AmesosGenOp (const Teuchos::RCP<Amesos_BaseSolver>& solver,

             const Teuchos::RCP<Epetra_Operator>& massMtx,

             const bool useTranspose)

  : solver_ (solver),

    massMtx_ (massMtx),

    problem_ (NULL),

    useTranspose_ (useTranspose)

{

  if (solver.is_null ()) {

    throw std::invalid_argument ("AmesosGenOp constructor: The 'solver' "

                                 "input argument is null.");

  }

  if (massMtx.is_null ()) {

    throw std::invalid_argument ("AmesosGenOp constructor: The 'massMtx' "

                                 "input argument is null.");

  }


  Epetra_LinearProblem* problem = const_cast<Epetra_LinearProblem*> (solver->GetProblem ());

  if (problem == NULL) {

    throw std::invalid_argument ("The solver's GetProblem() method returned "

                                 "NULL.  This probably means that its "

                                 "LinearProblem has not yet been set.");

  }

  problem_ = problem;


  if (solver_->UseTranspose ()) {

    solver_->SetUseTranspose (! useTranspose);

  } else {

    solver_->SetUseTranspose (useTranspose);

  }


  if (massMtx_->UseTranspose ()) {

    massMtx_->SetUseTranspose (! useTranspose);

  } else {

    massMtx_->SetUseTranspose (useTranspose);

  }

}


int

AmesosGenOp::SetUseTranspose (bool useTranspose)

{

  int err = 0;


  if (problem_ == NULL) {

    throw std::logic_error ("AmesosGenOp::SetUseTranspose: problem_ is NULL");

  }

  if (massMtx_.is_null ()) {

    throw std::logic_error ("AmesosGenOp::SetUseTranspose: massMtx_ is null");

  }

  if (solver_.is_null ()) {

    throw std::logic_error ("AmesosGenOp::SetUseTranspose: solver_ is null");

  }


  const bool solverUsesTranspose = solver_->UseTranspose ();

  if (solverUsesTranspose) {

    err = solver_->SetUseTranspose (! useTranspose);

  } else {

    err = solver_->SetUseTranspose (useTranspose);

  }


  // If SetUseTranspose returned zero above, then the Amesos solver

  // doesn't know how to change the transpose state.

  if (err != 0) {

    return err;

  }


  if (massMtx_->UseTranspose ()) {

    err = massMtx_->SetUseTranspose (! useTranspose);

  } else {

    err = massMtx_->SetUseTranspose (useTranspose);

  }


  // If SetUseTranspose returned zero above, then the mass matrix

  // doesn't know how to change the transpose state.

  if (err != 0) {

    // Put the solver back like we found it.

    (void) solver_->SetUseTranspose (solverUsesTranspose);

    return err;

  }


  useTranspose_ = useTranspose;

  return 0; // the method completed correctly

}


int

AmesosGenOp::Apply (const Epetra_MultiVector& X, Epetra_MultiVector& Y) const

{

  if (problem_ == NULL) {

    throw std::logic_error ("AmesosGenOp::Apply: problem_ is NULL");

  }

  if (massMtx_.is_null ()) {

    throw std::logic_error ("AmesosGenOp::Apply: massMtx_ is null");

  }

  if (solver_.is_null ()) {

    throw std::logic_error ("AmesosGenOp::Apply: solver_ is null");

  }


  if (! useTranspose_) {

    // Storage for M*X

    Epetra_MultiVector MX (X.Map (), X.NumVectors ());


    // Apply M*X

    massMtx_->Apply (X, MX);

    Y.PutScalar (0.0);


    // Set the LHS and RHS

    problem_->SetRHS (&MX);

    problem_->SetLHS (&Y);


    // Solve the linear system A*Y = MX

    solver_->Solve ();

  }

  else { // apply the transposed operator

    // Storage for A^{-T}*X

    Epetra_MultiVector ATX (X.Map (), X.NumVectors ());

    Epetra_MultiVector tmpX = const_cast<Epetra_MultiVector&> (X);


    // Set the LHS and RHS

    problem_->SetRHS (&tmpX);

    problem_->SetLHS (&ATX);


    // Solve the linear system A^T*Y = X

    solver_->Solve ();


    // Apply M*ATX

    massMtx_->Apply (ATX, Y);

  }


  return 0; // the method completed correctly

}