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 SerialTridiagLinearOp with the example linear ANA implementation sillyCgSolve(). More...

Classes

class  SerialTridiagLinearOp< 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 SerialTridiagLinearOp with the example linear ANA implementation sillyCgSolve().

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

template<class Scalar>
class SerialTridiagLinearOp : public Thyra::SerialLinearOpBase<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::SerialLinearOpBase<Scalar>::euclideanApply;

  SerialTridiagLinearOp() : dim_(0) {}

  SerialTridiagLinearOp( 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->setDimensions(dim,dim); // We must tell the base class our dimension 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("SerialTridiagLinearOp<") + 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 SerialLinearOpBase

  void euclideanApply(
    const Thyra::ETransp                         M_trans
    ,const RTOpPack::SubVectorT<Scalar>          &x_in
    ,const RTOpPack::MutableSubVectorT<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 SerialTridiagLinearOp

The above serial 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 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                                                    verbose
  ,const typename Teuchos::ScalarTraits<Scalar>::magnitudeType   tolerance
  ,const int                                                     maxNumIters
  )
{
  using Teuchos::RefCountPtr; using Teuchos::rcp;
  typedef Teuchos::ScalarTraits<Scalar> ST;
  using Thyra::multiply; using Thyra::scale;
  typedef typename ST::magnitudeType    ScalarMag;
  bool success = true;
  bool result;
  if(verbose)
    std::cout << "\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) std::cout << "\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 SerialTridiagLinearOp<Scalar>(dim,&lower[0],&diag[0],&upper[0]));
  // (A.2) Testing the linear operator constructed linear operator
  if(verbose) std::cout << "\nTesting the constructed linear operator A ...\n";
  Thyra::LinearOpTester<Scalar> linearOpTester;
  linearOpTester.set_all_error_tol(tolerance);
  linearOpTester.set_all_warning_tol(ScalarMag(ScalarMag(1e-2)*tolerance));
  linearOpTester.show_all_tests(true);
  result = linearOpTester.check(*A,verbose?&std::cout: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() );
  // (A.5) Create the final linear system
  if(!symOp) {
    if(verbose) std::cout << "\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) std::cout << "\nTesting the linear operator used with the solve ...\n";
  linearOpTester.check_for_symmetry(true);
  result = linearOpTester.check(*A,verbose?&std::cout:0);
  if(!result) success = false;
  //
  // (B) Solve the linear system with the silly CG solver
  //
  result = sillyCgSolve(*A,*b,maxNumIters,tolerance,&*x,verbose?&std::cout: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)
    std::cout
      << "\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) std::cout << "\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;

  try {

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

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

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

    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( "tol", &tolerance, "Relative tolerance for linear system solve." );
    clp.setOption( "max-num-iters", &maxNumIters, "Maximum of CG iterations." );

    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,verbose,tolerance,maxNumIters);
    if(!result) success = false;

    // Run using double
    result = runCgSolveExample<double>(dim,diagScale,symOp,verbose,tolerance,maxNumIters);
    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,verbose,tolerance,maxNumIters);
    if(!result) success = false;

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

#endif

#ifdef HAVE_TEUCHOS_GNU_MP

    // Run using mpf_class
    result = runCgSolveExample<mpf_class>(dim,diagScale,symOp,verbose,tolerance,maxNumIters);
    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,verbose,mpf_class(tolerance),maxNumIters);
    //if(!result) success = false;
    //The above commented-out code throws a floating-point exception?

#endif

