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     typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType Magnitude;
00016     typedef Tpetra::Vector<Scalar,Ordinal> Vector;
00017
00018     if ( A->getRangeMap() != A->getDomainMap() ) {
00019       throw std::runtime_error("TpetraExamples::powerMethod(): operator must have domain and range maps that are equivalent.");
00020     }
00021     // create three vectors, fill z with random numbers
00022     Teuchos::RCP<Vector> z, q, r;
00023     q = Tpetra::createVector<Scalar>(A->getRangeMap());
00024     r = Tpetra::createVector<Scalar>(A->getRangeMap());
00025     z = Tpetra::createVector<Scalar>(A->getRangeMap());
00026     z->randomize();
00027     //
00028     Scalar lambda = 0.0;
00029     Magnitude normz, residual = 0.0;
00030     // power iteration
00031     for (int iter = 0; iter < niters; ++iter) {
00032       normz = z->norm2();                             // Compute 2-norm of z
00033       q->scale(1.0/normz, *z);                        // Set q = z / normz
00034       A->apply(*q, *z);                               // Compute z = A*q
00035       lambda = q->dot(*z);                            // Approximate maximum eigenvalue: lamba = dot(q,z)
00036       if ( iter % 100 == 0 || iter + 1 == niters ) {
00037         r->update(1.0, *z, -lambda, *q, 0.0);         // Compute A*q - lambda*q
00038         residual = Teuchos::ScalarTraits<Scalar>::magnitude(r->norm2() / lambda);
00039         if (verbose) {
00040           std::cout << "Iter = " << iter
00041                     << "  Lambda = " << lambda
00042                     << "  Residual of A*q - lambda*q = " << residual
00043                     << std::endl;
00044         }
00045       }
00046       if (residual < tolerance) {
00047         break;
00048       }
00049     }
00050     return lambda;
00051   }
00052
00053 } // end of namespace TpetraExamples
00054
00063 #endif
```