test_bl_cg_complex_hb.cpp

Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 //
00004 //                 Belos: Block Linear Solvers Package
00005 //                 Copyright (2004) Sandia Corporation
00006 //
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 //
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00025 //
00026 // ***********************************************************************
00027 // @HEADER
00028 //
00029 // This driver reads a problem from a Harwell-Boeing (HB) file.
00030 // The right-hand-side from the HB file is used instead of random vectors.
00031 // The initial guesses are all set to zero. 
00032 //
00033 // NOTE: No preconditioner is used in this case. 
00034 //
00035 #include "BelosConfigDefs.hpp"
00036 #include "BelosLinearProblem.hpp"
00037 #include "BelosBlockCGSolMgr.hpp"
00038 #include "Teuchos_CommandLineProcessor.hpp"
00039 #include "Teuchos_ParameterList.hpp"
00040 
00041 #ifdef HAVE_MPI
00042 #include <mpi.h>
00043 #endif
00044 
00045 // I/O for Harwell-Boeing files
00046 #ifdef HAVE_BELOS_TRIUTILS
00047 #include "iohb.h"
00048 #endif
00049 
00050 #include "MyMultiVec.hpp"
00051 #include "MyBetterOperator.hpp"
00052 #include "MyOperator.hpp"
00053 
00054 namespace Belos {
00055   class MPIFinalize {
00056   public:
00057     ~MPIFinalize() {
00058 #ifdef HAVE_MPI
00059       MPI_Finalize();
00060 #endif
00061     }
00062   };
00063 }
00064 
00065 using namespace Teuchos;
00066 
00067 int main(int argc, char *argv[]) {
00068   //
00069 #ifdef HAVE_COMPLEX
00070   typedef std::complex<double> ST;
00071 #elif HAVE_COMPLEX_H
00072   typedef std::complex<double> ST;
00073 #else
00074   std::cout << "Not compiled with std::complex support." << std::endl;
00075   std::cout << "End Result: TEST FAILED" << std::endl;
00076   return -1;
00077 #endif
00078 
00079   typedef ScalarTraits<ST>                 SCT;
00080   typedef SCT::magnitudeType                MT;
00081   typedef Belos::MultiVec<ST>               MV;
00082   typedef Belos::Operator<ST>               OP;
00083   typedef Belos::MultiVecTraits<ST,MV>     MVT;
00084   typedef Belos::OperatorTraits<ST,MV,OP>  OPT;
00085   ST one  = SCT::one();
00086   ST zero = SCT::zero();  
00087 
00088   int info = 0;
00089   int MyPID = 0;
00090   bool norm_failure = false;
00091 
00092 #ifdef HAVE_MPI 
00093   // Initialize MPI
00094   MPI_Init(&argc,&argv);
00095   Belos::MPIFinalize mpiFinalize; // Will call finalize with *any* return
00096   (void)mpiFinalize;
00097 #endif
00098   //
00099   using Teuchos::RCP;
00100   using Teuchos::rcp;
00101 
00102   bool verbose = false, proc_verbose = false;
00103   int frequency = -1;  // how often residuals are printed by solver
00104   int blocksize = 1;
00105   int numrhs = 1;
00106   int maxrestarts = 15;
00107   int length = 50;
00108   std::string filename("mhd1280b.cua");
00109   MT tol = 1.0e-5;  // relative residual tolerance
00110 
00111   CommandLineProcessor cmdp(false,true);
00112   cmdp.setOption("verbose","quiet",&verbose,"Print messages and results.");
00113   cmdp.setOption("frequency",&frequency,"Solvers frequency for printing residuals (#iters).");
00114   cmdp.setOption("filename",&filename,"Filename for Harwell-Boeing test matrix.");
00115   cmdp.setOption("tol",&tol,"Relative residual tolerance used by CG solver.");
00116   cmdp.setOption("num-rhs",&numrhs,"Number of right-hand sides to be solved for.");
00117   cmdp.setOption("num-restarts",&maxrestarts,"Maximum number of restarts allowed for the CG solver.");
00118   cmdp.setOption("blocksize",&blocksize,"Block size used by CG .");
00119   cmdp.setOption("subspace-length",&length,"Maximum dimension of block-subspace used by CG solver.");
00120   if (cmdp.parse(argc,argv) != CommandLineProcessor::PARSE_SUCCESSFUL) {
00121     return -1;
00122   }
00123   
00124   proc_verbose = verbose && (MyPID==0);  /* Only print on the zero processor */
00125   if (proc_verbose) {
00126     std::cout << Belos::Belos_Version() << std::endl << std::endl;
00127   }
00128   if (!verbose)
00129     frequency = -1;  // reset frequency if test is not verbose
00130   
00131   
00132 #ifndef HAVE_BELOS_TRIUTILS
00133   std::cout << "This test requires Triutils. Please configure with --enable-triutils." << std::endl;
00134   if (MyPID==0) {
00135     std::cout << "End Result: TEST FAILED" << std::endl;  
00136   }
00137   return -1;
00138 #endif
00139   
00140   // Get the data from the HB file
00141   int dim,dim2,nnz;
00142   MT *dvals;
00143   int *colptr,*rowind;
00144   ST *cvals;
00145   nnz = -1;
00146   info = readHB_newmat_double(filename.c_str(),&dim,&dim2,&nnz,
00147                               &colptr,&rowind,&dvals);
00148   if (info == 0 || nnz < 0) {
00149     if (MyPID==0) {
00150       std::cout << "Error reading '" << filename << "'" << std::endl;
00151       std::cout << "End Result: TEST FAILED" << std::endl;
00152     }
00153     return -1;
00154   }
00155   // Convert interleaved doubles to std::complex values
00156   cvals = new ST[nnz];
00157   for (int ii=0; ii<nnz; ii++) {
00158     cvals[ii] = ST(dvals[ii*2],dvals[ii*2+1]);
00159   }
00160   // Build the problem matrix
00161   RCP< MyBetterOperator<ST> > A 
00162     = rcp( new MyBetterOperator<ST>(dim,colptr,nnz,rowind,cvals) );
00163   //
00164   // ********Other information used by block solver***********
00165   // *****************(can be user specified)******************
00166   //
00167   int maxits = dim/blocksize; // maximum number of iterations to run
00168   //
00169   ParameterList belosList;
00170   belosList.set( "Num Blocks", length );                 // Maximum number of blocks in Krylov factorization
00171   belosList.set( "Block Size", blocksize );              // Blocksize to be used by iterative solver
00172   belosList.set( "Maximum Iterations", maxits );         // Maximum number of iterations allowed
00173   belosList.set( "Maximum Restarts", maxrestarts );      // Maximum number of restarts allowed
00174   belosList.set( "Convergence Tolerance", tol );         // Relative convergence tolerance requested
00175   if (verbose) {
00176     belosList.set( "Verbosity", Belos::Errors + Belos::Warnings + 
00177        Belos::TimingDetails + Belos::FinalSummary + Belos::StatusTestDetails );
00178     if (frequency > 0)
00179       belosList.set( "Output Frequency", frequency );
00180   }
00181   else
00182     belosList.set( "Verbosity", Belos::Errors + Belos::Warnings );
00183   //
00184   // Construct the right-hand side and solution multivectors.
00185   // NOTE:  The right-hand side will be constructed such that the solution is
00186   // a vectors of one.
00187   //
00188   RCP<MyMultiVec<ST> > soln = rcp( new MyMultiVec<ST>(dim,numrhs) );
00189   RCP<MyMultiVec<ST> > rhs = rcp( new MyMultiVec<ST>(dim,numrhs) );
00190   MVT::MvRandom( *soln );
00191   OPT::Apply( *A, *soln, *rhs );
00192   MVT::MvInit( *soln, zero );
00193   //
00194   //  Construct an unpreconditioned linear problem instance.
00195   //
00196   RCP<Belos::LinearProblem<ST,MV,OP> > problem = 
00197     rcp( new Belos::LinearProblem<ST,MV,OP>( A, soln, rhs ) );
00198   bool set = problem->setProblem();
00199   if (set == false) {
00200     if (proc_verbose)
00201       std::cout << std::endl << "ERROR:  Belos::LinearProblem failed to set up correctly!" << std::endl;
00202     return -1;
00203   }
00204   //
00205   // *******************************************************************
00206   // *************Start the block CG iteration***********************
00207   // *******************************************************************
00208   //
00209   Belos::BlockCGSolMgr<ST,MV,OP> solver( problem, rcp(&belosList,false) );
00210 
00211   //
00212   // **********Print out information about problem*******************
00213   //
00214   if (proc_verbose) {
00215     std::cout << std::endl << std::endl;
00216     std::cout << "Dimension of matrix: " << dim << std::endl;
00217     std::cout << "Number of right-hand sides: " << numrhs << std::endl;
00218     std::cout << "Block size used by solver: " << blocksize << std::endl;
00219     std::cout << "Max number of CG iterations: " << maxits << std::endl; 
00220     std::cout << "Relative residual tolerance: " << tol << std::endl;
00221     std::cout << std::endl;
00222   }
00223   //
00224   // Perform solve
00225   //
00226   Belos::ReturnType ret = solver.solve();
00227   //
00228   // Compute actual residuals.
00229   //
00230   RCP<MyMultiVec<ST> > temp = rcp( new MyMultiVec<ST>(dim,numrhs) );
00231   OPT::Apply( *A, *soln, *temp );
00232   MVT::MvAddMv( one, *rhs, -one, *temp, *temp );
00233   std::vector<MT> norm_num(numrhs), norm_denom(numrhs);
00234   MVT::MvNorm( *temp, norm_num );
00235   MVT::MvNorm( *rhs, norm_denom );
00236   for (int i=0; i<numrhs; ++i) {
00237     if (proc_verbose) 
00238       std::cout << "Relative residual "<<i<<" : " << norm_num[i] / norm_denom[i] << std::endl;
00239     if ( norm_num[i] / norm_denom[i] > tol ) {
00240       norm_failure = true;
00241     }
00242   }
00243   
00244   // Clean up.
00245   delete [] dvals;
00246   delete [] colptr;
00247   delete [] rowind;
00248   delete [] cvals;
00249 
00250   if ( ret!=Belos::Converged || norm_failure ) {
00251     if (proc_verbose)
00252       std::cout << "End Result: TEST FAILED" << std::endl;  
00253     return -1;
00254   }
00255   //
00256   // Default return value
00257   //
00258   if (proc_verbose)
00259     std::cout << "End Result: TEST PASSED" << std::endl;
00260   return 0;
00261 
00262   //
00263 } // end test_bl_gmres_complex_hb.cpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
Generated on Wed Apr 13 10:05:15 2011 for Belos Package Browser (Single Doxygen Collection) by  doxygen 1.6.3