belos_tpetra_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 // Tpetra::CisMatrix.  This matrix is then converted into a Thyra linear operator
00032 // through the Thyra-Tpetra adapters.
00033 //
00034 // The right-hand-side from the HB file is used instead of 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 // I/O for Harwell-Boeing files
00050 #ifdef HAVE_BELOS_TRIUTILS
00051 #include "iohb.h"
00052 #endif
00053 
00054 // Teuchos includes
00055 #include "Teuchos_ParameterList.hpp"
00056 #include "Teuchos_VerboseObject.hpp"
00057 #include "Teuchos_GlobalMPISession.hpp"
00058 
00059 #ifdef HAVE_MPI
00060 #  include "Tpetra_MpiPlatform.hpp"
00061 #else
00062 #  include "Tpetra_SerialPlatform.hpp"
00063 #endif
00064 
00065 // My Tpetra::Operator include
00066 #include "MyOperator.hpp"
00067 
00068 int main(int argc, char* argv[])
00069 {
00070   //
00071   // Get a default output stream from the Teuchos::VerboseObjectBase
00072   //
00073   Teuchos::RCP<Teuchos::FancyOStream>
00074     out = Teuchos::VerboseObjectBase::getDefaultOStream();
00075   
00076   Teuchos::GlobalMPISession mpiSession(&argc,&argv);
00077 
00078 #ifdef HAVE_COMPLEX
00079   typedef std::complex<double> ST;  // Scalar-type typedef
00080 #elif HAVE_COMPLEX_H
00081   typedef std::complex<double> ST;     // Scalar-type typedef
00082 #else
00083   typedef double ST;                // Scalar-type typedef
00084 #endif
00085   
00086   typedef Teuchos::ScalarTraits<ST>::magnitudeType MT;  // Magnitude-type typedef
00087   typedef int OT;                   // Ordinal-type typedef
00088   ST one = Teuchos::ScalarTraits<ST>::one(); 
00089   ST zero = Teuchos::ScalarTraits<ST>::zero(); 
00090   
00091 #ifdef HAVE_MPI
00092   MPI_Comm mpiComm = MPI_COMM_WORLD;
00093   const Tpetra::MpiPlatform<OT,OT>  ordinalPlatform(mpiComm);
00094   const Tpetra::MpiPlatform<OT,ST>   scalarPlatform(mpiComm);
00095 #else
00096   const Tpetra::SerialPlatform<OT,OT>  ordinalPlatform;
00097   const Tpetra::SerialPlatform<OT,ST>   scalarPlatform;
00098 #endif
00099   
00100   //
00101   // Get the data from the HB file
00102   //
00103   
00104   // Name of input matrix file
00105   std::string matrixFile = "mhd1280b.cua";
00106   
00107   int info=0;
00108   int dim,dim2,nnz;
00109   MT *dvals;
00110   int *colptr,*rowind;
00111   ST *cvals;
00112   nnz = -1;
00113   info = readHB_newmat_double(matrixFile.c_str(),&dim,&dim2,&nnz,
00114                               &colptr,&rowind,&dvals);
00115 
00116   if (info == 0 || nnz < 0) {
00117     *out << "Error reading '" << matrixFile << "'" << std::endl;
00118   }
00119 #ifdef HAVE_MPI
00120   MPI_Finalize();
00121 #endif
00122 
00123   // Convert interleaved doubles to std::complex values
00124   cvals = new ST[nnz];
00125   for (int ii=0; ii<nnz; ii++) {
00126     cvals[ii] = ST(dvals[ii*2],dvals[ii*2+1]);
00127   }
00128   
00129   // Declare global dimension of the linear operator
00130   OT globalDim = dim;
00131   
00132   // Create the element space and std::vector space
00133   const Tpetra::ElementSpace<OT> elementSpace(globalDim,0,ordinalPlatform);
00134   const Tpetra::VectorSpace<OT,ST> vectorSpace(elementSpace,scalarPlatform);
00135   
00136   // Create my implementation of a Tpetra::Operator
00137   RCP<Tpetra::Operator<OT,ST> >
00138     tpetra_A = rcp( new MyOperator<OT,ST>(vectorSpace,dim,colptr,nnz,rowind,cvals) );
00139 
00140   // Create a Thyra linear operator (A) using the Tpetra::CisMatrix (tpetra_A).
00141   RCP<Thyra::LinearOpBase<ST> >
00142     A = Teuchos::rcp( new Thyra::TpetraLinearOp<OT,ST>(tpetra_A) );
00143 
00144   //
00145   // Set the parameters for the Belos LOWS Factory and create a parameter list.
00146   //
00147   int             blockSize              = 1;
00148   int             maxIterations          = globalDim;
00149   int             maxRestarts            = 15;
00150   int             gmresKrylovLength      = 50;
00151   int             outputFrequency        = 100;
00152   bool            outputMaxResOnly       = true;
00153   MT              maxResid               = 1e-5;
00154 
00155   Teuchos::RCP<Teuchos::ParameterList>
00156     belosLOWSFPL = Teuchos::rcp( new Teuchos::ParameterList() );
00157  
00158   belosLOWSFPL->set("Solver Type","Block GMRES");
00159 
00160   Teuchos::ParameterList& belosLOWSFPL_solver =
00161     belosLOWSFPL->sublist("Solver Types");
00162 
00163   Teuchos::ParameterList& belosLOWSFPL_gmres =
00164     belosLOWSFPL_solver.sublist("Block GMRES");
00165 
00166   belosLOWSFPL_gmres.set("Maximum Iterations",int(maxIterations));
00167   belosLOWSFPL_gmres.set("Convergence Tolerance",MT(maxResid));
00168   belosLOWSFPL_gmres.set("Maximum Restarts",int(maxRestarts));
00169   belosLOWSFPL_gmres.set("Block Size",int(blockSize));
00170   belosLOWSFPL_gmres.set("Num Blocks",int(gmresKrylovLength));
00171   belosLOWSFPL_gmres.set("Output Frequency",int(outputFrequency));
00172   belosLOWSFPL_gmres.set("Show Maximum Residual Norm Only",bool(outputMaxResOnly));
00173  
00174   // Whether the linear solver succeeded.
00175   // (this will be set during the residual check at the end)
00176   bool success = true;
00177 
00178   // Number of random right-hand sides we will be solving for.
00179   int numRhs = 1;
00180 
00181   // Get the domain space for the Thyra linear operator 
00182   Teuchos::RCP<const Thyra::VectorSpaceBase<ST> > domain = A->domain();
00183 
00184   // Create the Belos LOWS factory.
00185   Teuchos::RCP<Thyra::LinearOpWithSolveFactoryBase<ST> >
00186     belosLOWSFactory = Teuchos::rcp(new Thyra::BelosLinearOpWithSolveFactory<ST>());
00187 
00188   // Set the parameter list to specify the behavior of the factory.
00189   belosLOWSFactory->setParameterList( belosLOWSFPL );
00190 
00191   // Set the output stream and the verbosity level (prints to std::cout by defualt)
00192   // NOTE:  Set to VERB_NONE for no output from the solver.
00193   belosLOWSFactory->setVerbLevel(Teuchos::VERB_LOW);
00194 
00195   // Create a BelosLinearOpWithSolve object from the Belos LOWS factory.
00196   Teuchos::RCP<Thyra::LinearOpWithSolveBase<ST> >
00197     nsA = belosLOWSFactory->createOp();
00198 
00199   // Initialize the BelosLinearOpWithSolve object with the Thyra linear operator.
00200   Thyra::initializeOp<ST>( *belosLOWSFactory, A, &*nsA );
00201 
00202   // Create a right-hand side with numRhs vectors in it.
00203   Teuchos::RCP< Thyra::MultiVectorBase<ST> > 
00204     b = Thyra::createMembers(domain, numRhs);
00205 
00206   // Create an initial std::vector with numRhs vectors in it and initialize it to one.
00207   Teuchos::RCP< Thyra::MultiVectorBase<ST> >
00208     x = Thyra::createMembers(domain, numRhs);
00209   Thyra::assign(&*x, one);
00210 
00211   // Initialize the right-hand side so that the solution is a std::vector of ones.
00212   A->apply( Thyra::NONCONJ_ELE, *x, &*b );
00213   Thyra::assign(&*x, zero);
00214 
00215   // Perform solve using the linear operator to get the approximate solution of Ax=b,
00216   // where b is the right-hand side and x is the left-hand side.
00217   Thyra::SolveStatus<ST> solveStatus;
00218   solveStatus = Thyra::solve( *nsA, Thyra::NONCONJ_ELE, *b, &*x );
00219 
00220   // Print out status of solve.
00221   *out << "\nBelos LOWS Status: "<< solveStatus << std::endl;
00222 
00223   //
00224   // Compute residual and ST check convergence.
00225   //
00226   std::vector<MT> norm_b(numRhs), norm_res(numRhs);
00227   Teuchos::RCP< Thyra::MultiVectorBase<ST> >
00228     y = Thyra::createMembers(domain, numRhs);
00229 
00230   // Compute the column norms of the right-hand side b.
00231   Thyra::norms_2( *b, &norm_b[0] );
00232 
00233   // Compute y=A*x, where x is the solution from the linear solver.
00234   A->apply( Thyra::NONCONJ_ELE, *x, &*y );
00235   
00236   // Compute A*x-b = y-b
00237   Thyra::update( -one, *b, &*y );
00238 
00239   // Compute the column norms of A*x-b.
00240   Thyra::norms_2( *y, &norm_res[0] );
00241 
00242   // Print out the final relative residual norms.
00243   MT rel_res = 0.0;
00244   *out << "Final relative residual norms" << std::endl;  
00245   for (int i=0; i<numRhs; ++i) {
00246     rel_res = norm_res[i]/norm_b[i];
00247     if (rel_res > maxResid)
00248       success = false;
00249     *out << "RHS " << i+1 << " : " 
00250    << std::setw(16) << std::right << rel_res << std::endl;
00251   }
00252 
00253   return ( success ? 0 : 1 );
00254 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends
Generated on Wed Apr 13 10:19:51 2011 for Stratimikos Package Browser (Single Doxygen Collection) by  doxygen 1.6.3