![\[ A= \left[\begin{array}{rrrrrrrrrr} 2 & -1 \\ -1 & 2 & -1 \\ & \ddots & \ddots & \ddots \\ & & -1 & 2 & -1 \\ & & & -1 & 2 \end{array}\right]. \]](form_17.png)
ExampleTridiagSerialLinearOp with the example linear ANA implementation sillyPowerMethod().
More...Classes | |
| class | ExampleTridiagSerialLinearOp< Scalar > |
| Simple example subclass for serial tridiagonal matrices. More... | |
ExampleTridiagSerialLinearOp with the example linear ANA implementation sillyPowerMethod().
This example program is contained in the source file:
./example/operator_vector/sillyPowerMethod_serial.cpp
where ./ is the base source directory for Thyra (i.e. ???/Trilinos/packages/Thyra).
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 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 ) { using Teuchos::OSTab; typedef Teuchos::ScalarTraits<Scalar> ST; 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 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) *out << "\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<ExampleTridiagSerialLinearOp<Scalar> > A = Teuchos::rcp( new ExampleTridiagSerialLinearOp<Scalar>(dim,&lower[0],&diag[0],&upper[0]) ); if( verbose && dumpAll ) *out << "\nA =\n" << *A; // // (2) Run the power method ANA // if(verbose) *out << "\n(2) Running the power method on matrix A ...\n"; Scalar lambda = ST::zero(); ScalarMag tolerance = ScalarMag(1e-3)*Teuchos::ScalarTraits<ScalarMag>::one(); if(1){ OSTab tab(out); result = sillyPowerMethod(*A,maxNumIters,tolerance,&lambda,out.get()); if(!result) success = false; if(verbose) *out << "\nEstimate of dominate eigenvalue lambda = " << lambda << std::endl; } // // (3) Increase dominance of first eigenvalue // if(verbose) *out << "\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 ) *out << "A =\n" << *A; // // (4) Run the power method ANA // if(verbose) *out << "\n(4) Running the power method again on matrix A ...\n"; if(1){ OSTab tab(out); result = sillyPowerMethod(*A,maxNumIters,tolerance,&lambda,out.get()); if(!result) success = false; if(verbose) *out << "\nEstimate 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; Teuchos::RefCountPtr<Teuchos::FancyOStream> out = Teuchos::VerboseObjectBase::getDefaultOStream(); try { // // Read in command-line options // int dim = 4; bool dumpAll = false; CommandLineProcessor clp; clp.throwExceptions(false); clp.addOutputSetupOptions(true); clp.setOption( "verbose", "quiet", &verbose, "Determines if any output is printed or not." ); clp.setOption( "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 } 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/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 produces the following output:
*** *** 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 (target tolerrance = 1.000000e-03) ... Iter = 0, lambda = 5.745785e-01, ||A*q-lambda*q|| = 4.269211e-01 Iter = 4, lambda = 3.301265e+00, ||A*q-lambda*q|| = 5.727251e-01 Iter = 8, lambda = 3.595284e+00, ||A*q-lambda*q|| = 1.493326e-01 Iter = 12, lambda = 3.616292e+00, ||A*q-lambda*q|| = 4.170860e-02 Iter = 16, lambda = 3.617903e+00, ||A*q-lambda*q|| = 1.145323e-02 Iter = 20, lambda = 3.618024e+00, ||A*q-lambda*q|| = 3.140391e-03 Iter = 24, lambda = 3.618033e+00, ||A*q-lambda*q|| = 8.610450e-04 Estimate of dominate eigenvalue lambda = 3.618033e+00 (3) Increasing first diagonal entry by factor of 10 ... (4) Running the power method again on matrix A ... Starting power method (target tolerrance = 1.000000e-03) ... Iter = 0, lambda = 8.132596e+00, ||A*q-lambda*q|| = 9.248439e+00 Iter = 4, lambda = 2.005555e+01, ||A*q-lambda*q|| = 2.941549e-03 Iter = 8, lambda = 2.005556e+01, ||A*q-lambda*q|| = 2.368338e-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 (target tolerrance = 1.000000e-03) ... Iter = 0, lambda = 5.745785e-01, ||A*q-lambda*q|| = 4.269210e-01 Iter = 4, lambda = 3.301264e+00, ||A*q-lambda*q|| = 5.727251e-01 Iter = 8, lambda = 3.595283e+00, ||A*q-lambda*q|| = 1.493328e-01 Iter = 12, lambda = 3.616291e+00, ||A*q-lambda*q|| = 4.170859e-02 Iter = 16, lambda = 3.617903e+00, ||A*q-lambda*q|| = 1.145327e-02 Iter = 20, lambda = 3.618024e+00, ||A*q-lambda*q|| = 3.140456e-03 Iter = 24, lambda = 3.618033e+00, ||A*q-lambda*q|| = 8.610081e-04 Estimate of dominate eigenvalue lambda = 3.618033e+00 (3) Increasing first diagonal entry by factor of 10 ... (4) Running the power method again on matrix A ... Starting power method (target tolerrance = 1.000000e-03) ... Iter = 0, lambda = 8.132596e+00, ||A*q-lambda*q|| = 9.248439e+00 Iter = 4, lambda = 2.005556e+01, ||A*q-lambda*q|| = 2.941630e-03 Iter = 8, lambda = 2.005556e+01, ||A*q-lambda*q|| = 2.259420e-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 (target tolerrance = 1.000000e-03) ... Iter = 0, lambda = (1.617793e+00,1.490116e-08), ||A*q-lambda*q|| = 1.333844e+00 Iter = 4, lambda = (3.550094e+00,-5.820766e-09), ||A*q-lambda*q|| = 2.528435e-01 Iter = 8, lambda = (3.612626e+00,7.217750e-09), ||A*q-lambda*q|| = 7.334734e-02 Iter = 12, lambda = (3.617626e+00,6.402843e-10), ||A*q-lambda*q|| = 2.020937e-02 Iter = 16, lambda = (3.618003e+00,-1.698209e-08), ||A*q-lambda*q|| = 5.542850e-03 Iter = 20, lambda = (3.618032e+00,-8.076313e-09), ||A*q-lambda*q|| = 1.519599e-03 Iter = 24, lambda = (3.618034e+00,-1.025001e-08), ||A*q-lambda*q|| = 4.167181e-04 Estimate of dominate eigenvalue lambda = (3.618034e+00,-1.025001e-08) (3) Increasing first diagonal entry by factor of 10 ... (4) Running the power method again on matrix A ... Starting power method (target tolerrance = 1.000000e-03) ... Iter = 0, lambda = (6.482740e+00,-1.490116e-08), ||A*q-lambda*q|| = 7.724599e+00 Iter = 4, lambda = (2.005554e+01,-1.837986e-08), ||A*q-lambda*q|| = 1.717740e-02 Iter = 8, lambda = (2.005555e+01,-3.240486e-07), ||A*q-lambda*q|| = 1.432346e-05 Estimate of dominate eigenvalue lambda = (2.005555e+01,-3.240486e-07) *** *** 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 (target tolerrance = 1.000000e-03) ... Iter = 0, lambda = (1.617793e+00,2.775558e-17), ||A*q-lambda*q|| = 1.333844e+00 Iter = 4, lambda = (3.550094e+00,1.040834e-17), ||A*q-lambda*q|| = 2.528435e-01 Iter = 8, lambda = (3.612625e+00,1.040834e-17), ||A*q-lambda*q|| = 7.334728e-02 Iter = 12, lambda = (3.617625e+00,-1.116728e-17), ||A*q-lambda*q|| = 2.020930e-02 Iter = 16, lambda = (3.618003e+00,1.081492e-17), ||A*q-lambda*q|| = 5.542757e-03 Iter = 20, lambda = (3.618032e+00,-3.577867e-18), ||A*q-lambda*q|| = 1.519668e-03 Iter = 24, lambda = (3.618034e+00,9.903509e-18), ||A*q-lambda*q|| = 4.166393e-04 Estimate of dominate eigenvalue lambda = (3.618034e+00,9.903509e-18) (3) Increasing first diagonal entry by factor of 10 ... (4) Running the power method again on matrix A ... Starting power method (target tolerrance = 1.000000e-03) ... Iter = 0, lambda = (6.482739e+00,1.942890e-16), ||A*q-lambda*q|| = 7.724598e+00 Iter = 4, lambda = (2.005554e+01,6.367994e-17), ||A*q-lambda*q|| = 1.717743e-02 Iter = 8, lambda = (2.005556e+01,3.011987e-16), ||A*q-lambda*q|| = 1.414814e-05 Estimate of dominate eigenvalue lambda = (2.005556e+01,3.011987e-16) 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 returned from ./sillyPowerMethod_serial.exe --echo-command-line --help are:
Echoing the command-line:
../example/operator_vector/sillyPowerMethod_serial.exe --echo-command-line --help
Usage: ../example/operator_vector/sillyPowerMethod_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=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