Stratimikos Package Browser (Single Doxygen Collection) Version of the Day
belos_tpetra_thyra_lowsf_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 // This matrix is then converted into a Thyra linear operator 
00032 // 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_TpetraLinearOp.hpp"
00043 
00044 // Tpetra includes
00045 #include "Tpetra_ElementSpace.hpp"
00046 #include "Tpetra_VectorSpace.hpp"
00047 #include "Tpetra_CisMatrix.hpp"
00048 
00049 // Teuchos includes
00050 #include "Teuchos_ParameterList.hpp"
00051 #include "Teuchos_VerboseObject.hpp"
00052 #include "Teuchos_GlobalMPISession.hpp"
00053 
00054 #ifdef HAVE_MPI
00055 #  include "Tpetra_MpiPlatform.hpp"
00056 #else
00057 #  include "Tpetra_SerialPlatform.hpp"
00058 #endif
00059 
00060 int main(int argc, char* argv[])
00061 {
00062   //
00063   // Get a default output stream from the Teuchos::VerboseObjectBase
00064   //
00065   Teuchos::RCP<Teuchos::FancyOStream>
00066     out = Teuchos::VerboseObjectBase::getDefaultOStream();
00067   
00068   Teuchos::GlobalMPISession mpiSession(&argc,&argv);
00069 
00070   // --------------------------------------------------------------------------------
00071   //
00072   // Create the scalar-type typedefs
00073   //
00074   // --------------------------------------------------------------------------------
00075 
00076 #ifdef HAVE_COMPLEX
00077   typedef std::complex<double> ST;
00078 #elif HAVE_COMPLEX_H
00079   typedef std::complex<double> ST;
00080 #else
00081   typedef double ST;
00082 #endif
00083   typedef Teuchos::ScalarTraits<ST>::magnitudeType MT;
00084   typedef int OT;
00085 
00086   // --------------------------------------------------------------------------------
00087   //
00088   // Create 2D Laplacian using Tpetra::CisMatrix
00089   //
00090   // --------------------------------------------------------------------------------
00091 
00092 #ifdef HAVE_MPI
00093   MPI_Comm mpiComm = MPI_COMM_WORLD;
00094   const Tpetra::MpiPlatform<OT,OT>  ordinalPlatform(mpiComm);
00095   const Tpetra::MpiPlatform<OT,ST>   scalarPlatform(mpiComm);
00096 #else
00097   const Tpetra::SerialPlatform<OT,OT>  ordinalPlatform;
00098   const Tpetra::SerialPlatform<OT,ST>   scalarPlatform;
00099 #endif
00100 
00101   // Declare global dimension of the linear operator
00102   OT globalDim = 500;
00103 
00104   // Create the element space and std::vector space
00105   const Tpetra::ElementSpace<OT> elementSpace(globalDim,0,ordinalPlatform);
00106   const Tpetra::VectorSpace<OT,ST> vectorSpace(elementSpace,scalarPlatform);
00107   
00108   // Allocate the Tpetra::CisMatrix object.
00109   RCP<Tpetra::CisMatrix<OT,ST> >
00110     tpetra_A = rcp(new Tpetra::CisMatrix<OT,ST>(vectorSpace));
00111 
00112   // Get the indexes of the rows on this processor
00113   const int numMyElements = vectorSpace.getNumMyEntries();
00114   const std::vector<int> &myGlobalElements = vectorSpace.elementSpace().getMyGlobalElements();
00115 
00116   // Fill the local matrix entries one row at a time.
00117   const ST offDiag = -1.0, diag = 2.0;
00118   int numEntries; ST values[3]; int indexes[3];
00119   for( int k = 0; k < numMyElements; ++k ) {
00120     const int rowIndex = myGlobalElements[k];
00121     if( rowIndex == 0 ) {                     // First row
00122       numEntries = 2;
00123       values[0]  = diag;             values[1]  = offDiag;
00124       indexes[0] = 0;                indexes[1] = 1; 
00125     }
00126     else if( rowIndex == globalDim - 1 ) {    // Last row
00127       numEntries = 2;
00128       values[0]  = offDiag;         values[1]  = diag;
00129       indexes[0] = globalDim-2;     indexes[1] = globalDim-1; 
00130     }
00131     else {                                    // Middle rows
00132       numEntries = 3;
00133       values[0]  = offDiag;         values[1]  = diag;          values[2]  = offDiag;
00134       indexes[0] = rowIndex-1;      indexes[1] = rowIndex;      indexes[2] = rowIndex+1; 
00135     }
00136     tpetra_A->submitEntries(Tpetra::Insert,rowIndex,numEntries,values,indexes);
00137   }
00138 
00139   // Finish the construction of the Tpetra::CisMatrix
00140   tpetra_A->fillComplete();
00141 
00142   // Create a Thyra linear operator (A) using the Tpetra::CisMatrix (tpetra_A).
00143   RCP<Thyra::LinearOpBase<ST> >
00144     A = Teuchos::rcp( new Thyra::TpetraLinearOp<OT,ST>(
00145                       Teuchos::rcp_implicit_cast<Tpetra::Operator<OT,ST> >(tpetra_A) ) );
00146 
00147   // Set the parameters for the Belos LOWS Factory and create a parameter list.
00148   // NOTE:  All the code below only uses Thyra and is independent of Tpetra.
00149   //
00150   int             blockSize              = 1;  // can only be 1 right now.
00151   int             maxIterations          = globalDim;
00152   int             maxRestarts            = 0;
00153   int             gmresKrylovLength      = globalDim;
00154   int             outputFrequency        = 100;
00155   bool            outputMaxResOnly       = true;
00156   MT              maxResid               = 1e-3;
00157 
00158   ST one = Teuchos::ScalarTraits<ST>::one();
00159   ST zero = Teuchos::ScalarTraits<ST>::zero();
00160   
00161   Teuchos::RCP<Teuchos::ParameterList>
00162     belosLOWSFPL = Teuchos::rcp( new Teuchos::ParameterList() );
00163  
00164   belosLOWSFPL->set("Solver Type","Block GMRES");
00165 
00166   Teuchos::ParameterList& belosLOWSFPL_solver =
00167     belosLOWSFPL->sublist("Solver Types");
00168 
00169   Teuchos::ParameterList& belosLOWSFPL_gmres =
00170     belosLOWSFPL_solver.sublist("Block GMRES");
00171 
00172   belosLOWSFPL_gmres.set("Maximum Iterations",int(maxIterations));
00173   belosLOWSFPL_gmres.set("Convergence Tolerance",double(maxResid));
00174   belosLOWSFPL_gmres.set("Maximum Restarts",int(maxRestarts));
00175   belosLOWSFPL_gmres.set("Block Size",int(blockSize));
00176   belosLOWSFPL_gmres.set("Num Blocks",int(gmresKrylovLength));
00177   belosLOWSFPL_gmres.set("Output Frequency",int(outputFrequency));
00178   belosLOWSFPL_gmres.set("Show Maximum Residual Norm Only",bool(outputMaxResOnly));
00179  
00180   // Whether the linear solver succeeded.
00181   // (this will be set during the residual check at the end)
00182   bool success = true;
00183 
00184   // --------------------------------------------------------------------------------
00185   //
00186   // Create the right-hand sides and solution vectors
00187   //
00188   // --------------------------------------------------------------------------------
00189 
00190   // Number of random right-hand sides we will be solving for.
00191   int numRhs = 5;
00192 
00193   // Create smart pointer to right-hand side and solution std::vector to be filled in below.
00194   Teuchos::RCP<Thyra::MultiVectorBase<ST> > x, b;
00195 
00196   if (numRhs==1) {
00197     //
00198     // In this case we can construct vectors using Tpetra and just "wrap" them in Thyra objects.
00199     //
00200 
00201     // Create RHS std::vector
00202     Teuchos::RCP<Tpetra::Vector<OT,ST> > tpetra_b =
00203       Teuchos::rcp( new Tpetra::Vector<OT,ST>(vectorSpace) );
00204 
00205     // Randomize RHS std::vector
00206     tpetra_b->setAllToRandom();
00207 
00208     // Wrap Tpetra std::vector as Thyra std::vector
00209     b = Thyra::create_Vector(tpetra_b);
00210     
00211     // Create solution (LHS) std::vector
00212     Teuchos::RCP<Tpetra::Vector<OT,ST> > tpetra_x =
00213       Teuchos::rcp( new Tpetra::Vector<OT,ST>(vectorSpace) );
00214 
00215     // Initialize solution to zero
00216     tpetra_x->setAllToScalar( zero );
00217 
00218     // Wrap Tpetra std::vector as Thyra std::vector
00219     x = Thyra::create_Vector(tpetra_x);
00220 
00221   } 
00222   else {
00223     //
00224     // In this case we can construct the multivector using Thyra and extract columns of
00225     // the multivector as Tpetra std::vector, which can then be filled.  This is because
00226     // Tpetra does not have a multivector object and Thyra will emulate the multivector.
00227     //
00228     
00229     // Get the domain space for the Thyra linear operator 
00230     Teuchos::RCP<const Thyra::VectorSpaceBase<ST> > domain = A->domain();
00231     
00232     // Create a right-hand side and solution vectors with numRhs vectors in it.
00233     x = Thyra::createMembers(domain, numRhs);
00234     b = Thyra::createMembers(domain, numRhs);
00235 
00236     // Extract the Tpetra std::vector from the columns of the multivector and fill them with 
00237     // random numbers (b) or zeros (x).
00238 
00239     for ( int j=0; j<numRhs; ++j ) {
00240       //
00241       // Get the j-th column from b as a Tpetra std::vector and randomize it.
00242       // 
00243       Teuchos::RCP<Tpetra::Vector<OT,ST> > 
00244   tpetra_b_j = Thyra::get_Tpetra_Vector(vectorSpace,b->col(j));
00245       tpetra_b_j->setAllToRandom();
00246       // 
00247       // Get the j-th column from x as a Tpetra std::vector and set it to zero.
00248       //
00249       Teuchos::RCP<Tpetra::Vector<OT,ST> > 
00250   tpetra_x_j = Thyra::get_Tpetra_Vector(vectorSpace,x->col(j));
00251       tpetra_x_j->setAllToScalar( zero );
00252       //
00253       // NOTE: Tpetra vectors have element access via the [] operator.
00254       //       So and additional approach for filling in a Tpetra std::vector is:
00255       //
00256       for ( int i=0; i<numMyElements; ++i ) {
00257   //
00258   // Get the global index.
00259   //
00260   const int rowIndex = myGlobalElements[ i ];
00261 
00262   // Set the entry to zero.
00263   (*tpetra_x_j)[ rowIndex ] = zero;
00264       }
00265     }
00266   }
00267   
00268   // --------------------------------------------------------------------------------
00269   //
00270   // Create the linear operator with solve factory/object and solve the linear
00271   // system using the iterative solvers in Belos.
00272   //
00273   // NOTE:  This is the only part of the code that solely uses Thyra.
00274   //
00275   // --------------------------------------------------------------------------------
00276 
00277   // Create the Belos LOWS factory.
00278   Teuchos::RCP<Thyra::LinearOpWithSolveFactoryBase<ST> >
00279     belosLOWSFactory = Teuchos::rcp(new Thyra::BelosLinearOpWithSolveFactory<ST>());
00280   
00281   // Set the output stream and the verbosity level (prints to std::cout by defualt)
00282   // NOTE:  Set to VERB_NONE for no output from the solver.
00283   belosLOWSFactory->setVerbLevel(Teuchos::VERB_LOW);
00284   
00285   // Set the parameter list to specify the behavior of the factory.
00286   belosLOWSFactory->setParameterList( belosLOWSFPL );
00287   
00288   // Create a BelosLinearOpWithSolve object from the Belos LOWS factory.
00289   Teuchos::RCP<Thyra::LinearOpWithSolveBase<ST> >
00290     nsA = belosLOWSFactory->createOp();
00291   
00292   // Initialize the BelosLinearOpWithSolve object with the Thyra linear operator.
00293   Thyra::initializeOp<ST>( *belosLOWSFactory, A, &*nsA );
00294   
00295   // Perform solve using the linear operator to get the approximate solution of Ax=b,
00296   // where b is the right-hand side and x is the left-hand side.
00297   Thyra::SolveStatus<ST> solveStatus;
00298   solveStatus = Thyra::solve( *nsA, Thyra::NONCONJ_ELE, *b, &*x );
00299   
00300   // Print out status of solve.
00301   *out << "\nBelos LOWS Status: "<< solveStatus << std::endl;
00302   
00303   // --------------------------------------------------------------------------------
00304   // Compute residual and check convergence.
00305   // NOTE: We will do this by extracting the Tpetra vectors from the multivector.
00306   //
00307   // NOTE 2: We are using scalar traits here instead of the magnitude type
00308   //         because Tpetra defines its norm methods to return the scalar type.
00309   // --------------------------------------------------------------------------------
00310   std::vector<ST> norm_b(numRhs), norm_res(numRhs);
00311   
00312   for ( int j=0; j<numRhs; ++j ) {
00313     //
00314     // Get the j-th columns from x and b as Tpetra vectors.
00315     // 
00316     Teuchos::RCP<Tpetra::Vector<OT,ST> > 
00317       tpetra_x_j = Thyra::get_Tpetra_Vector(vectorSpace,x->col(j));
00318 
00319     Teuchos::RCP<Tpetra::Vector<OT,ST> > 
00320       tpetra_b_j = Thyra::get_Tpetra_Vector(vectorSpace,b->col(j));
00321     
00322     // Compute the column norms of the right-hand side b.
00323     norm_b[j] = tpetra_b_j->norm2();
00324       
00325     // Compute y=A*x, where x is the solution from the linear solver.
00326     Tpetra::Vector<OT,ST> y(vectorSpace);
00327     tpetra_A->apply( *tpetra_x_j, y );
00328 
00329     // Compute b-A*x = b-y
00330     y.update( one, *tpetra_b_j, -one );
00331 
00332     // Compute the column norms of A*x-b.
00333     norm_res[j] = y.norm2();
00334   }
00335 
00336   // Print out the final relative residual norms.
00337   MT rel_res = 0.0;
00338   *out << "Final relative residual norms" << std::endl;  
00339   for (int i=0; i<numRhs; ++i) {
00340     rel_res = Teuchos::ScalarTraits<ST>::real(norm_res[i]/norm_b[i]);
00341     if (rel_res > maxResid)
00342       success = false;
00343     *out << "RHS " << i+1 << " : " 
00344          << std::setw(16) << std::right << rel_res << std::endl;
00345   }
00346   
00347   return ( success ? 0 : 1 );
00348 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines