Stratimikos Package Browser (Single Doxygen Collection) Version of the Day
belos_tpetra_thyra_lowsf_prec_lap.cpp
Go to the documentation of this file.
00001 
00002 // @HEADER
00003 // ***********************************************************************
00004 //
00005 //                 Belos: Block Linear Solvers Package
00006 //                 Copyright (2004) Sandia Corporation
00007 //
00008 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00009 // license for use of this work by or on behalf of the U.S. Government.
00010 //
00011 // This library is free software; you can redistribute it and/or modify
00012 // it under the terms of the GNU Lesser General Public License as
00013 // published by the Free Software Foundation; either version 2.1 of the
00014 // License, or (at your option) any later version.
00015 //
00016 // This library is distributed in the hope that it will be useful, but
00017 // WITHOUT ANY WARRANTY; without even the implied warranty of
00018 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00019 // Lesser General Public License for more details.
00020 //
00021 // You should have received a copy of the GNU Lesser General Public
00022 // License along with this library; if not, write to the Free Software
00023 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00024 // USA
00025 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00026 //
00027 // ***********************************************************************
00028 // @HEADER
00029 //
00030 // This driver creates a 2D Laplacian operator as a Tpetra::CisMatrix
00031 // and a preconditioner.  These matrices are then converted into Thyra 
00032 // linear operators through the Thyra-Tpetra adapters.
00033 //
00034 // The right-hand sides are all random vectors.
00035 // The initial guesses are all set to zero. 
00036 //
00037 //
00038 
00039 // Thyra includes
00040 #include "Thyra_BelosLinearOpWithSolveFactory.hpp"
00041 #include "Thyra_LinearOpWithSolveFactoryHelpers.hpp"
00042 #include "Thyra_PreconditionerFactoryHelpers.hpp"
00043 #include "Thyra_TpetraLinearOp.hpp"
00044 #include "Thyra_DefaultPreconditioner.hpp"
00045 
00046 // Tpetra includes
00047 #include "Tpetra_ElementSpace.hpp"
00048 #include "Tpetra_VectorSpace.hpp"
00049 #include "Tpetra_CisMatrix.hpp"
00050 
00051 // Teuchos includes
00052 #include "Teuchos_ParameterList.hpp"
00053 #include "Teuchos_VerboseObject.hpp"
00054 #include "Teuchos_GlobalMPISession.hpp"
00055 
00056 #ifdef HAVE_MPI
00057 #  include "Tpetra_MpiPlatform.hpp"
00058 #else
00059 #  include "Tpetra_SerialPlatform.hpp"
00060 #endif
00061 
00062 int main(int argc, char* argv[])
00063 {
00064   //
00065   // Get a default output stream from the Teuchos::VerboseObjectBase
00066   //
00067   Teuchos::RCP<Teuchos::FancyOStream>
00068     out = Teuchos::VerboseObjectBase::getDefaultOStream();
00069   
00070   Teuchos::GlobalMPISession mpiSession(&argc,&argv);
00071 
00072   //
00073   // Create the scalar-type typedefs
00074   //
00075 
00076 #ifdef HAVE_COMPLEX
00077   typedef std::complex<double> ST;  // Scalar-type typedef
00078 #elif HAVE_COMPLEX_H
00079   typedef std::complex<double> ST;     // Scalar-type typedef
00080 #else
00081   typedef double ST;                // Scalar-type typedef
00082 #endif
00083   typedef Teuchos::ScalarTraits<ST>::magnitudeType MT;  // Magnitude-type typdef
00084   typedef int OT;                   // Ordinal-type typdef
00085 
00086   ST one = Teuchos::ScalarTraits<ST>::one();
00087   ST zero = Teuchos::ScalarTraits<ST>::zero();
00088   
00089   //
00090   // Create 2D Laplacian using Tpetra::CisMatrix
00091   //
00092 
00093 #ifdef HAVE_MPI
00094   MPI_Comm mpiComm = MPI_COMM_WORLD;
00095   const Tpetra::MpiPlatform<OT,OT>  ordinalPlatform(mpiComm);
00096   const Tpetra::MpiPlatform<OT,ST>   scalarPlatform(mpiComm);
00097 #else
00098   const Tpetra::SerialPlatform<OT,OT>  ordinalPlatform;
00099   const Tpetra::SerialPlatform<OT,ST>   scalarPlatform;
00100 #endif
00101 
00102   // Declare global dimension of the linear operator
00103   OT globalDim = 500;
00104 
00105   // Create the element space and std::vector space
00106   const Tpetra::ElementSpace<OT> elementSpace(globalDim,0,ordinalPlatform);
00107   const Tpetra::VectorSpace<OT,ST> vectorSpace(elementSpace,scalarPlatform);
00108   
00109   // Allocate the Tpetra::CisMatrix object.
00110   RCP<Tpetra::CisMatrix<OT,ST> >
00111     tpetra_A = rcp(new Tpetra::CisMatrix<OT,ST>(vectorSpace));
00112 
00113   // Allocate the Tpetra::CisMatrix preconditioning object.
00114   RCP<Tpetra::CisMatrix<OT,ST> >
00115     tpetra_Prec = rcp(new Tpetra::CisMatrix<OT,ST>(vectorSpace));
00116 
00117   // Get the indexes of the rows on this processor
00118   const int numMyElements = vectorSpace.getNumMyEntries();
00119   const std::vector<int> &myGlobalElements = vectorSpace.elementSpace().getMyGlobalElements();
00120 
00121   // Fill the local matrix entries one row at a time.
00122   const ST offDiag = -1.0, diag = 2.0;
00123   int numEntries; ST values[3]; int indexes[3];
00124   ST prec_value[1]; int prec_index[1];
00125   prec_value[0] = one/diag;   // Diagonal preconditioning
00126   for( int k = 0; k < numMyElements; ++k ) {
00127     const int rowIndex = myGlobalElements[k];
00128     prec_index[0] = rowIndex;
00129     if( rowIndex == 0 ) {                     // First row
00130       numEntries = 2;
00131       values[0]  = diag;             values[1]  = offDiag;
00132       indexes[0] = 0;                indexes[1] = 1; 
00133     }
00134     else if( rowIndex == globalDim - 1 ) {    // Last row
00135       numEntries = 2;
00136       values[0]  = offDiag;         values[1]  = diag;
00137       indexes[0] = globalDim-2;     indexes[1] = globalDim-1; 
00138     }
00139     else {                                    // Middle rows
00140       numEntries = 3;
00141       values[0]  = offDiag;         values[1]  = diag;          values[2]  = offDiag;
00142       indexes[0] = rowIndex-1;      indexes[1] = rowIndex;      indexes[2] = rowIndex+1; 
00143     }
00144     tpetra_A->submitEntries(Tpetra::Insert,rowIndex,numEntries,values,indexes);
00145     tpetra_Prec->submitEntries(Tpetra::Insert,rowIndex,1,prec_value,prec_index);
00146   }
00147 
00148   // Finish the construction of the Tpetra::CisMatrix
00149   tpetra_A->fillComplete();
00150   tpetra_Prec->fillComplete();
00151 
00152   // Create a Thyra linear operator (A) using the Tpetra::CisMatrix (tpetra_A).
00153   RCP<Thyra::LinearOpBase<ST> >
00154     A = Teuchos::rcp( new Thyra::TpetraLinearOp<OT,ST>(
00155                       Teuchos::rcp_implicit_cast<Tpetra::Operator<OT,ST> >(tpetra_A) ) );
00156 
00157   // Create a Thyra linear operator (Prec) using the Tpetra::CisMatrix (tpetra_Prec)
00158   RCP<Thyra::LinearOpBase<ST> >
00159     Prec = Teuchos::rcp( new Thyra::TpetraLinearOp<OT,ST>(
00160                     Teuchos::rcp_implicit_cast<Tpetra::Operator<OT,ST> >(tpetra_Prec) ) );
00161 
00162   // Create a Thyra default preconditioner (DefPrec) using the Thyra linear operator (Prec)
00163   RCP<const Thyra::DefaultPreconditioner<ST> >
00164     DefPrec = Thyra::unspecifiedPrec<ST>(Prec);
00165 
00166   //
00167   // Set the parameters for the Belos LOWS Factory and create a parameter list.
00168   // NOTE:  All the code below only uses Thyra and is independent of Tpetra.
00169   //
00170   int             blockSize              = 1;  // can only be 1 right now.
00171   int             maxIterations          = globalDim;
00172   int             maxRestarts            = 0;
00173   int             gmresKrylovLength      = globalDim;
00174   int             outputFrequency        = 100;
00175   bool            outputMaxResOnly       = true;
00176   MT              maxResid               = 1e-3;
00177 
00178   Teuchos::RCP<Teuchos::ParameterList>
00179     belosLOWSFPL = Teuchos::rcp( new Teuchos::ParameterList() );
00180  
00181   belosLOWSFPL->set("Solver Type","Block GMRES");
00182 
00183   Teuchos::ParameterList& belosLOWSFPL_solver =
00184     belosLOWSFPL->sublist("Solver Types");
00185 
00186   Teuchos::ParameterList& belosLOWSFPL_gmres =
00187     belosLOWSFPL_solver.sublist("Block GMRES");
00188 
00189   belosLOWSFPL_gmres.set("Maximum Iterations",int(maxIterations));
00190   belosLOWSFPL_gmres.set("Convergence Tolerance",double(maxResid));
00191   belosLOWSFPL_gmres.set("Maximum Restarts",int(maxRestarts));
00192   belosLOWSFPL_gmres.set("Block Size",int(blockSize));
00193   belosLOWSFPL_gmres.set("Num Blocks",int(gmresKrylovLength));
00194   belosLOWSFPL_gmres.set("Output Frequency",int(outputFrequency));
00195   belosLOWSFPL_gmres.set("Show Maximum Residual Norm Only",bool(outputMaxResOnly));
00196  
00197   // Whether the linear solver succeeded.
00198   // (this will be set during the residual check at the end)
00199   bool success = true;
00200 
00201   // Number of random right-hand sides we will be solving for.
00202   int numRhs = 5;
00203 
00204   // Get the domain space for the Thyra linear operator 
00205   Teuchos::RCP<const Thyra::VectorSpaceBase<ST> > domain = A->domain();
00206 
00207   // Create the Belos LOWS factory.
00208   Teuchos::RCP<Thyra::LinearOpWithSolveFactoryBase<ST> >
00209     belosLOWSFactory = Teuchos::rcp(new Thyra::BelosLinearOpWithSolveFactory<ST>());
00210 
00211   // Set the parameter list to specify the behavior of the factory.
00212   belosLOWSFactory->setParameterList( belosLOWSFPL );
00213 
00214   // Set the output stream and the verbosity level (prints to std::cout by defualt)
00215   // NOTE:  Set to VERB_NONE for no output from the solver.
00216   belosLOWSFactory->setVerbLevel(Teuchos::VERB_LOW);
00217 
00218   // Create a BelosLinearOpWithSolve object from the Belos LOWS factory.
00219   Teuchos::RCP<Thyra::LinearOpWithSolveBase<ST> >
00220     nsA = belosLOWSFactory->createOp();
00221 
00222   // Initialize the BelosLinearOpWithSolve object with the Thyra linear operator (A)
00223   // and preconditioner (DefPrec).
00224   Thyra::initializeOp<ST>( *belosLOWSFactory, A, &*nsA );
00225 
00226   // Create a right-hand side with numRhs vectors in it and randomize it.
00227   Teuchos::RCP< Thyra::MultiVectorBase<ST> > 
00228     b = Thyra::createMembers(domain, numRhs);
00229   Thyra::seed_randomize<ST>(0);
00230   Thyra::randomize(-one, one, &*b);
00231 
00232   // Create an initial std::vector with numRhs vectors in it and initialize it to zero.
00233   Teuchos::RCP< Thyra::MultiVectorBase<ST> >
00234     x = Thyra::createMembers(domain, numRhs);
00235   Thyra::assign(&*x, zero);
00236 
00237   // Perform solve using the linear operator to get the approximate solution of Ax=b,
00238   // where b is the right-hand side and x is the left-hand side.
00239   Thyra::SolveStatus<ST> solveStatus;
00240   solveStatus = Thyra::solve( *nsA, Thyra::NONCONJ_ELE, *b, &*x );
00241 
00242   // Print out status of solve.
00243   *out << "\nBelos LOWS Status: "<< solveStatus << std::endl;
00244 
00245   //
00246   // Compute residual and ST check convergence.
00247   //
00248   std::vector<MT> norm_b(numRhs), norm_res(numRhs);
00249   Teuchos::RCP< Thyra::MultiVectorBase<ST> >
00250     y = Thyra::createMembers(domain, numRhs);
00251 
00252   // Compute the column norms of the right-hand side b.
00253   Thyra::norms_2( *b, &norm_b[0] );
00254 
00255   // Compute y=A*x, where x is the solution from the linear solver.
00256   A->apply( Thyra::NONCONJ_ELE, *x, &*y );
00257   
00258   // Compute A*x-b = y-b
00259   Thyra::update( -one, *b, &*y );
00260 
00261   // Compute the column norms of A*x-b.
00262   Thyra::norms_2( *y, &norm_res[0] );
00263 
00264   // Print out the final relative residual norms.
00265   MT rel_res = 0.0;
00266   *out << "Final relative residual norms" << std::endl;  
00267   for (int i=0; i<numRhs; ++i) {
00268     rel_res = norm_res[i]/norm_b[i];
00269     if (rel_res > maxResid)
00270       success = false;
00271     *out << "RHS " << i+1 << " : " 
00272          << std::setw(16) << std::right << rel_res << std::endl;
00273   }
00274 
00275   return ( success ? 0 : 1 );
00276 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines