Templated MPI Implementation of the CG Method
[CG Examples]

Here is an example program that shows the use of the example MPI templated matrix class MPITridiagLinearOp with the example linear ANA implementation sillyCgSolve(). More...

Classes
class	MPITridiagLinearOp< Scalar >
	Simple example subclass for MPI tridiagonal matrices. More...

---

Detailed Description

Here is an example program that shows the use of the example MPI templated matrix class MPITridiagLinearOp with the example linear ANA implementation sillyCgSolve().

The class MPITridiagLinearOp that derives from the base class Thyra::MPILinearOpBase is quite simple and its implementation (minus the communication() function) looks like:

template<class Scalar>
class MPITridiagLinearOp : public Thyra::MPILinearOpBase<Scalar> {
private:

  MPI_Comm             mpiComm_;
  int                  procRank_;
  int                  numProc_;
  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::MPILinearOpBase<Scalar>::euclideanApply;

  MPITridiagLinearOp() : mpiComm_(MPI_COMM_NULL), procRank_(0), numProc_(0) {}

  MPITridiagLinearOp( MPI_Comm mpiComm, const Thyra::Index localDim, const Scalar lower[], const Scalar diag[], const Scalar upper[] )
    { this->initialize(mpiComm,localDim,lower,diag,upper);  }
  
