Epetra Package Browser (Single Doxygen Collection) Development
power_method.cpp
Go to the documentation of this file.
00001 //@HEADER
00002 // ************************************************************************
00003 // 
00004 //               Epetra: Linear Algebra Services Package 
00005 //                 Copyright 2001 Sandia Corporation
00006 // 
00007 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00008 // the U.S. Government retains certain rights in this software.
00009 //
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00038 // 
00039 // ************************************************************************
00040 //@HEADER
00041 
00042 
00043 #include <stdio.h>
00044 #include <stdlib.h>
00045 #include "Epetra_Comm.h"
00046 #include "Epetra_Map.h"
00047 #include "Epetra_Vector.h"
00048 #include "Epetra_CrsMatrix.h"
00049 // Simple Power method algorithm
00050 double power_method(const Epetra_CrsMatrix& A) {  
00051   // variable needed for iteration
00052   double lambda = 0.0;
00053   int niters = A.RowMap().NumGlobalElements()*10;
00054   double tolerance = 1.0e-10;
00055   // Create vectors
00056   Epetra_Vector q(A.RowMap());
00057   Epetra_Vector z(A.RowMap());
00058   Epetra_Vector resid(A.RowMap());
00059   // Fill z with random Numbers
00060   z.Random();
00061   // variable needed for iteration
00062   double normz;
00063   double residual = 0;
00064   int iter = 0;
00065   while (iter==0 || (iter < niters && residual > tolerance)) {
00066     z.Norm2(&normz); // Compute 2-norm of z
00067     q.Scale(1.0/normz, z);
00068     A.Multiply(false, q, z); // Compute z = A*q
00069     q.Dot(z, &lambda); // Approximate maximum eigenvalue
00070     if (iter%10==0 || iter+1==niters) {
00071       // Compute A*q - lambda*q every 10 iterations
00072       resid.Update(1.0, z, -lambda, q, 0.0);
00073       resid.Norm2(&residual);
00074       if (q.Map().Comm().MyPID()==0)
00075   cout << "Iter = " << iter << "  Lambda = " << lambda 
00076        << "  Two-norm of A*q - lambda*q = " 
00077        << residual << endl;
00078     } 
00079     iter++;
00080   }
00081   return(lambda);
00082 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines