Tpetra Matrix/Vector Services Version of the Day
Tpetra_Power_Method.hpp
00001 #ifndef TPETRA_POWER_METHOD_HPP
00002 #define TPETRA_POWER_METHOD_HPP
00003 
00004 #include <Tpetra_Operator.hpp>
00005 #include <Tpetra_Vector.hpp>
00006 #include <Teuchos_ScalarTraits.hpp>
00007 
00008 namespace TpetraExamples {
00009 
00012   template <class Scalar, class Ordinal>
00013   Scalar powerMethod(const Teuchos::RCP<const Tpetra::Operator<Scalar,Ordinal> > &A, int niters, typename Teuchos::ScalarTraits<Scalar>::magnitudeType tolerance, bool verbose) 
00014   {
00015     if ( A->getRangeMap() != A->getDomainMap() ) throw std::runtime_error("TpetraExamples::powerMethod(): operator must have domain and range maps that are equivalent.");
00016     typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType Magnitude;
00017     const bool NO_INITIALIZE_TO_ZERO = false;
00018     // create three vectors; do not bother initializing q to zero, as we will fill it with random below
00019     Tpetra::Vector<Scalar,Ordinal> z(A->getRangeMap(), NO_INITIALIZE_TO_ZERO),
00020                                    q(A->getRangeMap(), NO_INITIALIZE_TO_ZERO),
00021                                    r(A->getRangeMap(), NO_INITIALIZE_TO_ZERO);
00022     // Fill z with random numbers
00023     z.randomize();
00024     // Variables needed for iteration
00025     const Scalar ONE  = Teuchos::ScalarTraits<Scalar>::one();
00026     const Scalar ZERO = Teuchos::ScalarTraits<Scalar>::zero();
00027     Scalar lambda = static_cast<Scalar>(0.0);
00028     Magnitude normz, residual = static_cast<Magnitude>(0.0);
00029     // power iteration
00030     for (int iter = 0; iter < niters; ++iter) {
00031       normz = z.norm2();                            // Compute 2-norm of z
00032       q.scale(ONE/normz, z);                        // Set q = z / normz
00033       A->apply(q, z);                               // Compute z = A*q
00034       lambda = q.dot(z);                            // Approximate maximum eigenvalue: lamba = dot(q,z)
00035       if ( iter % 100 == 0 || iter + 1 == niters ) {
00036         r.update(ONE, z, -lambda, q, ZERO);     // Compute A*q - lambda*q
00037         residual = Teuchos::ScalarTraits<Scalar>::magnitude(r.norm2() / lambda);
00038         if (verbose) {
00039           std::cout << "Iter = " << iter << "  Lambda = " << lambda 
00040                     << "  Residual of A*q - lambda*q = " 
00041                     << residual << std::endl;
00042         }
00043       } 
00044       if (residual < tolerance) {
00045         break;
00046       }
00047     }
00048     return lambda;
00049   }
00050 
00051 } // end of namespace TpetraExamples
00052 
00061 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines