Main Page | Modules | Class Hierarchy | Alphabetical List | Class List | Directories | File List | Class Members | Examples | Search for 

Templated Serial Implementation of the CG Method
[CG Examples]

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

Classes

class  ExampleTridiagSerialLinearOp< Scalar >
 Simple example subclass for serial tridiagonal matrices. More...

Detailed Description

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

The class ExampleTridiagSerialLinearOp that derives from the base class Thyra::SpmdLinearOpBase is quite simple and its complete implementation looks like:

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

  Thyra::Index         dim_;
  std::vector<Scalar>  lower_;   // size = dim - 1    
  std::vector<Scalar>  diag_;    // size = dim
  std::vector<Scalar>  upper_;   // size = dim - 1

public:

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

  ExampleTridiagSerialLinearOp() : dim_(0) {}

  ExampleTridiagSerialLinearOp( const Thyra::Index dim, const Scalar lower[], const Scalar diag[], const Scalar upper[] )
    { this->initialize(dim,lower,diag,upper); }
  
  void initialize(
    const Thyra::Index              dim      // >= 2
    ,const Scalar                   lower[]  // size == dim - 1
    ,const Scalar                   diag[]   // size == dim
    ,const Scalar                   upper[]  // size == dim - 1
    )
    {
      TEST_FOR_EXCEPT( dim < 2 );
      this->setLocalDimensions(Teuchos::null,dim,dim); // Needed to setup range() and domain()
      dim_ = dim;
      lower_.resize(dim-1);  for( int k = 0; k < dim-1; ++k ) lower_[k] = lower[k];
      diag_.resize(dim);     for( int k = 0; k < dim;   ++k ) diag_[k]  = diag[k];
      upper_.resize(dim-1);  for( int k = 0; k < dim-1; ++k ) upper_[k] = upper[k];
    }

  // Overridden form Teuchos::Describable */

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

protected:

  // Overridden from SingleScalarEuclideanLinearOpBase

  bool opSupported(Thyra::ETransp M_trans) const { return true; }  // This class supports everything!

  // Overridden from SpmdLinearOpBase

  void euclideanApply(
    const Thyra::ETransp                         M_trans
    ,const RTOpPack::ConstSubVectorView<Scalar>  &x_in
    ,const RTOpPack::SubVectorView<Scalar>       *y_out
    ,const Scalar                                alpha
    ,const Scalar                                beta
    ) const
    {
      typedef Teuchos::ScalarTraits<Scalar> ST;
      // Get raw pointers to the values
      const Scalar *x = x_in.values();
      Scalar       *y = y_out->values();
      // Perform y = beta*y (being careful to set y=0 if beta=0 in case y is uninitialized on input!)
      Thyra::Index k = 0;
      if( beta == ST::zero() ) {
        for( k = 0; k < dim_; ++k ) y[k] = ST::zero();
      }
      else if( beta != ST::one() ) {
        for( k = 0; k < dim_; ++k ) y[k] *= beta;
      }
      // Perform y = alpha*op(M)*x 
      k = 0;
      if( M_trans == Thyra::NOTRANS ) {
        y[k] += alpha * ( diag_[k]*x[k] + upper_[k]*x[k+1] );                         // First row
        for( k = 1; k < dim_ - 1; ++k )
          y[k] += alpha * ( lower_[k-1]*x[k-1] + diag_[k]*x[k] + upper_[k]*x[k+1] );  // Middle rows
        y[k] += alpha * ( lower_[k-1]*x[k-1] + diag_[k]*x[k] );                       // Last row
      }
      else if( M_trans == Thyra::CONJ ) {
        y[k] += alpha * ( ST::conjugate(diag_[k])*x[k] + ST::conjugate(upper_[k])*x[k+1] );
        for( k = 1; k < dim_ - 1; ++k )
          y[k] += alpha * ( ST::conjugate(lower_[k-1])*x[k-1] + ST::conjugate(diag_[k])*x[k] + ST::conjugate(upper_[k])*x[k+1] );
        y[k] += alpha * ( ST::conjugate(lower_[k-1])*x[k-1] + ST::conjugate(diag_[k])*x[k] );
      }
      else if( M_trans == Thyra::TRANS ) {
        y[k] += alpha * ( diag_[k]*x[k] + lower_[k]*x[k+1] );
        for( k = 1; k < dim_ - 1; ++k )
          y[k] += alpha * ( upper_[k-1]*x[k-1] + diag_[k]*x[k] + lower_[k]*x[k+1] );
        y[k] += alpha * ( upper_[k-1]*x[k-1] + diag_[k]*x[k] );
      }
      else if( M_trans == Thyra::CONJTRANS ) {
        y[k] += alpha * ( ST::conjugate(diag_[k])*x[k] + ST::conjugate(lower_[k])*x[k+1] );
        for( k = 1; k < dim_ - 1; ++k )
          y[k] += alpha * ( ST::conjugate(upper_[k-1])*x[k-1] + ST::conjugate(diag_[k])*x[k] + ST::conjugate(lower_[k])*x[k+1] );
        y[k] += alpha * ( ST::conjugate(upper_[k-1])*x[k-1] + ST::conjugate(diag_[k])*x[k] );
      }
      else {
        TEST_FOR_EXCEPT(true); // Throw exception if we get here!
      }
    }

};  // end class ExampleTridiagSerialLinearOp

The above serial matrix class is used in an example program (see runCgSolveExample() below) that calls sillyCgSolve() or sillierCgSolve(). In this example program, the matrix constructed and used is the following tridiagonal matrix

\[ A= \left[\begin{array}{rrrrrrrrrr} 2 a & -1 \\ -r(1) & 2 a & -1 \\ & \ddots & \ddots & \ddots \\ & & -r(n-2) & 2 a & -1 \\ & & & -r(n-1) & 2 a \end{array}\right] \]

where $a$ is an adjustable diagonal scale factories that makes the matrix $A$ more or less well conditioned and $r(i)$ is either $-1$ for a symmetric operator or $rand()$ for an unsymmetric operator.

If a symmetric operator is used, then CG is run using $A$ directly. If $A$ is unsymmetric, then the normal equations

\[ A^H A x = A^H b \]

are solved and the operator used is

\[ A \Rightarrow A^H A \]

The CG method is then run on the matrix $A$ or $A^H 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 int                                                      dim
  ,const Scalar                                                  diagScale
  ,const bool                                                    symOp
  ,const bool                                                    showAllTests
  ,const bool                                                    verbose
  ,const typename Teuchos::ScalarTraits<Scalar>::magnitudeType   tolerance
  ,const int                                                     maxNumIters
  ,const bool                                                    useSillierCg
  )
{
  using Teuchos::RefCountPtr; using Teuchos::rcp; using Teuchos::null;
  using Teuchos::OSTab;
  typedef Teuchos::ScalarTraits<Scalar> ST;
  using Thyra::multiply; using Thyra::scale;
  typedef typename ST::magnitudeType    ScalarMag;
  bool success = true;
  bool result;
  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                         ]
  //       [ -r(1)  a*2       -1                ]
  //  A =  [          .        .        .       ]
  //       [             -r(n-2)      a*2    -1 ]
  //       [                      -r(n-1)   a*2 ]
  //
  // (A.1) Create the tridiagonal matrix operator
  if(verbose) *out << "\nConstructing tridiagonal matrix A of dimension = " << dim << " and diagonal multiplier = " << diagScale << " ...\n";
  std::vector<Scalar> lower(dim-1), diag(dim), upper(dim-1);
  const Scalar up = -ST::one(), diagTerm = Scalar(2)*diagScale*ST::one(), low = -(symOp?ST::one():ST::random());
  int k = 0;
  diag[k] = diagTerm; upper[k] = up;                          // First row
  for( k = 1; k < dim - 1; ++k ) {
    lower[k-1] = low; diag[k] = diagTerm; upper[k] = up;      // Middle rows
  }
  lower[k-1] = low; diag[k] = diagTerm;                       // Last row
  RefCountPtr<const Thyra::LinearOpBase<Scalar> >
    A = rcp(new ExampleTridiagSerialLinearOp<Scalar>(dim,&lower[0],&diag[0],&upper[0]));
  // (A.2) Testing the linear operator constructed linear operator
  if(verbose) *out << "\nTesting the constructed linear operator A ...\n";
  Thyra::LinearOpTester<Scalar> linearOpTester;
  linearOpTester.enable_all_tests(false);            // DEBUG ONLY
  linearOpTester.check_linear_properties(true);      // DEBUG ONLY
  linearOpTester.set_all_error_tol(tolerance);
  linearOpTester.set_all_warning_tol(ScalarMag(ScalarMag(1e-2)*tolerance));
  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() );
  // (A.5) Create the final linear system
  if(!symOp) {
    if(verbose) *out << "\nSetting up normal equations for unsymmetric system A^H*(A*x-b) => new A*x = b ...\n";
    RefCountPtr<const Thyra::LinearOpBase<Scalar> > AtA = multiply(adjoint(A),A);      // A^H*A
    RefCountPtr<Thyra::VectorBase<Scalar> >         nb = createMember(AtA->range());   // A^H*b
    Thyra::apply(*A,Thyra::CONJTRANS,*b,&*nb);
    A = AtA;
    b = nb;
  }
  // (A.6) Testing the linear operator used with the solve
  if(verbose) *out << "\nTesting the linear operator used with the solve ...\n";
  linearOpTester.check_for_symmetry(true);
  result = linearOpTester.check(*A,out.get());
  if(!result) success = false;
  //
  // (B) Solve the linear system with the silly CG solver
  //
  if(verbose) *out << "\nSolving the linear system with sillyCgSolve(...) ...\n";
  if(useSillierCg)
    result = sillierCgSolve(*A,*b,maxNumIters,tolerance,&*x,OSTab(out).getOStream().get());
  else
    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::RefCountPtr<Teuchos::FancyOStream>
    out = Teuchos::VerboseObjectBase::getDefaultOStream();

  try {

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

    int    dim          = 500;
    double diagScale    = 1.001;
    bool   symOp        = true;
    bool   showAllTests = false;
    double tolerance    = 1e-4;
    int    maxNumIters  = 300;
    bool   useSillierCg = false;

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

    clp.setOption( "verbose", "quiet", &verbose, "Determines if any output is printed or not." );
    clp.setOption( "dim", &dim, "Dimension of the linear system." );
    clp.setOption( "diag-scale", &diagScale, "Scaling of the diagonal to improve conditioning." );
    clp.setOption( "sym-op", "unsym-op", &symOp, "Determines if the operator is symmetric or not." );
    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( "use-sillier-cg", "use-silly-cg", &useSillierCg,
                   "Use the handle-based sillerCgSolve() function or the nonhandle-based sillyCgSolve() function");

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

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

    // Run using float
    result = runCgSolveExample<float>(dim,diagScale,symOp,showAllTests,verbose,tolerance,maxNumIters,useSillierCg);
    if(!result) success = false;

    // Run using double
    result = runCgSolveExample<double>(dim,diagScale,symOp,showAllTests,verbose,tolerance,maxNumIters,useSillierCg);
    if(!result) success = false;

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

    // Run using std::complex<float>
    result = runCgSolveExample<std::complex<float> >(dim,diagScale,symOp,showAllTests,verbose,tolerance,maxNumIters,useSillierCg);
    if(!result) success = false;

    // Run using std::complex<double>
    result = runCgSolveExample<std::complex<double> >(dim,diagScale,symOp,showAllTests,verbose,tolerance,maxNumIters,useSillierCg);
    if(!result) success = false;

#endif

#ifdef HAVE_TEUCHOS_GNU_MP

    // Run using mpf_class
    result = runCgSolveExample<mpf_class>(dim,diagScale,symOp,showAllTests,verbose,tolerance,maxNumIters,useSillierCg);
    if(!result) success = false;

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

    // Run using std::complex<mpf_class>
    //result = runCgSolveExample<std::complex<mpf_class> >(dim,mpf_class(diagScale),symOp,showAllTests,verbose,mpf_class(tolerance),maxNumIters);
    //if(!result) success = false;
    //The above commented-out code throws a floating-point exception?

#endif

#endif

  }
  TEUCHOS_STANDARD_CATCH_STATEMENTS(true,*out,success)

  if (verbose) {
    if(success)   *out << "\nCongratulations! All of the tests checked out!\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_serial.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 and produces the following output:

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

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

Testing the constructed linear operator A ...
  
  *** Entering LinearOpTester<float,float>::check(op,...) ...
  
  describe op:
    ExampleTridiagSerialLinearOp<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()==false: Skipping check of symmetry ...
  
  Congratulations, this LinearOpBase object seems to check out!
  
  *** Leaving LinearOpTester<float,float>::check(...)

Testing the linear operator used with the solve ...
  
  *** Entering LinearOpTester<float,float>::check(op,...) ...
  
  describe op:
    ExampleTridiagSerialLinearOp<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:
    ExampleTridiagSerialLinearOp<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 = 8.771700e-02 sec

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

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

Testing the constructed linear operator A ...
  
  *** Entering LinearOpTester<double,double>::check(op,...) ...
  
  describe op:
    ExampleTridiagSerialLinearOp<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()==false: Skipping check of symmetry ...
  
  Congratulations, this LinearOpBase object seems to check out!
  
  *** Leaving LinearOpTester<double,double>::check(...)

Testing the linear operator used with the solve ...
  
  *** Entering LinearOpTester<double,double>::check(op,...) ...
  
  describe op:
    ExampleTridiagSerialLinearOp<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:
    ExampleTridiagSerialLinearOp<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 = 8.638900e-02 sec

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

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

Testing the constructed linear operator A ...
  
  *** Entering LinearOpTester<std::complex<float>,std::complex<float>>::check(op,...) ...
  
  describe op:
    ExampleTridiagSerialLinearOp<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()==false: Skipping check of symmetry ...
  
  Congratulations, this LinearOpBase object seems to check out!
  
  *** Leaving LinearOpTester<std::complex<float>,std::complex<float>>::check(...)

Testing the linear operator used with the solve ...
  
  *** Entering LinearOpTester<std::complex<float>,std::complex<float>>::check(op,...) ...
  
  describe op:
    ExampleTridiagSerialLinearOp<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:
    ExampleTridiagSerialLinearOp<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.703310e-01 sec

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

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

Testing the constructed linear operator A ...
  
  *** Entering LinearOpTester<std::complex<double>,std::complex<double>>::check(op,...) ...
  
  describe op:
    ExampleTridiagSerialLinearOp<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()==false: Skipping check of symmetry ...
  
  Congratulations, this LinearOpBase object seems to check out!
  
  *** Leaving LinearOpTester<std::complex<double>,std::complex<double>>::check(...)

Testing the linear operator used with the solve ...
  
  *** Entering LinearOpTester<std::complex<double>,std::complex<double>>::check(op,...) ...
  
  describe op:
    ExampleTridiagSerialLinearOp<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:
    ExampleTridiagSerialLinearOp<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 = 2.138060e-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. The command-line options for this program returned from ./sillyCgSolve_serial.exe --echo-command-line --help are:

Echoing the command-line:

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

Usage: ../example/operator_vector/sillyCgSolve_serial.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)
  --dim                         int     Dimension of the linear system.
                                        (default: --dim=500)
  --diag-scale                  double  Scaling of the diagonal to improve conditioning.
                                        (default: --diag-scale=1.001)
  --sym-op                      bool    Determines if the operator is symmetric or not.
  --unsym-op                            (default: --sym-op)
  --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)
  --use-sillier-cg              bool    Use the handle-based sillerCgSolve() function or the nonhandle-based sillyCgSolve() function
  --use-silly-cg                        (default: --use-silly-cg)

When the option --unsym-op is selected, the normal equations are solved which is shown in the following example:

Echoing the command-line:

../example/operator_vector/sillyCgSolve_serial.exe --echo-command-line --unsym-op 


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

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

Testing the constructed linear operator A ...
  
  *** Entering LinearOpTester<float,float>::check(op,...) ...
  
  describe op:
    ExampleTridiagSerialLinearOp<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()==false: Skipping check of symmetry ...
  
  Congratulations, this LinearOpBase object seems to check out!
  
  *** Leaving LinearOpTester<float,float>::check(...)

Setting up normal equations for unsymmetric system A^H*(A*x-b) => new A*x = b ...

Testing the linear operator used with the solve ...
  
  *** Entering LinearOpTester<float,float>::check(op,...) ...
  
  describe op:
    Thyra::DefaultMultipliedLinearOp<float>{numOps = 2}
  
  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:
    Thyra::DefaultMultipliedLinearOp<float>{rangeDim=500,domainDim=500}
      numOps=2
      Constituent LinearOpBase objects for M = Op[0]*...*Op[numOps-1]:
        Op[0] =
          DefaultScaledAdjointLinearOp<float>{rangeDim=500,domainDim=500}
            overallScalar=1.000000e+00
            overallTransp=CONJTRANS
            Constituent transformations:
              transp=CONJTRANS
                origOp =
                  ExampleTridiagSerialLinearOp<float>{rangeDim=500,domainDim=500}
        Op[1] =
          ExampleTridiagSerialLinearOp<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|| = 1.969999e-03
  Iter = 48, ||b-A*x||/||b-A*x0|| = 9.033530e-05

Checking the residual ourselves ...
  
  ||b-A*x||/||b|| = 2.678777e-03/2.966171e+01 = 9.031095e-05 <= 10.0*tolerance = 9.999999e-04: passed

Total time = 5.998500e-02 sec

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

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

Testing the constructed linear operator A ...
  
  *** Entering LinearOpTester<double,double>::check(op,...) ...
  
  describe op:
    ExampleTridiagSerialLinearOp<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()==false: Skipping check of symmetry ...
  
  Congratulations, this LinearOpBase object seems to check out!
  
  *** Leaving LinearOpTester<double,double>::check(...)

Setting up normal equations for unsymmetric system A^H*(A*x-b) => new A*x = b ...

Testing the linear operator used with the solve ...
  
  *** Entering LinearOpTester<double,double>::check(op,...) ...
  
  describe op:
    Thyra::DefaultMultipliedLinearOp<double>{numOps = 2}
  
  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:
    Thyra::DefaultMultipliedLinearOp<double>{rangeDim=500,domainDim=500}
      numOps=2
      Constituent LinearOpBase objects for M = Op[0]*...*Op[numOps-1]:
        Op[0] =
          DefaultScaledAdjointLinearOp<double>{rangeDim=500,domainDim=500}
            overallScalar=1.000000e+00
            overallTransp=CONJTRANS
            Constituent transformations:
              transp=CONJTRANS
                origOp =
                  ExampleTridiagSerialLinearOp<double>{rangeDim=500,domainDim=500}
        Op[1] =
          ExampleTridiagSerialLinearOp<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|| = 9.488548e-05

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

Total time = 5.385100e-02 sec

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

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

Testing the constructed linear operator A ...
  
  *** Entering LinearOpTester<std::complex<float>,std::complex<float>>::check(op,...) ...
  
  describe op:
    ExampleTridiagSerialLinearOp<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()==false: Skipping check of symmetry ...
  
  Congratulations, this LinearOpBase object seems to check out!
  
  *** Leaving LinearOpTester<std::complex<float>,std::complex<float>>::check(...)

Setting up normal equations for unsymmetric system A^H*(A*x-b) => new A*x = b ...

Testing the linear operator used with the solve ...
  
  *** Entering LinearOpTester<std::complex<float>,std::complex<float>>::check(op,...) ...
  
  describe op:
    Thyra::DefaultMultipliedLinearOp<std::complex<float>>{numOps = 2}
  
  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:
    Thyra::DefaultMultipliedLinearOp<std::complex<float>>{rangeDim=500,domainDim=500}
      numOps=2
      Constituent LinearOpBase objects for M = Op[0]*...*Op[numOps-1]:
        Op[0] =
          DefaultScaledAdjointLinearOp<std::complex<float>>{rangeDim=500,domainDim=500}
            overallScalar=(1.000000e+00,0.000000e+00)
            overallTransp=CONJTRANS
            Constituent transformations:
              transp=CONJTRANS
                origOp =
                  ExampleTridiagSerialLinearOp<std::complex<float>>{rangeDim=500,domainDim=500}
        Op[1] =
          ExampleTridiagSerialLinearOp<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 = 22, ||b-A*x||/||b-A*x0|| = 8.391360e-05

