Stratimikos Package Browser (Single Doxygen Collection) Version of the Day
belos_epetra_thyra_lowsf_hb.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 reads a problem from a Harwell-Boeing (HB) file into an 
00031 // Epetra_CrsMatrix.  This matrix is then converted into a Thyra linear operator
00032 // through the Thyra-Epetra adapters.
00033 // The right-hand-side from the HB file is used instead of random vectors.
00034 // The initial guesses are all set to zero. 
00035 //
00036 //
00037 
00038 // Thyra includes
00039 #include "Thyra_BelosLinearOpWithSolveFactory.hpp"
00040 #include "Thyra_LinearOpWithSolveFactoryHelpers.hpp"
00041 #include "Thyra_EpetraLinearOp.hpp"
00042 
00043 // Epetra includes
00044 #include "Epetra_SerialComm.h"
00045 #include "EpetraExt_readEpetraLinearSystem.h"
00046 
00047 // Ifpack includes
00048 #ifdef HAVE_BELOS_IFPACK
00049   #include "Thyra_IfpackPreconditionerFactory.hpp"
00050 #endif
00051 
00052 // Teuchos includes
00053 #include "Teuchos_ParameterList.hpp"
00054 #include "Teuchos_VerboseObject.hpp"
00055 
00056 int main(int argc, char* argv[])
00057 {
00058   //
00059   // Get a default output stream from the Teuchos::VerboseObjectBase
00060   //
00061   Teuchos::RCP<Teuchos::FancyOStream>
00062     out = Teuchos::VerboseObjectBase::getDefaultOStream();
00063   //
00064   // Set the parameters for the Belos LOWS Factory and create a parameter list.
00065   //
00066   int             blockSize              = 2;
00067   int             maxIterations          = 400;
00068   int             maxRestarts            = 25;
00069   int             gmresKrylovLength      = 25;
00070   int             outputFrequency        = 1;
00071   bool            outputMaxResOnly       = true;
00072   double          maxResid               = 1e-6;
00073   
00074   Teuchos::RCP<Teuchos::ParameterList>
00075     belosLOWSFPL = Teuchos::rcp( new Teuchos::ParameterList() );
00076   
00077   belosLOWSFPL->set("Solver Type","Block GMRES");
00078   
00079   Teuchos::ParameterList& belosLOWSFPL_solver = 
00080     belosLOWSFPL->sublist("Solver Types");
00081 
00082   Teuchos::ParameterList& belosLOWSFPL_gmres = 
00083     belosLOWSFPL_solver.sublist("Block GMRES");
00084 
00085   belosLOWSFPL_gmres.set("Maximum Iterations",int(maxIterations));
00086   belosLOWSFPL_gmres.set("Convergence Tolerance",double(maxResid));
00087   belosLOWSFPL_gmres.set("Maximum Restarts",int(maxRestarts));
00088   belosLOWSFPL_gmres.set("Block Size",int(blockSize));
00089   belosLOWSFPL_gmres.set("Num Blocks",int(gmresKrylovLength));
00090   belosLOWSFPL_gmres.set("Output Frequency",int(outputFrequency));
00091   belosLOWSFPL_gmres.set("Show Maximum Residual Norm Only",bool(outputMaxResOnly));
00092 
00093 #ifdef HAVE_BELOS_IFPACK
00094   //
00095   // Set the parameters for the Ifpack Preconditioner Factory and create parameter list
00096   //
00097   Teuchos::ParameterList &ifpackPFSL = belosLOWSFPL->sublist("IfpackPreconditionerFactory");
00098   
00099   ifpackPFSL.set("Overlap",int(2));
00100   ifpackPFSL.set("Prec Type","ILUT");
00101 #endif
00102 
00103   // Whether the linear solver succeeded.
00104   // (this will be set during the residual check at the end)
00105   bool success = true;
00106 
00107   // Number of random right-hand sides we will be solving for.
00108   int numRhs = 5;
00109 
00110   // Name of input matrix file
00111   std::string matrixFile = "orsirr1.hb";
00112 
00113   // Read in the matrix file (can be *.mtx, *.hb, etc.)
00114   Epetra_SerialComm comm;
00115   Teuchos::RCP<Epetra_CrsMatrix> epetra_A;
00116   EpetraExt::readEpetraLinearSystem( matrixFile, comm, &epetra_A );
00117   
00118   // Create a Thyra linear operator (A) using the Epetra_CrsMatrix (epetra_A).
00119   Teuchos::RCP<const Thyra::LinearOpBase<double> > 
00120     A = Thyra::epetraLinearOp(epetra_A);
00121 
00122   // Get the domain space for the Thyra linear operator 
00123   Teuchos::RCP<const Thyra::VectorSpaceBase<double> > domain = A->domain();
00124 
00125   // Create the Belos LOWS factory.
00126   Teuchos::RCP<Thyra::LinearOpWithSolveFactoryBase<double> >
00127     belosLOWSFactory = Teuchos::rcp(new Thyra::BelosLinearOpWithSolveFactory<double>());
00128 
00129 #ifdef HAVE_BELOS_IFPACK
00130 
00131   // Set the preconditioner factory for the LOWS factory.
00132   belosLOWSFactory->setPreconditionerFactory(
00133                Teuchos::rcp(new Thyra::IfpackPreconditionerFactory())
00134                ,"IfpackPreconditionerFactory"
00135                );
00136 #endif
00137 
00138   // Set the parameter list to specify the behavior of the factory.
00139   belosLOWSFactory->setParameterList( belosLOWSFPL );
00140 
00141   // Set the output stream and the verbosity level (prints to std::cout by defualt)
00142   belosLOWSFactory->setVerbLevel(Teuchos::VERB_LOW);
00143 
00144   // Create a BelosLinearOpWithSolve object from the Belos LOWS factory.
00145   Teuchos::RCP<Thyra::LinearOpWithSolveBase<double> >
00146     nsA = belosLOWSFactory->createOp();
00147 
00148   // Initialize the BelosLinearOpWithSolve object with the Thyra linear operator.
00149   Thyra::initializeOp<double>( *belosLOWSFactory, A, &*nsA );
00150 
00151   // Create a right-hand side with numRhs vectors in it and randomize it.
00152   Teuchos::RCP< Thyra::MultiVectorBase<double> > 
00153     b = Thyra::createMembers(domain, numRhs);
00154   Thyra::seed_randomize<double>(0);
00155   Thyra::randomize(-1.0, 1.0, &*b);
00156 
00157   // Create an initial std::vector with numRhs vectors in it and initialize it to zero.
00158   Teuchos::RCP< Thyra::MultiVectorBase<double> >
00159     x = Thyra::createMembers(domain, numRhs);
00160   Thyra::assign(&*x, 0.0);
00161 
00162   // Perform solve using the linear operator to get the approximate solution of Ax=b,
00163   // where b is the right-hand side and x is the left-hand side.
00164   Thyra::SolveStatus<double> solveStatus;
00165   solveStatus = Thyra::solve( *nsA, Thyra::NONCONJ_ELE, *b, &*x );
00166 
00167   // Print out status of solve.
00168   *out << "\nBelos LOWS Status: "<< solveStatus << std::endl;
00169 
00170   //
00171   // Compute residual and double check convergence.
00172   //
00173   std::vector<double> norm_b(numRhs), norm_res(numRhs);
00174   Teuchos::RCP< Thyra::MultiVectorBase<double> >
00175     y = Thyra::createMembers(domain, numRhs);
00176 
00177   // Compute the column norms of the right-hand side b.
00178   Thyra::norms_2( *b, &norm_b[0] );
00179 
00180   // Compute y=A*x, where x is the solution from the linear solver.
00181   A->apply( Thyra::NONCONJ_ELE, *x, &*y );
00182   
00183   // Compute A*x-b = y-b
00184   Thyra::update( -1.0, *b, &*y );
00185 
00186   // Compute the column norms of A*x-b.
00187   Thyra::norms_2( *y, &norm_res[0] );
00188 
00189   // Print out the final relative residual norms.
00190   double rel_res = 0.0;
00191   *out << "Final relative residual norms" << std::endl;  
00192   for (int i=0; i<numRhs; ++i) {
00193     rel_res = norm_res[i]/norm_b[i];
00194     if (rel_res > maxResid)
00195       success = false;
00196     *out << "RHS " << i+1 << " : " 
00197          << std::setw(16) << std::right << rel_res << std::endl;
00198   }
00199 
00200   return ( success ? 0 : 1 );
00201 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines