doc/html/BlockKrylovSchur_2BlockKrylovSchurEpetraExGenBelos_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 the

//  discretized 2D Laplacian operator using the block Krylov-Schur method.

//  This problem shows the construction of an inner-outer iteration using

//  Belos as the linear solver within Anasazi.  An Ifpack preconditioner

//  is constructed to precondition the linear solver.  This operator is

//  discretized using linear finite elements and constructed as an Epetra

//  matrix, then passed into the Belos solver to perform the shift-invert

//  operation to be used within the Krylov decomposition.  The specifics

//  of the block Krylov-Schur method can be set by the user.


// Include autoconfigured header

#include "AnasaziConfigDefs.hpp"


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

#include "AnasaziEpetraAdapter.hpp"


// Include header for Epetra compressed-row storage matrix and linear problem

#include "Epetra_CrsMatrix.h"

#include "Epetra_LinearProblem.h"

#include "Epetra_InvOperator.h"


// Include header for Belos solver and solver interface for Epetra_Operator

#include "BelosEpetraOperator.h"

#include "BelosEpetraAdapter.hpp"


// Include header for Ifpack preconditioner factory

#include "Ifpack.h"


// Include header for Teuchos serial dense matrix

#include "Teuchos_SerialDenseMatrix.hpp"


// Include header for the problem definition

#include "ModeLaplace2DQ2.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"


int main(int argc, char *argv[]) {

  using std::cout;

  using std::endl;

  int i;


#ifdef EPETRA_MPI

  // Initialize MPI

  MPI_Init(&argc,&argv);

  Epetra_MpiComm Comm(MPI_COMM_WORLD);

#else

  Epetra_SerialComm Comm;

#endif


  int MyPID = Comm.MyPID();


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


  //

  // ************Construct preconditioner*************

  //

  Teuchos::ParameterList ifpackList;


  // allocates an IFPACK factory. No data is associated

  // to this object (only method Create()).

  Ifpack Factory;


  // create the preconditioner. For valid PrecType values,

  // please check the documentation

  std::string PrecType = "ICT"; // incomplete Cholesky

  int OverlapLevel = 0; // must be >= 0. If Comm.NumProc() == 1,

                        // it is ignored.


  Teuchos::RCP<Ifpack_Preconditioner> Prec = Teuchos::rcp( Factory.Create(PrecType, &*K, OverlapLevel) );

  assert(Prec != Teuchos::null);


  // specify parameters for ICT

  ifpackList.set("fact: drop tolerance", 1e-4);

  ifpackList.set("fact: ict level-of-fill", 0.);

  // the combine mode is on the following:

  // "Add", "Zero", "Insert", "InsertAdd", "Average", "AbsMax"

  // Their meaning is as defined in file Epetra_CombineMode.h

  ifpackList.set("schwarz: combine mode", "Add");

  // sets the parameters

  IFPACK_CHK_ERR(Prec->SetParameters(ifpackList));


  // initialize the preconditioner. At this point the matrix must

  // have been FillComplete()'d, but actual values are ignored.

  IFPACK_CHK_ERR(Prec->Initialize());


  // Builds the preconditioners, by looking for the values of

  // the matrix.

  IFPACK_CHK_ERR(Prec->Compute());


  //

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

  // Set up Belos Block CG operator for inner iteration

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

  //

  int blockSize = 3; // block size used by linear solver and eigensolver [ not required to be the same ]

  int maxits = K->NumGlobalRows(); // maximum number of iterations to run

  //

  // Create the Belos::LinearProblem

  //

  Teuchos::RCP<Belos::LinearProblem<double,Epetra_MultiVector,Epetra_Operator> >

    My_LP = Teuchos::rcp( new Belos::LinearProblem<double,Epetra_MultiVector,Epetra_Operator>() );

  My_LP->setOperator( K );


  // Create the Belos preconditioned operator from the Ifpack preconditioner.

  // NOTE:  This is necessary because Belos expects an operator to apply the

  //        preconditioner with Apply() NOT ApplyInverse().

  Teuchos::RCP<Epetra_Operator> belosPrec = Teuchos::rcp( new Epetra_InvOperator( Prec.get() ) );

  My_LP->setLeftPrec( belosPrec );

  //

  // Create the ParameterList for the Belos Operator

  //

  Teuchos::RCP<Teuchos::ParameterList> My_List = Teuchos::rcp( new Teuchos::ParameterList() );

  My_List->set( "Solver", "BlockCG" );

  My_List->set( "Maximum Iterations", maxits );

  My_List->set( "Block Size", 1 );

  My_List->set( "Convergence Tolerance", 1e-12 );

  //

  // Create the Belos::EpetraOperator

  //

  Teuchos::RCP<Belos::EpetraOperator> BelosOp =

    Teuchos::rcp( new Belos::EpetraOperator( My_LP, My_List ));

  //

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

  // Start the block Arnoldi iteration

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

  //

  //  Variables used for the Block Arnoldi Method

  //

  double tol = 1.0e-8;

  int nev = 10;

  int numBlocks = 3*nev/blockSize;

  int maxRestarts = 5;

  //int step = 5;

  std::string which = "LM";

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

  //

  // 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( "Num Blocks", numBlocks );

  MyPL.set( "Maximum Restarts", maxRestarts );

  MyPL.set( "Convergence Tolerance", tol );

  //MyPL.set( "Step Size", step );


  typedef Epetra_MultiVector MV;

  typedef Epetra_Operator OP;

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

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


  // Call the ctor that calls the petra ctor for a matrix

  Teuchos::RCP<Anasazi::EpetraGenOp> Aop = Teuchos::rcp( new Anasazi::EpetraGenOp(BelosOp, M, false) );


  Teuchos::RCP<Anasazi::BasicEigenproblem<double,MV,OP> > MyProblem =

    Teuchos::rcp( new Anasazi::BasicEigenproblem<double,MV,OP>(Aop, 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 != true) {

    if (MyPID == 0) {

      cout << "Anasazi::BasicEigenproblem::setProblem() returned with error." << endl;

    }

#ifdef HAVE_MPI

    MPI_Finalize() ;

#endif

    return -1;

  }


  // Initialize the Block Arnoldi solver

  Anasazi::BlockKrylovSchurSolMgr<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) {

    cout << "Anasazi::EigensolverMgr::solve() returned unconverged." << 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) {


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

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

    OPT::Apply( *K, *evecs, tempvec );

    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 (i=0; i<numev; i++) {

        compeval = dmatr(i,i);

        cout<<std::setw(16)<<compeval

          <<std::setw(16)<<Teuchos::ScalarTraits<double>::magnitude(compeval-1.0/evals[i].realpart)

          <<endl;

      }

      cout<<"------------------------------------------------------"<<endl;

    }


  }


#ifdef EPETRA_MPI

  MPI_Finalize();

#endif


  return 0;

}