![\[ A= \left[\begin{array}{rrrrrrrrrr} 2 & -1 \\ -1 & 2 & -1 \\ & \ddots & \ddots & \ddots \\ & & -1 & 2 & -1 \\ & & & -1 & 2 \end{array}\right]. \]](form_13.png)
SerialTridiagLinearOp with the example linear ANA implementation sillyPowerMethod().
More...Classes | |
| class | SerialTridiagLinearOp< Scalar > |
| Simple example subclass for serial tridiagonal matrices. More... | |
SerialTridiagLinearOp with the example linear ANA implementation sillyPowerMethod().
This example program is contained in the source file:
./example/Core/sillyPowerMethod_serial.cpp
where ./ is the base source directory for Thyra (i.e. ???/Trilinos/packages/Thyra).
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 runPowerMethodExample() below) that calls sillyPowerMethod(). In this example program, the matrix constructed and used is the well-known tridiagonal matrix
The power method is then run on the matrix 
run for a number of iterations (or until convergence to some tolerance).
After this, the first diagonal element 
is then scaled to 
and the power method is run again (which much faster convergence).
The following templated function implements the example described above:
template<class Scalar> bool runPowerMethodExample( const int dim ,const int maxNumIters ,const bool verbose ,const bool dumpAll ) { typedef Teuchos::ScalarTraits<Scalar> ST; typedef typename ST::magnitudeType ScalarMag; bool success = true; bool result; if(verbose) std::cout << "\n***\n*** Running power method example using scalar type = \'" << ST::name() << "\' ...\n***\n" << std::scientific; // // (1) Setup the initial tridiagonal operator // // [ 2 -1 ] // [ -1 2 -1 ] // A = [ . . . ] // [ -1 2 -1 ] // [ -1 2 ] // if(verbose) std::cout << "\n(1) Constructing tridiagonal matrix A of dimension = " << dim << " ...\n"; std::vector<Scalar> lower(dim-1), diag(dim), upper(dim-1); const Scalar one = ST::one(), two = Scalar(2)*one; int k = 0; diag[k] = two; upper[k] = -one; // First row for( k = 1; k < dim - 1; ++k ) { lower[k-1] = -one; diag[k] = two; upper[k] = -one; // Middle rows } lower[k-1] = -one; diag[k] = two; // Last row Teuchos::RefCountPtr<SerialTridiagLinearOp<Scalar> > A = Teuchos::rcp( new SerialTridiagLinearOp<Scalar>(dim,&lower[0],&diag[0],&upper[0]) ); if( verbose && dumpAll ) std::cout << "\nA =\n" << *A; // // (2) Run the power method ANA // if(verbose) std::cout << "\n(2) Running the power method on matrix A ...\n"; Scalar lambda = ST::zero(); ScalarMag tolerance = ScalarMag(1e-3)*Teuchos::ScalarTraits<ScalarMag>::one(); result = sillyPowerMethod(*A,maxNumIters,tolerance,&lambda,(verbose?&std::cout:NULL)); if(!result) success = false; if(verbose) std::cout << "\n Estimate of dominate eigenvalue lambda = " << lambda << std::endl; // // (3) Increase dominance of first eigenvalue // if(verbose) std::cout << "\n(3) Increasing first diagonal entry by factor of 10 ...\n"; diag[0] *= 10; A->initialize(dim,&lower[0],&diag[0],&upper[0]); if( verbose && dumpAll ) std::cout << "A =\n" << *A; // // (4) Run the power method ANA // if(verbose) std::cout << "\n(4) Running the power method again on matrix A ...\n"; result = sillyPowerMethod(*A,maxNumIters,tolerance,&lambda,(verbose?&std::cout:NULL)); if(!result) success = false; if(verbose) std::cout << "\n Estimate of dominate eigenvalue lambda = " << lambda << std::endl; return success; } // end runPowerMethodExample()
The above templated function runPowerMethodExample() is instantiated with the following scalar types:
float and double std::complex<float> and std::complex<double> (if --enable-teuchos-complex was used at configuration time) mpf_class (if --enable-teuchos-gmp was used at configuration time) std::complex<mpf_class> (if --enable-teuchos-complex and --enable-teuchos-gmp where used at configuration time)
and is called multiple times from within the following main() program function:
int main(int argc, char *argv[]) { using Teuchos::CommandLineProcessor; bool success = true; bool verbose = true; bool result; try { // // Read in command-line options // int dim = 4; bool dumpAll = false; CommandLineProcessor clp(false); // Don't throw exceptions clp.setOption( "verbose", "quiet", &verbose, "Determines if any output is printed or not." ); clp.setOption( "dim", &dim, "Dimension of the linear system." ); clp.setOption( "dump-all", "no-dump", &dumpAll, "Determines if quantities are dumped or not." ); 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!" ); int maxNumIters = 10*dim; // Run using float result = runPowerMethodExample<float>(dim,maxNumIters,verbose,dumpAll); if(!result) success = false; // Run using double result = runPowerMethodExample<double>(dim,maxNumIters,verbose,dumpAll); if(!result) success = false; #if defined(HAVE_COMPLEX) && defined(HAVE_TEUCHOS_COMPLEX) // Run using std::complex<float> result = runPowerMethodExample<std::complex<float> >(dim,maxNumIters,verbose,dumpAll); if(!result) success = false; // Run using std::complex<double> result = runPowerMethodExample<std::complex<double> >(dim,maxNumIters,verbose,dumpAll); if(!result) success = false; #endif #ifdef HAVE_TEUCHOS_GNU_MP // Run using mpf_class result = runPowerMethodExample<mpf_class >(dim,maxNumIters,verbose,dumpAll); if(!result) success = false; #if defined(HAVE_COMPLEX) && defined(HAVE_TEUCHOS_COMPLEX) // Run using std::complex<mpf_class> //result = runPowerMethodExample<std::complex<mpf_class> >(dim,maxNumIters,verbose,dumpAll); //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/sillyPowerMethod_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:
$ ./sillyPowerMethod_serial.exe *** *** Running power method example using scalar type = 'float' ... *** (1) Constructing tridiagonal matrix A of dimension = 4 ... (2) Running the power method on matrix A ... Starting power method ... Iter = 0, lambda = 2.634806e+00, ||A*q-lambda*q|| = 9.051834e-01 Iter = 4, lambda = 3.552294e+00, ||A*q-lambda*q|| = 2.496006e-01 Iter = 8, lambda = 3.612831e+00, ||A*q-lambda*q|| = 7.194513e-02 Iter = 12, lambda = 3.617640e+00, ||A*q-lambda*q|| = 1.981858e-02 Iter = 16, lambda = 3.618005e+00, ||A*q-lambda*q|| = 5.435463e-03 Iter = 20, lambda = 3.618032e+00, ||A*q-lambda*q|| = 1.490343e-03 Iter = 24, lambda = 3.618034e+00, ||A*q-lambda*q|| = 4.085834e-04 Estimate of dominate eigenvalue lambda = 3.618034e+00 (3) Increasing first diagonal entry by factor of 10 ... (4) Running the power method again on matrix A ... Starting power method ... Iter = 0, lambda = 1.791493e+01, ||A*q-lambda*q|| = 5.881544e+00 Iter = 4, lambda = 2.005555e+01, ||A*q-lambda*q|| = 3.251727e-03 Iter = 8, lambda = 2.005556e+01, ||A*q-lambda*q|| = 2.788594e-06 Estimate of dominate eigenvalue lambda = 2.005556e+01 *** *** Running power method example using scalar type = 'double' ... *** (1) Constructing tridiagonal matrix A of dimension = 4 ... (2) Running the power method on matrix A ... Starting power method ... Iter = 0, lambda = 2.634806e+00, ||A*q-lambda*q|| = 9.051834e-01 Iter = 4, lambda = 3.552293e+00, ||A*q-lambda*q|| = 2.496007e-01 Iter = 8, lambda = 3.612831e+00, ||A*q-lambda*q|| = 7.194514e-02 Iter = 12, lambda = 3.617641e+00, ||A*q-lambda*q|| = 1.981868e-02 Iter = 16, lambda = 3.618004e+00, ||A*q-lambda*q|| = 5.435543e-03 Iter = 20, lambda = 3.618032e+00, ||A*q-lambda*q|| = 1.490271e-03 Iter = 24, lambda = 3.618034e+00, ||A*q-lambda*q|| = 4.085797e-04 Estimate of dominate eigenvalue lambda = 3.618034e+00 (3) Increasing first diagonal entry by factor of 10 ... (4) Running the power method again on matrix A ... Starting power method ... Iter = 0, lambda = 1.791493e+01, ||A*q-lambda*q|| = 5.881543e+00 Iter = 4, lambda = 2.005555e+01, ||A*q-lambda*q|| = 3.251777e-03 Iter = 8, lambda = 2.005556e+01, ||A*q-lambda*q|| = 2.681092e-06 Estimate of dominate eigenvalue lambda = 2.005556e+01 *** *** Running power method example using scalar type = 'std::complex<float>' ... *** (1) Constructing tridiagonal matrix A of dimension = 4 ... (2) Running the power method on matrix A ... Starting power method ... Iter = 0, lambda = (1.885460e+00,2.235174e-08), ||A*q-lambda*q|| = 8.366854e-01 Iter = 4, lambda = (2.908303e+00,1.490116e-08), ||A*q-lambda*q|| = 4.634469e-01 Iter = 8, lambda = (3.465825e+00,-5.960464e-08), ||A*q-lambda*q|| = 3.592391e-01 Iter = 12, lambda = (3.604720e+00,-2.980232e-08), ||A*q-lambda*q|| = 1.146174e-01 Iter = 16, lambda = (3.617021e+00,-9.778887e-09), ||A*q-lambda*q|| = 3.181574e-02 Iter = 20, lambda = (3.617958e+00,9.010546e-08), ||A*q-lambda*q|| = 8.730849e-03 Iter = 24, lambda = (3.618028e+00,1.565786e-08), ||A*q-lambda*q|| = 2.393888e-03 Iter = 28, lambda = (3.618034e+00,9.727955e-09), ||A*q-lambda*q|| = 6.563177e-04 Estimate of dominate eigenvalue lambda = (3.618034e+00,9.727955e-09) (3) Increasing first diagonal entry by factor of 10 ... (4) Running the power method again on matrix A ... Starting power method ... Iter = 0, lambda = (1.316405e+01,-1.341105e-07), ||A*q-lambda*q|| = 8.828070e+00 Iter = 4, lambda = (2.005556e+01,1.859358e-07), ||A*q-lambda*q|| = 4.617708e-03 Iter = 8, lambda = (2.005555e+01,5.564779e-08), ||A*q-lambda*q|| = 4.349416e-06 Estimate of dominate eigenvalue lambda = (2.005555e+01,5.564779e-08) *** *** Running power method example using scalar type = 'std::complex<double>' ... *** (1) Constructing tridiagonal matrix A of dimension = 4 ... (2) Running the power method on matrix A ... Starting power method ... Iter = 0, lambda = (1.885459e+00,1.387779e-17), ||A*q-lambda*q|| = 8.366853e-01 Iter = 4, lambda = (2.908303e+00,-2.775558e-17), ||A*q-lambda*q|| = 4.634469e-01 Iter = 8, lambda = (3.465825e+00,-1.040834e-16), ||A*q-lambda*q|| = 3.592392e-01 Iter = 12, lambda = (3.604720e+00,1.249001e-16), ||A*q-lambda*q|| = 1.146174e-01 Iter = 16, lambda = (3.617021e+00,-1.058181e-16), ||A*q-lambda*q|| = 3.181573e-02 Iter = 20, lambda = (3.617958e+00,6.179952e-17), ||A*q-lambda*q|| = 8.730912e-03 Iter = 24, lambda = (3.618028e+00,-4.602438e-17), ||A*q-lambda*q|| = 2.393871e-03 Iter = 28, lambda = (3.618034e+00,-8.632960e-18), ||A*q-lambda*q|| = 6.563170e-04 Estimate of dominate eigenvalue lambda = (3.618034e+00,-8.632960e-18) (3) Increasing first diagonal entry by factor of 10 ... (4) Running the power method again on matrix A ... Starting power method ... Iter = 0, lambda = (1.316405e+01,-5.551115e-16), ||A*q-lambda*q|| = 8.828069e+00 Iter = 4, lambda = (2.005555e+01,-3.814166e-16), ||A*q-lambda*q|| = 4.617726e-03 Iter = 8, lambda = (2.005556e+01,-1.371085e-17), ||A*q-lambda*q|| = 3.719354e-06 Estimate of dominate eigenvalue lambda = (2.005556e+01,-1.371085e-17) *** *** Running power method example using scalar type = 'mpf_class' ... *** (1) Constructing tridiagonal matrix A of dimension = 4 ... (2) Running the power method on matrix A ... Starting power method ... Iter = 0, lambda = 9.402176e-01, ||A*q-lambda*q|| = 1.103226e+00 Iter = 4, lambda = 3.577075e+00, ||A*q-lambda*q|| = 1.999005e-01 Iter = 8, lambda = 3.614877e+00, ||A*q-lambda*q|| = 5.609685e-02 Iter = 12, lambda = 3.617796e+00, ||A*q-lambda*q|| = 1.542336e-02 Iter = 16, lambda = 3.618016e+00, ||A*q-lambda*q|| = 4.229461e-03 Iter = 20, lambda = 3.618033e+00, ||A*q-lambda*q|| = 1.159585e-03 Iter = 24, lambda = 3.618034e+00, ||A*q-lambda*q|| = 3.179170e-04 Estimate of dominate eigenvalue lambda = 3.618034e+00 (3) Increasing first diagonal entry by factor of 10 ... (4) Running the power method again on matrix A ... Starting power method ... Iter = 0, lambda = 1.102725e+01, ||A*q-lambda*q|| = 9.675451e+00 Iter = 4, lambda = 2.005556e+01, ||A*q-lambda*q|| = 1.920248e-03 Iter = 8, lambda = 2.005556e+01, ||A*q-lambda*q|| = 1.572132e-06 Estimate of dominate eigenvalue lambda = 2.005556e+01 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 ./sillyPowerMethod_serial.exe --help are:
$ ./sillyPowerMethod_serial.exe --help
Usage: ./sillyPowerMethod_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=4)
--dump-all bool Determines if quantities are dumped or not.
--no-dump (default: --no-dump)
To see the full listing of this example program click: sillyPowerMethod_serial.cpp
1.3.9.1