AnasaziBlockKrylovSchurSVDWithCustomATAWrapper

// This example shows how to use the block Krylov-Schur method to
// compute a few of the largest singular values (sigma) and
// corresponding right singular vectors (v) for the matrix A by
// solving the symmetric problem:
//
//                  (A^T * (A * v) = sigma * v
//
// where A is an m by n real matrix that is derived from the simplest
// finite difference discretization of the 2-dimensional kernel
// K(s,t)dt, where
//
//                  K(s,t) = s(t-1)   if 0 <= s <= t <= 1
//                           t(s-1)   if 0 <= t <= s <= 1
//
// NOTE:  This example came from the ARPACK SVD driver dsvd.f
//
// The main solver parameters are the number of singular values to
// compute (numSingularValues), and the number of starting vectors
// (blockSize).  The implementation of getParameterList() below
// includes other parameters, like the iteration tolerance and the
// maximum number of restart cycles.
//
// This example differs from AnasaziBlockKrylovSchurSVD because it
// uses a custom EpetraSymOp wrapper, rather than
// Anasazi::EpetraSymOp.
//
#include "AnasaziBlockKrylovSchurSolMgr.hpp"
#include "AnasaziBasicEigenproblem.hpp"
#include "AnasaziConfigDefs.hpp"
#include "AnasaziEpetraAdapter.hpp"

#include "Epetra_CrsMatrix.h"

// Include the appropriate communicator include files based on whether
// or not Trilinos was built with MPI.
#ifdef EPETRA_MPI
#  include "Epetra_MpiComm.h"
#  include <mpi.h>
#else
#  include "Epetra_SerialComm.h"
#endif
#include "Epetra_Map.h"

// Wraps an Epetra_Operator A.  Computes either y := A^T (A x) or y :=
// A (A^T x), depending on the "isTrans" argument of the constructor.
// While implementing Epetra_Operator requires implementing all its
// pure virtual methods, we only have to provide _working_
// implementations for a few of them.  The others admit
// exception-throwing stubs.
class EpetraSymOp : public Epetra_Operator {
private:
  Teuchos::RCP<const Epetra_Operator> op_;
  mutable Teuchos::RCP<Epetra_MultiVector> tempVec_;
  bool isTrans_;

public:
  EpetraSymOp (const Teuchos::RCP<const Epetra_Operator>& op, 
	       const bool isTrans = false) :
    op_ (op), isTrans_ (isTrans) {}

  virtual ~EpetraSymOp () {}

  int SetUseTranspose (bool useTranspose) {
    isTrans_ = useTranspose;
  }

  double NormInf() const {
    throw std::runtime_error ("*** NormInf() not implemented! ***");
  }

  bool HasNormInf() const {
    return false;
  }

  const char* Label() const {
    return "A^T * A";
  }

  bool UseTranspose() const {
    return isTrans_;
  }

  const Epetra_Comm& Comm() const {
    return op_->Comm();
  }

  const Epetra_Map& 
  OperatorDomainMap() const 
  {
    if (isTrans_) 
      // Applying A * A^T, so the domain is that of A^T, namely the
      // range of A.
      return op_->OperatorRangeMap();
    else
      // Applying A^T * A, so the domain is that of A.
      return op_->OperatorDomainMap();
  }

  const Epetra_Map& 
  OperatorRangeMap() const 
  {
    // The domain of A^T * A is the same as the range of A^T * A.
    // Ditto for A * A^T: the domain is the same as the range.
    if (isTrans_) 
      // Applying A * A^T, so the range is that of A.
      return op_->OperatorRangeMap();
    else
      // Applying A^T * A, so the range is that of A^T, namely the
      // domain of A.
      return op_->OperatorDomainMap();
  }

  int
  Apply (const Epetra_MultiVector &X, 
	 Epetra_MultiVector &Y) const
  {
    using Teuchos::rcp;
    using Teuchos::rcp_const_cast;
    typedef Epetra_MultiVector MV;
    int info = 0;

    // Keep the original transpose state of the operator A, so that we
    // can restore it on exit.
    const bool origTransState = op_->UseTranspose ();

    if (tempVec_.is_null()) {
      if (isTrans_)
	tempVec_ = rcp (new MV (op_->OperatorDomainMap(), Y.NumVectors()));
      else
	tempVec_ = rcp (new MV (op_->OperatorRangeMap(), Y.NumVectors()));
    }

    try {
      // Transpose the operator first if computing Y = A * (A^T * X).
      // We call SetUseTranspose() anyway because we need to check
      // whether the operator A allows applying the transpose.
      info = rcp_const_cast<Epetra_Operator> (op_)->SetUseTranspose (isTrans_); 
      if (info != 0) {
	throw std::runtime_error ("Failed to set the transpose state");
      }
      // tempVec_ := A*X or A^T*X 
      info = op_->Apply (X, *tempVec_);
      if (info != 0) { 
	throw std::runtime_error ("Failed to apply the first operator");
      }
      // (Un)transpose the operator based on value of isTrans_.
      info = rcp_const_cast<Epetra_Operator> (op_)->SetUseTranspose (! isTrans_);
      if (info!=0) { 
	throw std::runtime_error ("Failed to toggle the transpose state");
      }
      // Compute A^T*(A*X) or A*A^T.
      info = op_->Apply (*tempVec_, Y);
      if (info!=0) { 
	throw std::runtime_error ("Failed to apply the second operator");
      }
    } catch (std::runtime_error& e) {
      cerr << "EpetraSymOp::Apply: " << e.what() << endl;
    }
    // We've got a shallow copy of the operator that other parts of
    // the code may be using, so we need to be polite and restore the
    // original transpose state.
    const int info2 = 
      rcp_const_cast<Epetra_Operator> (op_)->SetUseTranspose (origTransState);
    return info || info2;
  }

  int 
  ApplyInverse (const Epetra_MultiVector &X, 
		Epetra_MultiVector &Y) const
  {
    throw std::runtime_error ("*** ApplyInverse() not implemented! ***");
  }
};

//
// Return a list of parameters to pass into Anasazi's Block
// Krylov-Schur solver.
//
// blockSize: Desired block size (number of starting vectors).
//
Teuchos::RCP<Teuchos::ParameterList> 
getParameterList (const int blockSize);

//
// Build an m by n (nonsquare) sparse matrix with entries
//
//          A(i,j) = k*(si)*(tj - 1) if i <= j
//                 = k*(tj)*(si - 1) if i  > j
//
// where si = i/(m+1) and tj = j/(n+1) and k = 1/(n+1).  We use this
// matrix to exercise Anasazi's Block Krylov-Schur solver.
//
Teuchos::RCP<Epetra_CrsMatrix> 
buildSparseMatrix (const Epetra_Comm& Comm,
                   const Epetra_Map& RowMap,
                   const Epetra_Map& ColMap,
                   const int m, 
                   const int n);

//
// The "main" driver routine.
//
int 
main (int argc, char *argv[]) 
{
  // "using" statements make the code more concise.
  using Teuchos::ParameterList;
  using Teuchos::RCP;
  using Teuchos::rcp;
  typedef Teuchos::ScalarTraits<double> STS;

  // Shorthand for a collection of operations on the Scalar=double type.
  typedef Teuchos::ScalarTraits<double> STS;

  // If Trilinos was built with MPI, initialize MPI, otherwise
  // initialize the serial "communicator" that stands in for MPI.
#ifdef EPETRA_MPI
  MPI_Init (&argc,&argv);
  Epetra_MpiComm Comm (MPI_COMM_WORLD);
#else
  Epetra_SerialComm Comm;
#endif

  // Number of singular values to compute.
  const int numSingularValues = 4; 

  // The number of columns in the starting (multi)vector.
  const int blockSize = 1; 

  // Dimensions of the matrix A: m rows and n columns.
  int m = 500;
  int n = 100;

  // Construct a Map that puts approximately the same number of rows
  // of the matrix A on each processor.
  Epetra_Map RowMap(m, 0, Comm);
  Epetra_Map ColMap(n, 0, Comm);

  // Create an Epetra_CrsMatrix using the above row and column maps,
  // and the given matrix dimensions m and n.
  RCP<Epetra_CrsMatrix> A = buildSparseMatrix (Comm, RowMap, ColMap, m, n);

  //////////////////////////////////////////////////////////////////////
  //
  // Start the Block Krylov-Schur iteration.
  //
  //////////////////////////////////////////////////////////////////////

  // Get solver parameters for Block Krylov-Schur.
  RCP<ParameterList> plist = getParameterList (blockSize);

  // Set the data type of the matrix operator and multivectors.
  //
  // Anasazi solvers work with many different operator and multivector
  // types.  You may use Epetra, Tpetra, or Thyra types directly, or
  // even your own type.  If you use your own MV and OP types, you
  // must either make them inherit from Anasazi::MultiVec
  // resp. Anasazi::Operator, or specialize the
  // Anasazi::MultiVecTraits and Anasazi::OperatorTraits traits
  // classes for your MV and OP types.  Either of these options is
  // about the same amount of effort.
  typedef Epetra_MultiVector MV; 
  typedef Epetra_Operator OP;

  // The MultiVecTraits traits class maps from the Scalar and MV (see
  // above) types to a common set of operations on MV objects.  It is
  // the preferred public interface of all MV objects that Anasazi
  // accepts.  There is a corresponding OperatorTraits traits class
  // for OP objects as well, that provides a minimal common set of
  // operations for OP objects.
  typedef Anasazi::MultiVecTraits<double, MV> MVT;
  typedef Anasazi::OperatorTraits<double, MV, OP> OPT;

  // Create a vector to be the solver's starting vector.
  //
  // The starting vector must have the same number of columns as the
  // "Block Size" parameter's value.  In this case, the initial vector
  // is in the domain of the matrix A, so it uses A's domain map as
  // its map.
  RCP<MV> initVec = rcp (new MV (A->OperatorDomainMap(), blockSize));

  // Fill the initial vector with random data.
  MVT::MvRandom (*initVec);

  // Create the operator that uses the matrix A to represent the A^T *
  // A operator.  ATA does not copy the matrix A; it just keeps a
  // pointer.  An EpetraSymOp "is-an" (inherits from) Epetra_Operator.
  RCP<OP> ATA = rcp (new EpetraSymOp (A));

  // Create the object that holds the eigenvalue problem to solve.
  // The problem object is templated on the Scalar, MultiVector (MV),
  // and Operator (OP) types.
  RCP<Anasazi::BasicEigenproblem<double, MV, OP> > problem =
    rcp (new Anasazi::BasicEigenproblem<double, MV, OP> (ATA, initVec));

  // Inform the eigenproblem that the operator (A^T * A) is symmetric.
  problem->setHermitian (true);

  // Set the number of eigenvalues (singular values, in this case) to
  // compute.
  problem->setNEV (numSingularValues);

  // Inform the eigenvalue problem that you are finished passing it
  // information.
  TEST_FOR_EXCEPTION( ! problem->setProblem(), 
                      std::runtime_error,
                      "Failed to set the eigenvalue problem." );

  // Initialize the Block Krylov-Schur solver.  The solver may fill in
  // the given parameter list with defaults for any parameters that
  // were not supplied.  Thus, you don't have to know all the
  // parameters, just those that matter to you.
  Anasazi::BlockKrylovSchurSolMgr<double, MV, OP> MySolverMgr (problem, *plist);

  // Solve the problem to the specified tolerance or number of iterations.
  Anasazi::ReturnType returnCode = MySolverMgr.solve();

  if (returnCode != Anasazi::Converged && Comm.MyPID() == 0) {
    cout << "The Anasazi solver's solve() routine returned Unconverged." 
         << endl;
  }

  // Get the eigenvalues and eigenvectors from the eigenproblem.
  Anasazi::Eigensolution<double, MV> sol = problem->getSolution();

  // Anasazi stores the eigenvalues of a real matrix as pairs of real
  // values, using the Value struct.  Our operator A^T * A is
  // symmetric, so the eigenvalues should all be real.  However,
  // storing eigenvalues as (real, imaginary) pairs allows us to solve
  // nonsymmetric problems as well.
  std::vector<Anasazi::Value<double> > evals = sol.Evals;
  int computedNumSingularValues = sol.numVecs;

  if (computedNumSingularValues > 0) 
    {
      const double one = 1.0;
      const double zero = 0.0;
      // My (MPI) process rank.
      int MyPID = Comm.MyPID();

      //////////////////////////////////////////////////////////////////////
      //
      // Compute the singular values, singular vectors, and direct residuals.
      //
      //////////////////////////////////////////////////////////////////////    

      // The singular values of the matrix A are the square roots of the
      // eigenvalues of A^T * A.
      if (MyPID == 0) {
	cout << "------------------------------------------------------" << endl;
	cout << "Computed Singular Values: " << endl;
	cout << "------------------------------------------------------" << endl;
      }
      for (int i = 0; i < computedNumSingularValues; ++i) {
	// The operator A^T * A is symmetric, so the eigenvalues
	// should all have zero imaginary parts.
	evals[i].realpart = STS::squareroot (evals[i].realpart); 
      }

      //////////////////////////////////////////////////////////////////////
      // Compute the left singular vectors: u_j = (Av)_j / sigma_j
      //////////////////////////////////////////////////////////////////////

      std::vector<double> tempnrm (computedNumSingularValues);
      std::vector<double> directnrm (computedNumSingularValues);

      // The multivectors Av and u are in the range of the operator A,
      // so we have to create them using the range map of A.
      Epetra_MultiVector Av (A->OperatorRangeMap(), computedNumSingularValues);
      Epetra_MultiVector u (A->OperatorRangeMap(), computedNumSingularValues);

      // Make a view of the eigenvectors of interest.
      //
      // The current (10.6) release of Trilinos refers to a range of
      // columns of a multivector using an std::vector of their
      // (zero-based) column indices.  The development branch, and the
      // next release, will also allow using a Teuchos::Range1D to refer
      // to a range of columns of a multivector.
      std::vector<int> index (computedNumSingularValues);
      for (int i = 0; i < computedNumSingularValues; ++i) { 
	index[i] = i; 
      }
      // CloneView() produces a non-modifying ("const") view, so we have
      // to store that view in a pointer that contains a "const MV"
      // object.
      RCP<const MV> evecs = MVT::CloneView (*(sol.Evecs), index);
      OPT::Apply (*A, *evecs, Av); // Av = A * evecs
      MVT::MvNorm (Av, tempnrm);   // tempnrm_j := || (Av)_j ||_2

      Teuchos::SerialDenseMatrix<int,double> S (computedNumSingularValues, 
						computedNumSingularValues);
      for (int i = 0; i < computedNumSingularValues; ++i) { 
	S(i,i) = one / tempnrm[i]; 
      }
      // u := Av * S + zero * u
      MVT::MvTimesMatAddMv (one, Av, S, zero, u);

      //////////////////////////////////////////////////////////////////////
      // Compute direct residuals: || (Av - sigma*u)_j ||_2
      //////////////////////////////////////////////////////////////////////

      for (int i = 0; i < computedNumSingularValues; ++i) { 
	S(i,i) = evals[i].realpart; 
      }
      // Av := -one * u * S + one * Av
      //     = Av - u*S.
      MVT::MvTimesMatAddMv (-one, u, S, one, Av);

      // directnrm_j = || (Av)_j ||_2
      MVT::MvNorm (Av, directnrm);

      //////////////////////////////////////////////////////////////////////
      // Print results to stdout on MPI process 0.
      //////////////////////////////////////////////////////////////////////

      if (MyPID == 0) {
	// It's rude to set the ostream flags without restoring their 
	// original values when you're done.
	std::ios_base::fmtflags originalFlags = cout.flags ();

	// Set the flags on cout for nice neat output.
	cout.setf (std::ios_base::right, std::ios_base::adjustfield);
	cout << std::setw(16) << "Singular Value"
	     << std::setw(20) << "Direct Residual"
	     << endl;
	cout << "------------------------------------------------------" << endl;
	for (int i = 0; i < computedNumSingularValues; ++i) {
	  cout << std::setw(16) << evals[i].realpart
	       << std::setw(20) << directnrm[i] 
	       << endl;
	}  
	cout << "------------------------------------------------------" << endl;

	// Restore cout's original flags.
	cout.flags (originalFlags);
      }
    }

#ifdef EPETRA_MPI
  MPI_Finalize() ;
#endif
  return 0;
}