  void initialize(
    MPI_Comm                        mpiComm
    ,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(mpiComm,localDim,localDim); // We must tell the base class our local dimensions to setup range() and domain()
      mpiComm_  = mpiComm;
      localDim_ = localDim;
      MPI_Comm_size( mpiComm, &numProc_ );
      MPI_Comm_rank( mpiComm, &procRank_ );
      const Thyra::Index
        lowerDim = ( procRank_ == 0          ? localDim - 1 : localDim ),
        upperDim = ( procRank_ == numProc_-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("MPITridiagLinearOp<") + 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::SubVectorT<Scalar>          &local_x_in
    ,const RTOpPack::MutableSubVectorT<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 processor we are
      const bool first = ( procRank_ == 0 ), last = ( procRank_ == numProc_-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 MPITridiagLinearOp

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

where is an adjustable diagonal scale factor that makes the matrix more or less well conditioned.

The CG method is then run on the matrix 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(
  MPI_Comm                                                       mpiComm
  ,const int                                                     procRank
  ,const int                                                     numProc
  ,const int                                                     localDim
  ,const Scalar                                                  diagScale
  ,const bool                                                    verbose
  ,const bool                                                    dumpAll
  ,const typename Teuchos::ScalarTraits<Scalar>::magnitudeType   tolerance
  ,const int                                                     maxNumIters
  )
{
  using Teuchos::RefCountPtr;
  using Teuchos::rcp;
  typedef Teuchos::ScalarTraits<Scalar> ST;
  typedef typename ST::magnitudeType    ScalarMag;
  bool success = true;
  bool result;
  // Setup the output stream (do output only on root process!)
  Teuchos::oblackholestream black_hole_out;
  std::ostream &out = ( procRank == 0 ? std::cout : black_hole_out );
  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 MPITridiagLinearOp<Scalar>(mpiComm,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);
  result = linearOpTester.check(*A,verbose?&out:0);
  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
  //
  result = sillyCgSolve(*A,*b,maxNumIters,tolerance,&*x,verbose?&out:0);
  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
      << "\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:

• Real types float and double
• Complex real types std::complex<float> and std::complex<double> (if --enable-teuchos-complex was used at configuration time)

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;

  MPI_Init(&argc,&argv);

  MPI_Comm mpiComm = MPI_COMM_WORLD;
  int procRank, numProc;
  MPI_Comm_size( mpiComm, &numProc );
  MPI_Comm_rank( mpiComm, &procRank );

  try {

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

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

    CommandLineProcessor  clp(false); // Don't throw exceptions

    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( "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>(mpiComm,procRank,numProc,localDim,diagScale,verbose,dumpAll,tolerance,maxNumIters);
    if(!result) success = false;

    // Run using double
    result = runCgSolveExample<double>(mpiComm,procRank,numProc,localDim,diagScale,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> >(mpiComm,procRank,numProc,localDim,diagScale,verbose,dumpAll,tolerance,maxNumIters);
    if(!result) success = false;

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

#endif    

  }
  catch( const std::exception &excpt ) {
    std::cerr << "*** p="<<procRank<<": Caught standard exception : " << excpt.what() << std::endl;
    success = false;
  }
  catch( ... ) {
    std::cerr << "*** p="<<procRank<<":Caught an unknown exception\n";
    success = false;
  }

  if( verbose && procRank==0 ) {
    if(success) std::cout << "\nAll of the tests seem to have run successfully!\n";
    else        std::cout << "\nOh no! at least one of the tests failed!\n";  
  }
  
  MPI_Finalize();
  
  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/Core/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, at the time of this writing, for one processor produces the following output:

$ ./sillyCgSolve_mpi.exe

***
*** 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<Scalar>::check(op,...) ...

op = 18MPITridiagLinearOpIfE

Checking the domain and range spaces ...
op.domain().get() != NULL ? passed
op.range().get() != NULL ? passed

Checking that the operator truly is linear:

  0.5*op*(v1 + v2) == 0.5*op*v1 + 0.5*op*v2
          \_____/         \___/
             v3            v5
  \_____________/     \___________________/
         v4                    v5

           sum(v4) == sum(v5)

Random vector tests = 1

v1 = randomize(-1,+1); ...

v2 = randomize(-1,+1); ...

v3 = v1 + v2 ...

v4 = 0.5*op*v3 ...

v5 = op*v1 ...

v5 = 0.5*op*v2 + 0.5*v5 ...

Check: rel_err(sum(v4),sum(v5))
       = rel_err(-0.0412178,-0.0412174) = 1.0936e-05
         <= linear_properties_error_tol() = 0.0001 : passed
Warning! rel_err(sum(v4),sum(v5))
       = rel_err(-0.0412178,-0.0412174) = 1.0936e-05
         >= linear_properties_warning_tol() = 1e-06!

Skipping adjoint check since op.opSupported(CONJTRANS)==false

Performing check of symmetry since check_for_symmetry()==true ...

op.domain()->isCompatible(*op.range()) == true : passed

Checking that the operator is symmetric as:

  <0.5*op*v2,v1> == <v2,0.5*op*v1>
   \_______/            \_______/
      v4                    v3

         <v4,v1> == <v2,v3>

Random vector tests = 1

v1 = randomize(-1,+1); ...

v2 = randomize(-1,+1); ...

v3 = 0.5*op*v1 ...

v4 = 0.5*op*v2 ...

Check: rel_err(<v4,v1>,<v2,v3>)
       = rel_err(-15.3963,-15.3963) = 3.09709e-07
         <= symmetry_error_tol() = 0.0001 : passed

*** Leaving LinearOpTester<Scalar>::check(...)

Starting CG solver ...

describe A:
  type = '18MPITridiagLinearOpIfE', rangeDim = 500, domainDim = 500

describe b:
  type = 'N5Thyra12MPIVectorStdIfEE', size = 500

describe x:
  type = 'N5Thyra12MPIVectorStdIfEE', size = 500

Iter = 0, ||b-A*x||/||b-A*x0|| = 1.000000e+00
Iter = 31, ||b-A*x||/||b-A*x0|| = 2.812529e-01
Iter = 62, ||b-A*x||/||b-A*x0|| = 6.589733e-02
Iter = 93, ||b-A*x||/||b-A*x0|| = 1.381910e-02
Iter = 124, ||b-A*x||/||b-A*x0|| = 3.546838e-03
Iter = 155, ||b-A*x||/||b-A*x0|| = 9.840684e-04
Iter = 186, ||b-A*x||/||b-A*x0|| = 1.891940e-04
Iter = 204, ||b-A*x||/||b-A*x0|| = 9.613070e-05

||b-A*x||/||b|| = 1.344796e-03/1.299397e+01 = 1.034939e-04 <= 10.0*tolerance = 9.999999e-04: passed

Total time = 5.490000e-01 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<Scalar>::check(op,...) ...

op = 18MPITridiagLinearOpIdE

Checking the domain and range spaces ...
op.domain().get() != NULL ? passed
op.range().get() != NULL ? passed

Checking that the operator truly is linear:

  0.5*op*(v1 + v2) == 0.5*op*v1 + 0.5*op*v2
          \_____/         \___/
             v3            v5
  \_____________/     \___________________/
         v4                    v5

           sum(v4) == sum(v5)

Random vector tests = 1

v1 = randomize(-1,+1); ...

v2 = randomize(-1,+1); ...

v3 = v1 + v2 ...

v4 = 0.5*op*v3 ...

v5 = op*v1 ...

v5 = 0.5*op*v2 + 0.5*v5 ...

Check: rel_err(sum(v4),sum(v5))
       = rel_err(2.165119e-01,2.165119e-01) = 4.871381e-15
         <= linear_properties_error_tol() = 1.000000e-04 : passed

Skipping adjoint check since op.opSupported(CONJTRANS)==false

Performing check of symmetry since check_for_symmetry()==true ...

op.domain()->isCompatible(*op.range()) == true : passed

Checking that the operator is symmetric as:

  <0.5*op*v2,v1> == <v2,0.5*op*v1>
   \_______/            \_______/
      v4                    v3

         <v4,v1> == <v2,v3>

Random vector tests = 1

v1 = randomize(-1,+1); ...

v2 = randomize(-1,+1); ...

v3 = 0.5*op*v1 ...

v4 = 0.5*op*v2 ...

Check: rel_err(<v4,v1>,<v2,v3>)
       = rel_err(-1.054952e+01,-1.054952e+01) = 5.051482e-16
         <= symmetry_error_tol() = 1.000000e-04 : passed

*** Leaving LinearOpTester<Scalar>::check(...)

Starting CG solver ...

describe A:
  type = '18MPITridiagLinearOpIdE', rangeDim = 500, domainDim = 500

describe b:
  type = 'N5Thyra12MPIVectorStdIdEE', size = 500

describe x:
  type = 'N5Thyra12MPIVectorStdIdEE', size = 500

Iter = 0, ||b-A*x||/||b-A*x0|| = 1.000000e+00
Iter = 31, ||b-A*x||/||b-A*x0|| = 2.812612e-01
Iter = 62, ||b-A*x||/||b-A*x0|| = 6.590140e-02
Iter = 93, ||b-A*x||/||b-A*x0|| = 1.382036e-02
Iter = 124, ||b-A*x||/||b-A*x0|| = 3.547269e-03
Iter = 155, ||b-A*x||/||b-A*x0|| = 9.842233e-04
Iter = 186, ||b-A*x||/||b-A*x0|| = 1.892296e-04
Iter = 204, ||b-A*x||/||b-A*x0|| = 9.615059e-05

||b-A*x||/||b|| = 1.249378e-03/1.299397e+01 = 9.615059e-05 <= 10.0*tolerance = 1.000000e-03: passed

Total time = 2.120000e-01 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<Scalar>::check(op,...) ...

op = 18MPITridiagLinearOpISt7complexIfEE

Checking the domain and range spaces ...
op.domain().get() != NULL ? passed
op.range().get() != NULL ? passed

Checking that the operator truly is linear:

  0.5*op*(v1 + v2) == 0.5*op*v1 + 0.5*op*v2
          \_____/         \___/
             v3            v5
  \_____________/     \___________________/
         v4                    v5

           sum(v4) == sum(v5)

Random vector tests = 1

v1 = randomize(-1,+1); ...

v2 = randomize(-1,+1); ...

v3 = v1 + v2 ...

v4 = 0.5*op*v3 ...

v5 = op*v1 ...

v5 = 0.5*op*v2 + 0.5*v5 ...

Check: rel_err(sum(v4),sum(v5))
       = rel_err((5.142615e-01,1.741409e-01),(5.142629e-01,1.741405e-01)) = 2.639293e-06
         <= linear_properties_error_tol() = 1.000000e-04 : passed
Warning! rel_err(sum(v4),sum(v5))
       = rel_err((5.142615e-01,1.741409e-01),(5.142629e-01,1.741405e-01)) = 2.639293e-06
         >= linear_properties_warning_tol() = 1.000000e-06!

Skipping adjoint check since op.opSupported(CONJTRANS)==false

Performing check of symmetry since check_for_symmetry()==true ...

op.domain()->isCompatible(*op.range()) == true : passed

Checking that the operator is symmetric as:

  <0.5*op*v2,v1> == <v2,0.5*op*v1>
   \_______/            \_______/
      v4                    v3

         <v4,v1> == <v2,v3>

Random vector tests = 1

v1 = randomize(-1,+1); ...

v2 = randomize(-1,+1); ...

v3 = 0.5*op*v1 ...

v4 = 0.5*op*v2 ...

Check: rel_err(<v4,v1>,<v2,v3>)
       = rel_err((1.242970e+01,1.590466e+01),(1.242970e+01,1.590467e+01)) = 3.809048e-07
         <= symmetry_error_tol() = 1.000000e-04 : passed

*** Leaving LinearOpTester<Scalar>::check(...)

Starting CG solver ...

describe A:
  type = '18MPITridiagLinearOpISt7complexIfEE', rangeDim = 500, domainDim = 500

describe b:
  type = 'N5Thyra12MPIVectorStdISt7complexIfEEE', size = 500

describe x:
  type = 'N5Thyra12MPIVectorStdISt7complexIfEEE', size = 500

Iter = 0, ||b-A*x||/||b-A*x0|| = 1.000000e+00
Iter = 31, ||b-A*x||/||b-A*x0|| = 2.721108e-01
Iter = 62, ||b-A*x||/||b-A*x0|| = 7.355135e-02
Iter = 93, ||b-A*x||/||b-A*x0|| = 1.733771e-02
Iter = 124, ||b-A*x||/||b-A*x0|| = 4.754171e-03
Iter = 155, ||b-A*x||/||b-A*x0|| = 9.684857e-04
Iter = 186, ||b-A*x||/||b-A*x0|| = 2.239717e-04
Iter = 208, ||b-A*x||/||b-A*x0|| = 9.721012e-05

||b-A*x||/||b|| = 1.887406e-03/1.800947e+01 = 1.048008e-04 <= 10.0*tolerance = 9.999999e-04: passed

Total time = 3.020000e-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<Scalar>::check(op,...) ...

op = 18MPITridiagLinearOpISt7complexIdEE

Checking the domain and range spaces ...
op.domain().get() != NULL ? passed
op.range().get() != NULL ? passed

Checking that the operator truly is linear:

  0.5*op*(v1 + v2) == 0.5*op*v1 + 0.5*op*v2
          \_____/         \___/
             v3            v5
  \_____________/     \___________________/
         v4                    v5

           sum(v4) == sum(v5)

Random vector tests = 1

v1 = randomize(-1,+1); ...

v2 = randomize(-1,+1); ...

v3 = v1 + v2 ...

v4 = 0.5*op*v3 ...

v5 = op*v1 ...

v5 = 0.5*op*v2 + 0.5*v5 ...

Check: rel_err(sum(v4),sum(v5))
       = rel_err((-9.677554e-01,6.222520e-01),(-9.677554e-01,6.222520e-01)) = 2.869021e-15
         <= linear_properties_error_tol() = 1.000000e-04 : passed

Skipping adjoint check since op.opSupported(CONJTRANS)==false

Performing check of symmetry since check_for_symmetry()==true ...

op.domain()->isCompatible(*op.range()) == true : passed

Checking that the operator is symmetric as:

  <0.5*op*v2,v1> == <v2,0.5*op*v1>
   \_______/            \_______/
      v4                    v3

         <v4,v1> == <v2,v3>

Random vector tests = 1

v1 = randomize(-1,+1); ...

v2 = randomize(-1,+1); ...

v3 = 0.5*op*v1 ...

v4 = 0.5*op*v2 ...

Check: rel_err(<v4,v1>,<v2,v3>)
       = rel_err((1.402400e+01,4.335187e+00),(1.402400e+01,4.335187e+00)) = 3.375439e-15
         <= symmetry_error_tol() = 1.000000e-04 : passed

*** Leaving LinearOpTester<Scalar>::check(...)

Starting CG solver ...

describe A:
  type = '18MPITridiagLinearOpISt7complexIdEE', rangeDim = 500, domainDim = 500

describe b:
  type = 'N5Thyra12MPIVectorStdISt7complexIdEEE', size = 500

describe x:
  type = 'N5Thyra12MPIVectorStdISt7complexIdEEE', size = 500

Iter = 0, ||b-A*x||/||b-A*x0|| = 1.000000e+00
Iter = 31, ||b-A*x||/||b-A*x0|| = 2.721187e-01
Iter = 62, ||b-A*x||/||b-A*x0|| = 7.355587e-02
Iter = 93, ||b-A*x||/||b-A*x0|| = 1.733925e-02
Iter = 124, ||b-A*x||/||b-A*x0|| = 4.754771e-03
Iter = 155, ||b-A*x||/||b-A*x0|| = 9.686419e-04
Iter = 186, ||b-A*x||/||b-A*x0|| = 2.240143e-04
Iter = 208, ||b-A*x||/||b-A*x0|| = 9.723106e-05

||b-A*x||/||b|| = 1.751081e-03/1.800948e+01 = 9.723106e-05 <= 10.0*tolerance = 1.000000e-03: passed

Total time = 3.110000e-01 sec

Congratulations! All of the tests checked out!

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

$ ./sillyCgSolve_mpi.exe --help
Usage: ./sillyCgSolve_mpi [options]
  options:
  --help                         Prints this help message
  --pause-for-debugging          Pauses for user input to allow attaching a debugger
  --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)
  --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)

To see the full listing of this example program click: sillyCgSolve_mpi.cpp

---