Thyra Operator/Vector Support : Templated SPMD Implementation of the CG Method

template<class Scalar>
class ExampleTridiagSpmdLinearOp : public Thyra::SpmdLinearOpBase<Scalar> {
private:

  Teuchos::RefCountPtr<const Teuchos::Comm<Thyra::Index> >  comm_;
  int                  procRank_;
  int                  numProcs_;
  Thyra::Index         localDim_;
  std::vector<Scalar>  lower_;   // size = ( procRank == 0         ? localDim - 1 : localDim )    
  std::vector<Scalar>  diag_;    // size = localDim
  std::vector<Scalar>  upper_;   // size = ( procRank == numProc-1 ? localDim - 1 : localDim )

  void communicate( const bool first, const bool last, const Scalar x[], Scalar *x_km1, Scalar *x_kp1 ) const;

public:

  using Thyra::SpmdLinearOpBase<Scalar>::euclideanApply;

  ExampleTridiagSpmdLinearOp() : procRank_(0), numProcs_(0) {}

  ExampleTridiagSpmdLinearOp(
    const Teuchos::RefCountPtr<const Teuchos::Comm<Thyra::Index> > &comm
    ,const Thyra::Index localDim, const Scalar lower[], const Scalar diag[], const Scalar upper[] )
    { this->initialize(comm,localDim,lower,diag,upper); }
  
  void initialize(
    const Teuchos::RefCountPtr<const Teuchos::Comm<Thyra::Index> > &comm
    ,const Thyra::Index             localDim // >= 2
    ,const Scalar                   lower[]  // size == ( procRank == 0         ? localDim - 1 : localDim )
    ,const Scalar                   diag[]   // size == localDim
    ,const Scalar                   upper[]  // size == ( procRank == numProc-1 ? localDim - 1 : localDim )
    )
    {
      TEST_FOR_EXCEPT( localDim < 2 );
      this->setLocalDimensions(comm,localDim,localDim); // Needed to set up range() and domain()
      comm_  = Teuchos::DefaultComm<Thyra::Index>::getDefaultSerialComm(comm);
      localDim_ = localDim;
      numProcs_ = comm->getSize();
      procRank_ = comm->getRank();
      const Thyra::Index
        lowerDim = ( procRank_ == 0           ? localDim - 1 : localDim ),
        upperDim = ( procRank_ == numProcs_-1 ? localDim - 1 : localDim );
      lower_.resize(lowerDim);  for( int k = 0; k < lowerDim; ++k ) lower_[k] = lower[k];
      diag_.resize(localDim);   for( int k = 0; k < localDim; ++k ) diag_[k]  = diag[k];
      upper_.resize(upperDim);  for( int k = 0; k < upperDim; ++k ) upper_[k] = upper[k];
    }

  // Overridden form Teuchos::Describable */

  std::string description() const
    {
      return (std::string("ExampleTridiagSpmdLinearOp<") + Teuchos::ScalarTraits<Scalar>::name() + std::string(">"));
    }

protected:


  // Overridden from SingleScalarEuclideanLinearOpBase

  bool opSupported( Thyra::ETransp M_trans ) const
    {
      typedef Teuchos::ScalarTraits<Scalar> ST;
      return (M_trans == Thyra::NOTRANS || (!ST::isComplex && M_trans == Thyra::CONJ) );
    }

  // Overridden from SerialLinearOpBase

  void euclideanApply(
    const Thyra::ETransp                         M_trans
    ,const RTOpPack::ConstSubVectorView<Scalar>  &local_x_in
    ,const RTOpPack::SubVectorView<Scalar>       *local_y_out
    ,const Scalar                                alpha
    ,const Scalar                                beta
    ) const
    {
      typedef Teuchos::ScalarTraits<Scalar> ST;
      TEST_FOR_EXCEPTION( M_trans != Thyra::NOTRANS, std::logic_error, "Error, can not handle transpose!" );
      // Get constants
      const Scalar zero = ST::zero();
      // Get raw pointers to vector data to make me feel better!
      const Scalar *x = local_x_in.values();
      Scalar       *y = local_y_out->values();
      // Determine what process we are
      const bool first = ( procRank_ == 0 ), last = ( procRank_ == numProcs_-1 );
      // Communicate ghost elements
      Scalar x_km1, x_kp1;
      communicate( first, last, x, &x_km1, &x_kp1 );
      // Perform operation (if beta==0 then we must be careful since y could be uninitialized on input!)
      Thyra::Index k = 0, lk = 0;
      if( beta == zero ) {
        y[k] = alpha * ( (first?zero:lower_[lk]*x_km1) + diag_[k]*x[k] + upper_[k]*x[k+1] ); if(!first) ++lk;             // First local row
        for( k = 1; k < localDim_ - 1; ++lk, ++k )
          y[k] = alpha * ( lower_[lk]*x[k-1] + diag_[k]*x[k] + upper_[k]*x[k+1] );                                        // Middle local rows
        y[k] = alpha * ( lower_[lk]*x[k-1] + diag_[k]*x[k] + (last?zero:upper_[k]*x_kp1) );                               // Last local row
      }
      else {
        y[k] = alpha * ( (first?zero:lower_[lk]*x_km1) + diag_[k]*x[k] + upper_[k]*x[k+1] ) + beta*y[k]; if(!first) ++lk; // First local row
        for( k = 1; k < localDim_ - 1; ++lk, ++k )
          y[k] = alpha * ( lower_[lk]*x[k-1] + diag_[k]*x[k] + upper_[k]*x[k+1] ) + beta*y[k];                            // Middle local rows
        y[k] = alpha * ( lower_[lk]*x[k-1] + diag_[k]*x[k] + (last?zero:upper_[k]*x_kp1) ) + beta*y[k];                   // Last local row
      }
      //std::cout << "\ny = ["; for(k=0;k<localDim_;++k) { std::cout << y[k]; if(k<localDim_-1) std::cout << ","; } std::cout << "]\n";
    }

};  // end class ExampleTridiagSpmdLinearOp

The above SPMD matrix class is used in an example program (see runCgSolveExample() below) that calls sillyCgSolve() or silliestCgSolve(). In this example program, the matrix constructed and used is the well-known tridiagonal matrix

$A= \left[\begin{array}{rrrrrrrrrr} 2 a & -1 \\ -1 & 2 a & -1 \\ & \ddots & \ddots & \ddots \\ & & -1 & 2 a & -1 \\ & & & -1 & 2 a \end{array}\right]$

where

is an adjustable diagonal scale factor that makes the matrix $A$

more or less well conditioned.

The CG method is then run on the matrix $A$

for a number of iterations or until convergence to some tolerance is achieved.

The following templated function runCgSolveExample() implements the example described above:

template<class Scalar>
bool runCgSolveExample(
  const Teuchos::RefCountPtr<const Teuchos::Comm<Thyra::Index> > &comm
  ,const int                                                     procRank
  ,const int                                                     numProc
  ,const int                                                     localDim
  ,const Scalar                                                  diagScale
  ,const bool                                                    showAllTests
  ,const bool                                                    verbose
  ,const bool                                                    dumpAll
  ,const typename Teuchos::ScalarTraits<Scalar>::magnitudeType   tolerance
  ,const int                                                     maxNumIters
  )
{
  using Teuchos::RefCountPtr;
  using Teuchos::rcp;
  using Teuchos::OSTab;
  typedef Teuchos::ScalarTraits<Scalar> ST;
  typedef typename ST::magnitudeType    ScalarMag;
  bool success = true;
  bool result;
  // ToDo: Get VerboseObjectBase to automatically setup for parallel
  Teuchos::RefCountPtr<Teuchos::FancyOStream>
    out = (verbose ? Teuchos::VerboseObjectBase::getDefaultOStream() : Teuchos::null);
  if(verbose)
    *out << "\n***\n*** Running silly CG solver using scalar type = \'" << ST::name() << "\' ...\n***\n";
  Teuchos::Time timer("");
  timer.start(true);
  //
  // (A) Setup a simple linear system with tridiagonal operator:
  //
  //       [  a*2   -1                ]
  //       [ -1    a*2  -1            ]
  //  A =  [         .   .    .       ]
  //       [            -1  a*2    -1 ]
  //       [                 -1   a*2 ]
  //
  // (A.1) Create the tridiagonal matrix operator
  if(verbose)
    *out << "\nConstructing tridiagonal matrix A of local dimension = " << localDim
         << " and diagonal multiplier = " << diagScale << " ...\n";
  const Thyra::Index
    lowerDim = ( procRank == 0         ? localDim - 1 : localDim ),
    upperDim = ( procRank == numProc-1 ? localDim - 1 : localDim );
  std::vector<Scalar> lower(lowerDim), diag(localDim), upper(upperDim);
  const Scalar one = ST::one(), diagTerm = Scalar(2)*diagScale*ST::one();
  int k = 0, kl = 0;
  if(procRank > 0) { lower[kl] = -one; ++kl; };  diag[k] = diagTerm; upper[k] = -one; // First local row
  for( k = 1; k < localDim - 1; ++k, ++kl ) {
    lower[kl] = -one; diag[k] = diagTerm; upper[k] = -one;                            // Middle local rows
  }
  lower[kl] = -one; diag[k] = diagTerm; if(procRank < numProc-1) upper[k] = -one;     // Last local row
  RefCountPtr<const Thyra::LinearOpBase<Scalar> >
    A = rcp(new ExampleTridiagSpmdLinearOp<Scalar>(comm,localDim,&lower[0],&diag[0],&upper[0]));
  if(verbose) *out << "\nGlobal dimension of A = " << A->domain()->dim() << std::endl;
  // (A.2) Testing the linear operator constructed linear operator
  if(verbose) *out << "\nTesting the constructed linear operator A ...\n";
  Thyra::LinearOpTester<Scalar> linearOpTester;
  linearOpTester.dump_all(dumpAll);
  linearOpTester.set_all_error_tol(tolerance);
  linearOpTester.set_all_warning_tol(ScalarMag(ScalarMag(1e-2)*tolerance));
  linearOpTester.show_all_tests(true);
  linearOpTester.check_adjoint(false);
  linearOpTester.check_for_symmetry(true);
  linearOpTester.show_all_tests(showAllTests);
  result = linearOpTester.check(*A,out.get());
  if(!result) success = false;
  // (A.3) Create RHS vector b and set to a random value
  RefCountPtr<Thyra::VectorBase<Scalar> > b = createMember(A->range());
  Thyra::seed_randomize<Scalar>(0);
  Thyra::randomize( Scalar(-ST::one()), Scalar(+ST::one()), &*b );
  // (A.4) Create LHS vector x and set to zero
  RefCountPtr<Thyra::VectorBase<Scalar> > x = createMember(A->domain());
  Thyra::assign( &*x, ST::zero() );
  //
  // (B) Solve the linear system with the silly CG solver
  //
  if(verbose) *out << "\nSolving the linear system with sillyCgSolve(...) ...\n";
  result = sillyCgSolve(*A,*b,maxNumIters,tolerance,&*x,OSTab(out).getOStream().get());
  if(!result) success = false;
  //
  // (C) Check that the linear system was solved to the specified tolerance
  //
  RefCountPtr<Thyra::VectorBase<Scalar> > r = createMember(A->range());                     
  Thyra::assign(&*r,*b);                                               // r = b
  Thyra::apply(*A,Thyra::NOTRANS,*x,&*r,Scalar(-ST::one()),ST::one()); // r = -A*x + r
  const ScalarMag r_nrm = Thyra::norm(*r), b_nrm = Thyra::norm(*b);
  const ScalarMag rel_err = r_nrm/b_nrm, relaxTol = ScalarMag(10.0)*tolerance;
  result = rel_err <= relaxTol;
  if(!result) success = false;
  if(verbose) {
    *out << "\nChecking the residual ourselves ...\n";
    if(1){
      OSTab tab(out);
      *out
        << "\n||b-A*x||/||b|| = "<<r_nrm<<"/"<<b_nrm<<" = "<<rel_err<<(result?" <= ":" > ")
        <<"10.0*tolerance = "<<relaxTol<<": "<<(result?"passed":"failed")<<std::endl;
    }
  }
  timer.stop();
  if(verbose) *out << "\nTotal time = " << timer.totalElapsedTime() << " sec\n";

  return success;
} // end runCgSolveExample()

The above templated function runCgSolveExample() is then instantiated with the following scalar types:

and is called multiple times from within the following main() program function:

int main(int argc, char *argv[])
{

  using Teuchos::CommandLineProcessor;

  bool success = true;
  bool verbose = true;
  bool result;

  Teuchos::GlobalMPISession mpiSession(&argc,&argv);
  const int procRank = Teuchos::GlobalMPISession::getRank();
  const int numProc = Teuchos::GlobalMPISession::getNProc();

  const Teuchos::RefCountPtr<const Teuchos::Comm<Thyra::Index> >
    comm = Teuchos::DefaultComm<Thyra::Index>::getComm();

  Teuchos::RefCountPtr<Teuchos::FancyOStream>
    out = Teuchos::VerboseObjectBase::getDefaultOStream();

  try {

    //
    // Read in command-line options
    //

    int    localDim    = 500;
    double diagScale   = 1.001;
    double tolerance   = 1e-4;
    bool   showAllTests = false;
    int    maxNumIters = 300;
    bool   dumpAll     = false;

    CommandLineProcessor  clp;
    clp.throwExceptions(false);
    clp.addOutputSetupOptions(true);

    clp.setOption( "verbose", "quiet", &verbose, "Determines if any output is printed or not." );
    clp.setOption( "local-dim", &localDim, "Local dimension of the linear system." );
    clp.setOption( "diag-scale", &diagScale, "Scaling of the diagonal to improve conditioning." );
    clp.setOption( "show-all-tests", "show-summary-only", &showAllTests, "Show all LinearOpTester tests or not" );
    clp.setOption( "tol", &tolerance, "Relative tolerance for linear system solve." );
    clp.setOption( "max-num-iters", &maxNumIters, "Maximum of CG iterations." );
    clp.setOption( "dump-all", "no-dump-all", &dumpAll, "Determines if vectors are printed or not." );

    CommandLineProcessor::EParseCommandLineReturn parse_return = clp.parse(argc,argv);
    if( parse_return != CommandLineProcessor::PARSE_SUCCESSFUL ) return parse_return;

    TEST_FOR_EXCEPTION( localDim < 2, std::logic_error, "Error, localDim=" << localDim << " < 2 is not allowed!" );

    // Run using float
    result = runCgSolveExample<float>(comm,procRank,numProc,localDim,diagScale,showAllTests,verbose,dumpAll,tolerance,maxNumIters);
    if(!result) success = false;

    // Run using double
    result = runCgSolveExample<double>(comm,procRank,numProc,localDim,diagScale,showAllTests,verbose,dumpAll,tolerance,maxNumIters);
    if(!result) success = false;

#if defined(HAVE_COMPLEX) && defined(HAVE_TEUCHOS_COMPLEX)

    // Run using std::complex<float>
    result = runCgSolveExample<std::complex<float> >(comm,procRank,numProc,localDim,diagScale,showAllTests,verbose,dumpAll,tolerance,maxNumIters);
    if(!result) success = false;

    // Run using std::complex<double>
    result = runCgSolveExample<std::complex<double> >(comm,procRank,numProc,localDim,diagScale,showAllTests,verbose,dumpAll,tolerance,maxNumIters);
    if(!result) success = false;

#endif    

  }
  TEUCHOS_STANDARD_CATCH_STATEMENTS(true,*out,success)

  if( verbose && procRank==0 ) {
    if(success) *out << "\nAll of the tests seem to have run successfully!\n";
    else        *out << "\nOh no! at least one of the tests failed!\n"; 
  }
  
  return success ? 0 : 1;

} // end main()

The above example program is built as part of the Thyra package (unless examples where disabled at configure time) and the executable can be found at:

 ./example/operator_vector/sillyCgSolve_mpi.exe

where ./ is the base build directory for Thyra (e.g. $TRILINOS_BUILD_DIR/packages/thyra).

This example program should run successfully with no arguments for any number of processors and for one processor produces the following output:

Teuchos::GlobalMPISession::GlobalMPISession(): started serial run

***
*** Running silly CG solver using scalar type = 'float' ...
***

Constructing tridiagonal matrix A of local dimension = 500 and diagonal multiplier = 1.001 ...

Global dimension of A = 500

Testing the constructed linear operator A ...
  
  *** Entering LinearOpTester<float,float>::check(op,...) ...
  
  describe op:
    ExampleTridiagSpmdLinearOp<float>{rangeDim=500,domainDim=500}
  
  Checking the domain and range spaces ... passed!
  
  this->check_linear_properties()==true: Checking the linear properties of the forward linear operator ... passed!
  
  (this->check_linear_properties()&&this->check_adjoint())==false: Skipping the check of the linear properties of the adjoint operator!
  
  this->check_adjoint()==false: Skipping check for the agreement of the adjoint and forward operators!
  
  this->check_for_symmetry()==true: Performing check of symmetry ... passed!
  
  Congratulations, this LinearOpBase object seems to check out!
  
  *** Leaving LinearOpTester<float,float>::check(...)

Solving the linear system with sillyCgSolve(...) ...
  
  Starting CG solver ...
  
  describe A:
    ExampleTridiagSpmdLinearOp<float>{rangeDim=500,domainDim=500}
  
  describe b:
    DefaultSpmdVector<float>{dim=500}
  
  describe x:
    DefaultSpmdVector<float>{dim=500}
  
  Iter = 0, ||b-A*x||/||b-A*x0|| = 1.000000e+00
  Iter = 31, ||b-A*x||/||b-A*x0|| = 3.172030e-01
  Iter = 62, ||b-A*x||/||b-A*x0|| = 4.875919e-02
  Iter = 93, ||b-A*x||/||b-A*x0|| = 1.090162e-02
  Iter = 124, ||b-A*x||/||b-A*x0|| = 2.477407e-03
  Iter = 155, ||b-A*x||/||b-A*x0|| = 8.997934e-04
  Iter = 186, ||b-A*x||/||b-A*x0|| = 1.621064e-04
  Iter = 202, ||b-A*x||/||b-A*x0|| = 9.563797e-05

Checking the residual ourselves ...
  
  ||b-A*x||/||b|| = 1.318619e-03/1.270762e+01 = 1.037660e-04 <= 10.0*tolerance = 9.999999e-04: passed

Total time = 7.228700e-02 sec

***
*** Running silly CG solver using scalar type = 'double' ...
***

Constructing tridiagonal matrix A of local dimension = 500 and diagonal multiplier = 1.001000e+00 ...

Global dimension of A = 500

Testing the constructed linear operator A ...
  
  *** Entering LinearOpTester<double,double>::check(op,...) ...
  
  describe op:
    ExampleTridiagSpmdLinearOp<double>{rangeDim=500,domainDim=500}
  
  Checking the domain and range spaces ... passed!
  
  this->check_linear_properties()==true: Checking the linear properties of the forward linear operator ... passed!
  
  (this->check_linear_properties()&&this->check_adjoint())==false: Skipping the check of the linear properties of the adjoint operator!
  
  this->check_adjoint()==false: Skipping check for the agreement of the adjoint and forward operators!
  
  this->check_for_symmetry()==true: Performing check of symmetry ... passed!
  
  Congratulations, this LinearOpBase object seems to check out!
  
  *** Leaving LinearOpTester<double,double>::check(...)

Solving the linear system with sillyCgSolve(...) ...
  
  Starting CG solver ...
  
  describe A:
    ExampleTridiagSpmdLinearOp<double>{rangeDim=500,domainDim=500}
  
  describe b:
    DefaultSpmdVector<double>{dim=500}
  
  describe x:
    DefaultSpmdVector<double>{dim=500}
  
  Iter = 0, ||b-A*x||/||b-A*x0|| = 1.000000e+00
  Iter = 31, ||b-A*x||/||b-A*x0|| = 3.172143e-01
  Iter = 62, ||b-A*x||/||b-A*x0|| = 4.876230e-02
  Iter = 93, ||b-A*x||/||b-A*x0|| = 1.090262e-02
  Iter = 124, ||b-A*x||/||b-A*x0|| = 2.477719e-03
  Iter = 155, ||b-A*x||/||b-A*x0|| = 8.999423e-04
  Iter = 186, ||b-A*x||/||b-A*x0|| = 1.621384e-04
  Iter = 202, ||b-A*x||/||b-A*x0|| = 9.565884e-05

Checking the residual ourselves ...
  
  ||b-A*x||/||b|| = 1.215596e-03/1.270762e+01 = 9.565884e-05 <= 10.0*tolerance = 1.000000e-03: passed

Total time = 7.735100e-02 sec

***
*** Running silly CG solver using scalar type = 'std::complex<float>' ...
***

Constructing tridiagonal matrix A of local dimension = 500 and diagonal multiplier = (1.001000e+00,0.000000e+00) ...

Global dimension of A = 500

Testing the constructed linear operator A ...
  
  *** Entering LinearOpTester<std::complex<float>,std::complex<float>>::check(op,...) ...
  
  describe op:
    ExampleTridiagSpmdLinearOp<std::complex<float>>{rangeDim=500,domainDim=500}
  
  Checking the domain and range spaces ... passed!
  
  this->check_linear_properties()==true: Checking the linear properties of the forward linear operator ... passed!
  
  (this->check_linear_properties()&&this->check_adjoint())==false: Skipping the check of the linear properties of the adjoint operator!
  
  this->check_adjoint()==false: Skipping check for the agreement of the adjoint and forward operators!
  
  this->check_for_symmetry()==true: Performing check of symmetry ... passed!
  
  Congratulations, this LinearOpBase object seems to check out!
  
  *** Leaving LinearOpTester<std::complex<float>,std::complex<float>>::check(...)

Solving the linear system with sillyCgSolve(...) ...
  
  Starting CG solver ...
  
  describe A:
    ExampleTridiagSpmdLinearOp<std::complex<float>>{rangeDim=500,domainDim=500}
  
  describe b:
    DefaultSpmdVector<std::complex<float>>{dim=500}
  
  describe x:
    DefaultSpmdVector<std::complex<float>>{dim=500}
  
  Iter = 0, ||b-A*x||/||b-A*x0|| = 1.000000e+00
  Iter = 31, ||b-A*x||/||b-A*x0|| = 3.578967e-01
  Iter = 62, ||b-A*x||/||b-A*x0|| = 6.437086e-02
  Iter = 93, ||b-A*x||/||b-A*x0|| = 1.508291e-02
  Iter = 124, ||b-A*x||/||b-A*x0|| = 3.664470e-03
  Iter = 155, ||b-A*x||/||b-A*x0|| = 9.545368e-04
  Iter = 186, ||b-A*x||/||b-A*x0|| = 2.053413e-04
  Iter = 202, ||b-A*x||/||b-A*x0|| = 9.788751e-05

Checking the residual ourselves ...
  
  ||b-A*x||/||b|| = 1.914953e-03/1.818107e+01 = 1.053268e-04 <= 10.0*tolerance = 9.999999e-04: passed

Total time = 1.590560e-01 sec

***
*** Running silly CG solver using scalar type = 'std::complex<double>' ...
***

Constructing tridiagonal matrix A of local dimension = 500 and diagonal multiplier = (1.001000e+00,0.000000e+00) ...

Global dimension of A = 500

Testing the constructed linear operator A ...
  
  *** Entering LinearOpTester<std::complex<double>,std::complex<double>>::check(op,...) ...
  
  describe op:
    ExampleTridiagSpmdLinearOp<std::complex<double>>{rangeDim=500,domainDim=500}
  
  Checking the domain and range spaces ... passed!
  
  this->check_linear_properties()==true: Checking the linear properties of the forward linear operator ... passed!
  
  (this->check_linear_properties()&&this->check_adjoint())==false: Skipping the check of the linear properties of the adjoint operator!
  
  this->check_adjoint()==false: Skipping check for the agreement of the adjoint and forward operators!
  
  this->check_for_symmetry()==true: Performing check of symmetry ... passed!
  
  Congratulations, this LinearOpBase object seems to check out!
  
  *** Leaving LinearOpTester<std::complex<double>,std::complex<double>>::check(...)

Solving the linear system with sillyCgSolve(...) ...
  
  Starting CG solver ...
  
  describe A:
    ExampleTridiagSpmdLinearOp<std::complex<double>>{rangeDim=500,domainDim=500}
  
  describe b:
    DefaultSpmdVector<std::complex<double>>{dim=500}
  
  describe x:
    DefaultSpmdVector<std::complex<double>>{dim=500}
  
  Iter = 0, ||b-A*x||/||b-A*x0|| = 1.000000e+00
  Iter = 31, ||b-A*x||/||b-A*x0|| = 3.579080e-01
  Iter = 62, ||b-A*x||/||b-A*x0|| = 6.437509e-02
  Iter = 93, ||b-A*x||/||b-A*x0|| = 1.508444e-02
  Iter = 124, ||b-A*x||/||b-A*x0|| = 3.664968e-03
  Iter = 155, ||b-A*x||/||b-A*x0|| = 9.546967e-04
  Iter = 186, ||b-A*x||/||b-A*x0|| = 2.053825e-04
  Iter = 202, ||b-A*x||/||b-A*x0|| = 9.790892e-05

Checking the residual ourselves ...
  
  ||b-A*x||/||b|| = 1.780089e-03/1.818107e+01 = 9.790892e-05 <= 10.0*tolerance = 1.000000e-03: passed

Total time = 1.790050e-01 sec

All of the tests seem to have run successfully!

This example program also takes a number of command-line options. To see what the command-line options are, use the --help option. The command-line options returned from ./sillyCgSolve_mpi.exe --echo-command-line --help are:

Echoing the command-line:

../example/operator_vector/sillyCgSolve_mpi.exe --echo-command-line --help 

Usage: ../example/operator_vector/sillyCgSolve_mpi.exe [options]
  options:
  --help                                Prints this help message
  --pause-for-debugging                 Pauses for user input to allow attaching a debugger
  --echo-command-line                   Echo the command-line but continue as normal
  --output-all-front-matter     bool    Set if all front matter is printed to the default FancyOStream or not
  --output-no-front-matter              (default: --output-no-front-matter)
  --output-show-line-prefix     bool    Set if the line prefix matter is printed to the default FancyOStream or not
  --output-no-show-line-prefix          (default: --output-no-show-line-prefix)
  --output-show-tab-count       bool    Set if the tab count is printed to the default FancyOStream or not
  --output-no-show-tab-count            (default: --output-no-show-tab-count)
  --output-show-proc-rank       bool    Set if the processor rank is printed to the default FancyOStream or not
  --output-no-show-proc-rank            (default: --output-no-show-proc-rank)
  --output-to-root-rank-only    int     Set which processor (the root) gets the output.  If < 0, then all processors get output.
                                        (default: --output-to-root-rank-only=0)
  --verbose                     bool    Determines if any output is printed or not.
  --quiet                               (default: --verbose)
  --local-dim                   int     Local dimension of the linear system.
                                        (default: --local-dim=500)
  --diag-scale                  double  Scaling of the diagonal to improve conditioning.
                                        (default: --diag-scale=1.001)
  --show-all-tests              bool    Show all LinearOpTester tests or not
  --show-summary-only                   (default: --show-summary-only)
  --tol                         double  Relative tolerance for linear system solve.
                                        (default: --tol=0.0001)
  --max-num-iters               int     Maximum of CG iterations.
                                        (default: --max-num-iters=300)
  --dump-all                    bool    Determines if vectors are printed or not.
  --no-dump-all                         (default: --no-dump-all)
Teuchos::GlobalMPISession::GlobalMPISession(): started serial run


Classes
class	ExampleTridiagSpmdLinearOp< Scalar >
	Simple example subclass for Spmd tridiagonal matrices. More...

Templated SPMD Implementation of the CG Method
[CG Examples]

Classes

Detailed Description

Templated SPMD Implementation of the CG Method [CG Examples]

Classes

Detailed Description

Templated SPMD Implementation of the CG Method
[CG Examples]