Teuchos::RCP<Epetra_CrsMatrix> 
buildSparseMatrix (const Epetra_Comm& Comm,
                   const Epetra_Map& RowMap,
                   const Epetra_Map& ColMap,
                   const int m, 
                   const int n)
{
  using Teuchos::RCP;
  using Teuchos::rcp;
  typedef Teuchos::ScalarTraits<double> STS;

  const double one = STS::one();
  const double zero = STS::zero();

  // My MPI process rank.
  const int MyPID = Comm.MyPID();

  // Get update list and number of local equations from newly created Map.
  const int NumMyRowElements = RowMap.NumMyElements ();
  std::vector<int> MyGlobalRowElements (NumMyRowElements);
  RowMap.MyGlobalElements (&MyGlobalRowElements[0]);

  // Create an Epetra_CrsMatrix using the given row map.
  RCP<Epetra_CrsMatrix> A = rcp (new Epetra_CrsMatrix (Copy, RowMap, n));

  // We use info to catch any errors that may have happened during
  // matrix assembly, and report them globally.  We do this so that
  // the MPI processes won't call FillComplete() unless they all
  // successfully filled their parts of the matrix.
  int info = 0;
  try {
    //
    // Compute coefficients for the discrete integral operator.
    //
    std::vector<double> Values (n);
    std::vector<int> Indices (n);
    const double inv_mp1 = one / (m+1);
    const double inv_np1 = one / (n+1);
    for (int i = 0; i < n; ++i) { 
      Indices[i] = i; 
    }
    for (int i = 0; i < NumMyRowElements; ++i) {
      for (int j = 0; j < n; ++j) 
        {
          if (MyGlobalRowElements[i] <= j) {
            Values[j] = inv_np1 * 
              ( (MyGlobalRowElements[i]+one)*inv_mp1 ) * 
              ( (j+one)*inv_np1 - one );  // k*(si)*(tj-1)
          }
          else {
            Values[j] = inv_np1 * 
              ( (j+one)*inv_np1 ) * 
              ( (MyGlobalRowElements[i]+one)*inv_mp1 - one );  // k*(tj)*(si-1)
          }
        }
      info = A->InsertGlobalValues (MyGlobalRowElements[i], n, 
                                    &Values[0], &Indices[0]);
      // Make sure that the insertion succeeded.  Teuchos'
      // TEST_FOR_EXCEPTION macro gives a nice error message if the
      // thrown exception isn't caught.  We'll report this on the
      // offending MPI process.
      TEST_FOR_EXCEPTION( info != 0, std::runtime_error, "Failed to insert n=" 
                          << n << " global value" << (n != 1 ? "s" : "") 
                          << " in row " << MyGlobalRowElements[i] 
                          << " of the matrix." );
    } // for i = 0...
    
    // Call FillComplete on the matrix.  Since the matrix isn't square,
    // we have to give FillComplete the domain and range maps, which in
    // this case are the column resp. row maps.
    info = A->FillComplete (ColMap, RowMap);
    TEST_FOR_EXCEPTION( info != 0, std::runtime_error, 
			"FillComplete failed with INFO = " << info << ".");
    info = A->OptimizeStorage();
    TEST_FOR_EXCEPTION( info != 0, std::runtime_error, 
			"OptimizeStorage failed with INFO = " << info << ".");
  } catch (std::runtime_error& e) {
    // If multiple MPI processes are reporting errors, sometimes
    // forming the error message as a string and then writing it to
    // the output stream prevents messages from different processes
    // from being interleaved.
    std::ostringstream os;
    os << "*** Error on MPI process " << MyPID << ": " << e.what();
    cerr << os.str() << endl;
    if (info == 0)
      info = -1; // All procs will share info later on.
  }

  // Do a reduction on the value of info, to ensure that all the MPI
  // processes successfully filled the sparse matrix.
  {
    int minInfo = 0;
    int maxInfo = 0;

    // Test both info < 0 and info > 0.
    Comm.MinAll (&info, &minInfo, 1);
    Comm.MaxAll (&info, &maxInfo, 1);
    TEST_FOR_EXCEPTION( minInfo != 0 || maxInfo != 0, std::runtime_error,
                        "Filling and assembling the sparse matrix failed." );
  }
  
  
  // Shut down Epetra Warning tracebacks.
  A->SetTracebackMode (1);

  return A;
}