#endif

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

  if (verbose) {
    if(success)   std::cout << "\nCongratulations! All of the tests checked out!\n";
    else          std::cout << "\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/Core/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, at the time of this writing, produces the following output:


$ ./sillyCgSolve_serial.exe

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

describe op:
  type = 'SerialTridiagLinearOp<float>', rangeDim = 500, domainDim = 500


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.03937e-05
         <= linear_properties_error_tol() = 0.0001 : passed
Warning! rel_err(sum(v4),sum(v5))
       = rel_err(-0.0412178,-0.0412174) = 1.03937e-05
         >= linear_properties_warning_tol() = 1e-06!

Checking that the adjoint agrees with the non-adjoint operator 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
         <= adjoint_error_tol() = 0.0001 : passed

Skipping check for symmetry since this->check_for_symmetry()==false

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

Testing the linear operator used with the solve ...

*** Entering LinearOpTester<float>::check(op,...) ...

describe op:
  type = 'SerialTridiagLinearOp<float>', rangeDim = 500, domainDim = 500


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.216513,0.216513) = 2.13352e-06
         <= linear_properties_error_tol() = 0.0001 : passed
Warning! rel_err(sum(v4),sum(v5))
       = rel_err(0.216513,0.216513) = 2.13352e-06
         >= linear_properties_warning_tol() = 1e-06!

Checking that the adjoint agrees with the non-adjoint operator 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(-10.5495,-10.5495) = 1.808e-07
         <= adjoint_error_tol() = 0.0001 : passed

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(18.106,18.106) = 3.1603e-07
         <= symmetry_error_tol() = 0.0001 : passed

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

Starting CG solver ...

describe A:
  type = 'SerialTridiagLinearOp<float>', rangeDim = 500, domainDim = 500

describe b:
  type = 'SerialVectorStd<float>', size = 500

describe x:
  type = 'SerialVectorStd<float>', 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.344797e-03/1.299397e+01 = 1.034939e-04 <= 10.0*tolerance = 9.999999e-04: passed

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

describe op:
  type = 'SerialTridiagLinearOp<double>', rangeDim = 500, domainDim = 500


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.394610e-01,2.394610e-01) = 1.101131e-14
         <= linear_properties_error_tol() = 1.000000e-04 : passed

Checking that the adjoint agrees with the non-adjoint operator 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(-5.913676e+00,-5.913676e+00) = 9.011435e-16
         <= adjoint_error_tol() = 1.000000e-04 : passed

Skipping check for symmetry since this->check_for_symmetry()==false

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

Testing the linear operator used with the solve ...

*** Entering LinearOpTester<double>::check(op,...) ...

describe op:
  type = 'SerialTridiagLinearOp<double>', rangeDim = 500, domainDim = 500


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.743187e-15
         <= linear_properties_error_tol() = 1.000000e-04 : passed

Checking that the adjoint agrees with the non-adjoint operator 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
         <= adjoint_error_tol() = 1.000000e-04 : passed

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.810601e+01,1.810601e+01) = 1.962174e-16
         <= symmetry_error_tol() = 1.000000e-04 : passed

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

Starting CG solver ...

describe A:
  type = 'SerialTridiagLinearOp<double>', rangeDim = 500, domainDim = 500

describe b:
  type = 'SerialVectorStd<double>', size = 500

describe x:
  type = 'SerialVectorStd<double>', 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.020000e-01 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>>::check(op,...) ...

describe op:
  type = 'SerialTridiagLinearOp<std::complex<float>>', rangeDim = 500, domainDim = 500


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((-3.456044e-01,-4.895731e-01),(-3.456052e-01,-4.895744e-01)) = 2.541641e-06
         <= linear_properties_error_tol() = 1.000000e-04 : passed
Warning! rel_err(sum(v4),sum(v5))
       = rel_err((-3.456044e-01,-4.895731e-01),(-3.456052e-01,-4.895744e-01)) = 2.541641e-06
         >= linear_properties_warning_tol() = 1.000000e-06!

Checking that the adjoint agrees with the non-adjoint operator 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.643561e+01,6.788700e+00),(1.643562e+01,6.788701e+00)) = 2.211224e-07
         <= adjoint_error_tol() = 1.000000e-04 : passed

Skipping check for symmetry since this->check_for_symmetry()==false

*** Leaving LinearOpTester<std::complex<float>>::check(...)

Testing the linear operator used with the solve ...

*** Entering LinearOpTester<std::complex<float>>::check(op,...) ...

describe op:
  type = 'SerialTridiagLinearOp<std::complex<float>>', rangeDim = 500, domainDim = 500


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.677564e-01,6.222531e-01),(-9.677570e-01,6.222509e-01)) = 1.950122e-06
         <= linear_properties_error_tol() = 1.000000e-04 : passed
Warning! rel_err(sum(v4),sum(v5))
       = rel_err((-9.677564e-01,6.222531e-01),(-9.677570e-01,6.222509e-01)) = 1.950122e-06
         >= linear_properties_warning_tol() = 1.000000e-06!

Checking that the adjoint agrees with the non-adjoint operator 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.335190e+00),(1.402399e+01,4.335191e+00)) = 3.911690e-07
         <= adjoint_error_tol() = 1.000000e-04 : passed

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.046669e+01,2.444866e+00),(1.046669e+01,2.444870e+00)) = 5.173613e-07
         <= symmetry_error_tol() = 1.000000e-04 : passed

*** Leaving LinearOpTester<std::complex<float>>::check(...)

Starting CG solver ...

describe A:
  type = 'SerialTridiagLinearOp<std::complex<float>>', rangeDim = 500, domainDim = 500

describe b:
  type = 'SerialVectorStd<std::complex<float>>', size = 500

describe x:
  type = 'SerialVectorStd<std::complex<float>>', 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.410000e-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>>::check(op,...) ...

describe op:
  type = 'SerialTridiagLinearOp<std::complex<double>>', rangeDim = 500, domainDim = 500


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((-4.443735e-01,-2.545424e-01),(-4.443735e-01,-2.545424e-01)) = 7.714438e-15
         <= linear_properties_error_tol() = 1.000000e-04 : passed

Checking that the adjoint agrees with the non-adjoint operator 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((6.554045e+00,7.729893e+00),(6.554045e+00,7.729893e+00)) = 8.314229e-16
         <= adjoint_error_tol() = 1.000000e-04 : passed

Skipping check for symmetry since this->check_for_symmetry()==false

*** Leaving LinearOpTester<std::complex<double>>::check(...)

Testing the linear operator used with the solve ...

*** Entering LinearOpTester<std::complex<double>>::check(op,...) ...

describe op:
  type = 'SerialTridiagLinearOp<std::complex<double>>', rangeDim = 500, domainDim = 500


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

Checking that the adjoint agrees with the non-adjoint operator 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
         <= adjoint_error_tol() = 1.000000e-04 : passed

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.046669e+01,2.444870e+00),(1.046669e+01,2.444870e+00)) = 1.652665e-16
         <= symmetry_error_tol() = 1.000000e-04 : passed

*** Leaving LinearOpTester<std::complex<double>>::check(...)

Starting CG solver ...

describe A:
  type = 'SerialTridiagLinearOp<std::complex<double>>', rangeDim = 500, domainDim = 500

describe b:
  type = 'SerialVectorStd<std::complex<double>>', size = 500

describe x:
  type = 'SerialVectorStd<std::complex<double>>', 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.120000e-01 sec

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

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

Testing the constructed linear operator A ...

*** Entering LinearOpTester<mpf_class>::check(op,...) ...

describe op:
  type = 'SerialTridiagLinearOp<mpf_class>', rangeDim = 500, domainDim = 500


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.554737e-01,9.554737e-01) = 3.054291e-19
         <= linear_properties_error_tol() = 1.000000e-04 : passed

Checking that the adjoint agrees with the non-adjoint operator 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(-6.665022e-01,-6.665022e-01) = 3.655873e-19
         <= adjoint_error_tol() = 1.000000e-04 : passed

Skipping check for symmetry since this->check_for_symmetry()==false

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

Testing the linear operator used with the solve ...

*** Entering LinearOpTester<mpf_class>::check(op,...) ...

describe op:
  type = 'SerialTridiagLinearOp<mpf_class>', rangeDim = 500, domainDim = 500


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(1.811994e+00,1.811994e+00) = 2.692563e-19
         <= linear_properties_error_tol() = 1.000000e-04 : passed

Checking that the adjoint agrees with the non-adjoint operator 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(-2.469042e+00,-2.469042e+00) = 3.952367e-20
         <= adjoint_error_tol() = 1.000000e-04 : passed

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(7.214884e-01,7.214884e-01) = 8.887522e-20
         <= symmetry_error_tol() = 1.000000e-04 : passed

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

Starting CG solver ...

describe A:
  type = 'SerialTridiagLinearOp<mpf_class>', rangeDim = 500, domainDim = 500

describe b:
  type = 'SerialVectorStd<mpf_class>', size = 500

describe x:
  type = 'SerialVectorStd<mpf_class>', size = 500

Iter = 0, ||b-A*x||/||b-A*x0|| = 1.000000e+00
Iter = 31, ||b-A*x||/||b-A*x0|| = 7.662633e-01
Iter = 62, ||b-A*x||/||b-A*x0|| = 1.869085e-01
Iter = 93, ||b-A*x||/||b-A*x0|| = 4.641508e-02
Iter = 124, ||b-A*x||/||b-A*x0|| = 1.001201e-02
Iter = 155, ||b-A*x||/||b-A*x0|| = 2.792430e-03
Iter = 186, ||b-A*x||/||b-A*x0|| = 5.392830e-04
Iter = 217, ||b-A*x||/||b-A*x0|| = 1.093640e-04
Iter = 219, ||b-A*x||/||b-A*x0|| = 9.685596e-05

||b-A*x||/||b|| = 1.255016e-03/1.295755e+01 = 9.685596e-05 <= 10.0*tolerance = 1.000000e-03: passed

Total time = 3.037000e+00 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_serial.exe --help are:


$ ./sillyCgSolve_serial.exe --help
Usage: ./sillyCgSolve_serial [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)
  --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)
  --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)

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


$ ./sillyCgSolve_serial.exe --unsym-op --max-num-iters=20

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

describe op:
  type = 'SerialTridiagLinearOp<float>', rangeDim = 500, domainDim = 500


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(19.2992,19.2992) = 3.95322e-07
         <= linear_properties_error_tol() = 0.0001 : passed

Checking that the adjoint agrees with the non-adjoint operator 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.72592,-1.72592) = 9.66979e-07
         <= adjoint_error_tol() = 0.0001 : passed

Skipping check for symmetry since this->check_for_symmetry()==false

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

describe op:
  type = 'MultiplicativeLinearOp<float>', rangeDim = 500, domainDim = 500
    numOps=2
    Constituent LinearOpBase objects for M = Op[0]*...*Op[numOps-1]:
      Op[0] =
        type = 'ScaledAdjointLinearOp<float>', rangeDim = 500, domainDim = 500
          overallScalar=1
          overallTransp=CONJTRANS
          Constituent transformations:
            transp=CONJTRANS
              origOp =
                type = 'SerialTridiagLinearOp<float>', rangeDim = 500, domainDim = 500
      Op[1] =
        type = 'SerialTridiagLinearOp<float>', rangeDim = 500, domainDim = 500


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(35.3825,35.3825) = 2.15626e-07
         <= linear_properties_error_tol() = 0.0001 : passed

Checking that the adjoint agrees with the non-adjoint operator 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(-36.7441,-36.744) = 2.07636e-07
         <= adjoint_error_tol() = 0.0001 : passed

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(37.9726,37.9726) = 0
         <= symmetry_error_tol() = 0.0001 : passed

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

Starting CG solver ...