Checking the residual ourselves ...
  
  ||b-A*x||/||b|| = 3.515154e-03/4.188961e+01 = 8.391472e-05 <= 10.0*tolerance = 9.999999e-04: passed

Total time = 6.419600e-02 sec

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

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

Testing the constructed linear operator A ...
  
  *** Entering LinearOpTester<std::complex<double>,std::complex<double>>::check(op,...) ...
  
  describe op:
    ExampleTridiagSerialLinearOp<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()==false: Skipping check of symmetry ...
  
  Congratulations, this LinearOpBase object seems to check out!
  
  *** Leaving LinearOpTester<std::complex<double>,std::complex<double>>::check(...)

Setting up normal equations for unsymmetric system A^H*(A*x-b) => new A*x = b ...

Testing the linear operator used with the solve ...
  
  *** Entering LinearOpTester<std::complex<double>,std::complex<double>>::check(op,...) ...
  
  describe op:
    Thyra::DefaultMultipliedLinearOp<std::complex<double>>{numOps = 2}
  
  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:
    Thyra::DefaultMultipliedLinearOp<std::complex<double>>{rangeDim=500,domainDim=500}
      numOps=2
      Constituent LinearOpBase objects for M = Op[0]*...*Op[numOps-1]:
        Op[0] =
          DefaultScaledAdjointLinearOp<std::complex<double>>{rangeDim=500,domainDim=500}
            overallScalar=(1.000000e+00,0.000000e+00)
            overallTransp=CONJTRANS
            Constituent transformations:
              transp=CONJTRANS
                origOp =
                  ExampleTridiagSerialLinearOp<std::complex<double>>{rangeDim=500,domainDim=500}
        Op[1] =
          ExampleTridiagSerialLinearOp<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 = 10, ||b-A*x||/||b-A*x0|| = 7.136897e-05

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

Total time = 5.536900e-02 sec

Congratulations! All of the tests checked out!

Note in the above example how the normal operator $A^H A$ is described. This aggregate operator is created by the function calls Thyra::adjoint() and Thyra::multiply() which create implicit Thyra::ScaledAdjointedLinearOp and Thyra::DefaultMultipliedLinearOp objects.

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

Generated on Thu Sep 18 12:32:32 2008 for Thyra Operator/Vector Support by doxygen 1.3.9.1