Teuchos::RCP<Teuchos::ParameterList> 
getParameterList (const int blockSize)
{
  using Teuchos::RCP;
  using Teuchos::rcp;
  using Teuchos::ParameterList;

  RCP<ParameterList> plist = rcp (new ParameterList ("Block Krylov-Schur"));

  // "Num Blocks" parameter is the Krylov basis length to use.
  const int numBlocks = 10;

  // "Maximum Restarts" parameter is the maximum number of times to
  // restart the Krylov methods.
  const int maxRestarts = 20;

  // For this test, we set the "Convergence Tolerance" parameter to
  // machine precision for double-precision floating-point values.
  // Teuchos' ScalarTraits traits class knows what machine precision
  // is for many different floating-point types.  We could also call
  // LAPACK's DLAMCH('E') to get machine precision.
  const double tol = Teuchos::ScalarTraits<double>::eps();

  // The "Which" parameter governs the order in which Anasazi computes
  // eigenvalues (or singular values, in this case).  "LM" means the
  // eigenvalues of Largest Magnitude.  Valid options are "SM", "LM",
  // "SR", and "LR".  These mean, respectively, "Smallest Magnitude,"
  // "Largest Magnitude," "Smallest Real", and "Largest Real."  
  //
  // These abbreviations are exactly the same as those used by
  // ARPACK's "which" parameter, and in turn by Matlab's "eigs"
  // function (which calls ARPACK internally).
  const std::string which = "LM";

  plist->set ("Which", which);
  plist->set ("Block Size", blockSize);
  plist->set ("Num Blocks", numBlocks);
  plist->set ("Maximum Restarts", maxRestarts);
  plist->set ("Convergence Tolerance", tol);
  // Tell Anasazi to print output for errors, warnings, timing
  // details, and the final summary of results.
  plist->set ("Verbosity",  
              Anasazi::Errors + 
              Anasazi::Warnings + 
              Anasazi::TimingDetails + 
              Anasazi::FinalSummary);

  return plist;
}