describe A:
  type = 'MultiplicativeLinearOp<float>', rangeDim = 500, domainDim = 500
    numOps=2
    Constituent LinearOpBase objects for M = Op[0]*...*Op[numOps-1]:
      Op[0] =
        type = 'ScaledAdjointLinearOp<float>', rangeDim = 500, domainDim = 500
          overallScalar=1.000000e+00
          overallTransp=CONJTRANS
          Constituent transformations:
            transp=CONJTRANS
              origOp =
                type = 'SerialTridiagLinearOp<float>', rangeDim = 500, domainDim = 500
      Op[1] =
        type = 'SerialTridiagLinearOp<float>', rangeDim = 500, domainDim = 500

describe b:
  type = 'SerialVectorStd<float>', size = 500

describe x:
  type = 'SerialVectorStd<float>', size = 500

Iter = 0, ||b-A*x||/||b-A*x0|| = 1.000000e+00
Iter = 3, ||b-A*x||/||b-A*x0|| = 6.660668e-03
Iter = 6, ||b-A*x||/||b-A*x0|| = 3.351218e-05

||b-A*x||/||b|| = 1.069767e-03/3.191998e+01 = 3.351403e-05 <= 10.0*tolerance = 9.999999e-04: passed

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

describe op:
  type = 'SerialTridiagLinearOp<double>', rangeDim = 500, domainDim = 500


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(-1.680574e+00,-1.680574e+00) = 1.321243e-16
         <= linear_properties_error_tol() = 1.000000e-04 : passed

Checking that the adjoint agrees with the non-adjoint operator 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(-6.383695e-01,-6.383695e-01) = 2.260901e-15
         <= adjoint_error_tol() = 1.000000e-04 : passed

Skipping check for symmetry since this->check_for_symmetry()==false

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

describe op:
  type = 'MultiplicativeLinearOp<double>', rangeDim = 500, domainDim = 500
    numOps=2
    Constituent LinearOpBase objects for M = Op[0]*...*Op[numOps-1]:
      Op[0] =
        type = 'ScaledAdjointLinearOp<double>', rangeDim = 500, domainDim = 500
          overallScalar=1.000000e+00
          overallTransp=CONJTRANS
          Constituent transformations:
            transp=CONJTRANS
              origOp =
                type = 'SerialTridiagLinearOp<double>', rangeDim = 500, domainDim = 500
      Op[1] =
        type = 'SerialTridiagLinearOp<double>', rangeDim = 500, domainDim = 500


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(6.938410e+00,6.938410e+00) = 1.024071e-15
         <= linear_properties_error_tol() = 1.000000e-04 : passed

Checking that the adjoint agrees with the non-adjoint operator 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(-2.814797e+01,-2.814797e+01) = 2.524312e-16
         <= adjoint_error_tol() = 1.000000e-04 : passed

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(4.438044e+01,4.438044e+01) = 1.601027e-16
         <= symmetry_error_tol() = 1.000000e-04 : passed

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

Starting CG solver ...

describe A:
  type = 'MultiplicativeLinearOp<double>', rangeDim = 500, domainDim = 500
    numOps=2
    Constituent LinearOpBase objects for M = Op[0]*...*Op[numOps-1]:
      Op[0] =
        type = 'ScaledAdjointLinearOp<double>', rangeDim = 500, domainDim = 500
          overallScalar=1.000000e+00
          overallTransp=CONJTRANS
          Constituent transformations:
            transp=CONJTRANS
              origOp =
                type = 'SerialTridiagLinearOp<double>', rangeDim = 500, domainDim = 500
      Op[1] =
        type = 'SerialTridiagLinearOp<double>', rangeDim = 500, domainDim = 500

describe b:
  type = 'SerialVectorStd<double>', size = 500

describe x:
  type = 'SerialVectorStd<double>', size = 500

Iter = 0, ||b-A*x||/||b-A*x0|| = 1.000000e+00
Iter = 3, ||b-A*x||/||b-A*x0|| = 1.233805e-01
Iter = 6, ||b-A*x||/||b-A*x0|| = 2.217761e-02
Iter = 9, ||b-A*x||/||b-A*x0|| = 4.338932e-03
Iter = 12, ||b-A*x||/||b-A*x0|| = 7.668738e-04
Iter = 15, ||b-A*x||/||b-A*x0|| = 1.317471e-04
Iter = 16, ||b-A*x||/||b-A*x0|| = 7.581160e-05

||b-A*x||/||b|| = 2.307234e-03/3.043379e+01 = 7.581160e-05 <= 10.0*tolerance = 1.000000e-03: passed

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

describe op:
  type = 'SerialTridiagLinearOp<std::complex<float>>', rangeDim = 500, domainDim = 500


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.723428e-01,2.028702e+00),(2.723447e-01,2.028702e+00)) = 9.390742e-07
         <= linear_properties_error_tol() = 1.000000e-04 : passed

Checking that the adjoint agrees with the non-adjoint operator 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.121651e+01,3.660773e+00),(1.121650e+01,3.660776e+00)) = 2.914302e-07
         <= adjoint_error_tol() = 1.000000e-04 : passed

Skipping check for symmetry since this->check_for_symmetry()==false

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

describe op:
  type = 'MultiplicativeLinearOp<std::complex<float>>', rangeDim = 500, domainDim = 500
    numOps=2
    Constituent LinearOpBase objects for M = Op[0]*...*Op[numOps-1]:
      Op[0] =
        type = 'ScaledAdjointLinearOp<std::complex<float>>', rangeDim = 500, domainDim = 500
          overallScalar=(1.000000e+00,0.000000e+00)
          overallTransp=CONJTRANS
          Constituent transformations:
            transp=CONJTRANS
              origOp =
                type = 'SerialTridiagLinearOp<std::complex<float>>', rangeDim = 500, domainDim = 500
      Op[1] =
        type = 'SerialTridiagLinearOp<std::complex<float>>', rangeDim = 500, domainDim = 500


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((-1.574269e+01,1.134796e+01),(-1.574271e+01,1.134796e+01)) = 8.368621e-07
         <= linear_properties_error_tol() = 1.000000e-04 : passed

Checking that the adjoint agrees with the non-adjoint operator 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((5.278718e+01,4.448534e+00),(5.278718e+01,4.448544e+00)) = 1.938938e-07
         <= adjoint_error_tol() = 1.000000e-04 : passed

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((2.419142e+01,1.037749e+01),(2.419140e+01,1.037749e+01)) = 6.680345e-07
         <= symmetry_error_tol() = 1.000000e-04 : passed

*** Leaving LinearOpTester<std::complex<float>>::check(...)

Starting CG solver ...

describe A:
  type = 'MultiplicativeLinearOp<std::complex<float>>', rangeDim = 500, domainDim = 500
    numOps=2
    Constituent LinearOpBase objects for M = Op[0]*...*Op[numOps-1]:
      Op[0] =
        type = 'ScaledAdjointLinearOp<std::complex<float>>', rangeDim = 500, domainDim = 500
          overallScalar=(1.000000e+00,0.000000e+00)
          overallTransp=CONJTRANS
          Constituent transformations:
            transp=CONJTRANS
              origOp =
                type = 'SerialTridiagLinearOp<std::complex<float>>', rangeDim = 500, domainDim = 500
      Op[1] =
        type = 'SerialTridiagLinearOp<std::complex<float>>', rangeDim = 500, domainDim = 500

describe b:
  type = 'SerialVectorStd<std::complex<float>>', size = 500

describe x:
  type = 'SerialVectorStd<std::complex<float>>', size = 500

Iter = 0, ||b-A*x||/||b-A*x0|| = 1.000000e+00
Iter = 3, ||b-A*x||/||b-A*x0|| = 1.431026e-01
Iter = 6, ||b-A*x||/||b-A*x0|| = 2.387369e-02
Iter = 9, ||b-A*x||/||b-A*x0|| = 3.633558e-03
Iter = 12, ||b-A*x||/||b-A*x0|| = 5.538159e-04
Iter = 15, ||b-A*x||/||b-A*x0|| = 8.121288e-05

||b-A*x||/||b|| = 3.305705e-03/4.070510e+01 = 8.121108e-05 <= 10.0*tolerance = 9.999999e-04: passed

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

describe op:
  type = 'SerialTridiagLinearOp<std::complex<double>>', rangeDim = 500, domainDim = 500


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((-1.970843e+00,-3.441251e+01),(-1.970843e+00,-3.441251e+01)) = 6.291975e-16
         <= linear_properties_error_tol() = 1.000000e-04 : passed

Checking that the adjoint agrees with the non-adjoint operator 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.372625e+00,1.850210e+01),(1.372625e+00,1.850210e+01)) = 8.360598e-16
         <= adjoint_error_tol() = 1.000000e-04 : passed

Skipping check for symmetry since this->check_for_symmetry()==false

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

describe op:
  type = 'MultiplicativeLinearOp<std::complex<double>>', rangeDim = 500, domainDim = 500
    numOps=2
    Constituent LinearOpBase objects for M = Op[0]*...*Op[numOps-1]:
      Op[0] =
        type = 'ScaledAdjointLinearOp<std::complex<double>>', rangeDim = 500, domainDim = 500
          overallScalar=(1.000000e+00,0.000000e+00)
          overallTransp=CONJTRANS
          Constituent transformations:
            transp=CONJTRANS
              origOp =
                type = 'SerialTridiagLinearOp<std::complex<double>>', rangeDim = 500, domainDim = 500
      Op[1] =
        type = 'SerialTridiagLinearOp<std::complex<double>>', rangeDim = 500, domainDim = 500


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.472665e+01,3.895074e+01),(-5.472665e+01,3.895074e+01)) = 9.520071e-16
         <= linear_properties_error_tol() = 1.000000e-04 : passed

Checking that the adjoint agrees with the non-adjoint operator 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.017174e+01,1.130154e+01),(-1.017174e+01,1.130154e+01)) = 1.263687e-15
         <= adjoint_error_tol() = 1.000000e-04 : passed

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((2.205033e+01,-4.284830e+01),(2.205033e+01,-4.284830e+01)) = 1.823840e-15
         <= symmetry_error_tol() = 1.000000e-04 : passed

*** Leaving LinearOpTester<std::complex<double>>::check(...)

Starting CG solver ...

describe A:
  type = 'MultiplicativeLinearOp<std::complex<double>>', rangeDim = 500, domainDim = 500
    numOps=2
    Constituent LinearOpBase objects for M = Op[0]*...*Op[numOps-1]:
      Op[0] =
        type = 'ScaledAdjointLinearOp<std::complex<double>>', rangeDim = 500, domainDim = 500
          overallScalar=(1.000000e+00,0.000000e+00)
          overallTransp=CONJTRANS
          Constituent transformations:
            transp=CONJTRANS
              origOp =
                type = 'SerialTridiagLinearOp<std::complex<double>>', rangeDim = 500, domainDim = 500
      Op[1] =
        type = 'SerialTridiagLinearOp<std::complex<double>>', rangeDim = 500, domainDim = 500

describe b:
  type = 'SerialVectorStd<std::complex<double>>', size = 500

describe x:
  type = 'SerialVectorStd<std::complex<double>>', size = 500

Iter = 0, ||b-A*x||/||b-A*x0|| = 1.000000e+00
Iter = 3, ||b-A*x||/||b-A*x0|| = 2.153769e-02
Iter = 6, ||b-A*x||/||b-A*x0|| = 3.394906e-04
Iter = 7, ||b-A*x||/||b-A*x0|| = 9.023713e-05

||b-A*x||/||b|| = 4.053765e-03/4.492347e+01 = 9.023713e-05 <= 10.0*tolerance = 1.000000e-03: passed

Total time = 5.900000e-02 sec

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

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

Testing the constructed linear operator A ...

*** Entering LinearOpTester<mpf_class>::check(op,...) ...

describe op:
  type = 'SerialTridiagLinearOp<mpf_class>', rangeDim = 500, domainDim = 500


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.352466e+02,2.352466e+02) = 1.221329e-20
         <= linear_properties_error_tol() = 1.000000e-04 : passed

Checking that the adjoint agrees with the non-adjoint operator 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(6.045283e+01,6.045283e+01) = 1.019303e-20
         <= adjoint_error_tol() = 1.000000e-04 : passed

Skipping check for symmetry since this->check_for_symmetry()==false

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

describe op:
  type = 'MultiplicativeLinearOp<mpf_class>', rangeDim = 500, domainDim = 500
    numOps=2
    Constituent LinearOpBase objects for M = Op[0]*...*Op[numOps-1]:
      Op[0] =
        type = 'ScaledAdjointLinearOp<mpf_class>', rangeDim = 500, domainDim = 500
          overallScalar=1.000000e+00
          overallTransp=CONJTRANS
          Constituent transformations:
            transp=CONJTRANS
              origOp =
                type = 'SerialTridiagLinearOp<mpf_class>', rangeDim = 500, domainDim = 500
      Op[1] =
        type = 'SerialTridiagLinearOp<mpf_class>', rangeDim = 500, domainDim = 500


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.173757e+02,2.173757e+02) = 1.945199e-20
         <= linear_properties_error_tol() = 1.000000e-04 : passed

Checking that the adjoint agrees with the non-adjoint operator 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(4.949277e+01,4.949277e+01) = 1.752502e-20
         <= adjoint_error_tol() = 1.000000e-04 : passed

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(5.540864e+01,5.540864e+01) = 2.152412e-20
         <= symmetry_error_tol() = 1.000000e-04 : passed

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

Starting CG solver ...

describe A:
  type = 'MultiplicativeLinearOp<mpf_class>', rangeDim = 500, domainDim = 500
    numOps=2
    Constituent LinearOpBase objects for M = Op[0]*...*Op[numOps-1]:
      Op[0] =
        type = 'ScaledAdjointLinearOp<mpf_class>', rangeDim = 500, domainDim = 500
          overallScalar=1.000000e+00
          overallTransp=CONJTRANS
          Constituent transformations:
            transp=CONJTRANS
              origOp =
                type = 'SerialTridiagLinearOp<mpf_class>', rangeDim = 500, domainDim = 500
      Op[1] =
        type = 'SerialTridiagLinearOp<mpf_class>', rangeDim = 500, domainDim = 500

describe b:
  type = 'SerialVectorStd<mpf_class>', size = 500

describe x:
  type = 'SerialVectorStd<mpf_class>', size = 500

Iter = 0, ||b-A*x||/||b-A*x0|| = 1.000000e+00
Iter = 3, ||b-A*x||/||b-A*x0|| = 2.692000e-01
Iter = 6, ||b-A*x||/||b-A*x0|| = 2.507806e-02
Iter = 9, ||b-A*x||/||b-A*x0|| = 3.479960e-03
Iter = 12, ||b-A*x||/||b-A*x0|| = 5.338185e-04
Iter = 15, ||b-A*x||/||b-A*x0|| = 7.988701e-05

||b-A*x||/||b|| = 1.412440e-03/1.768048e+01 = 7.988701e-05 <= 10.0*tolerance = 1.000000e-03: passed

Total time = 5.060000e-01 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::scale() and Thyra::multiply() which create implicit Thyra::ScaledAdjointedLinearOp and Thyra::MultiplicativeLinearOp objects.

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

Generated on Thu Sep 18 12:39:53 2008 for Thyra ANA Operator/VectorBase Interfaces and Related Software by doxygen 1.3.9.1