Epetra Package Browser (Single Doxygen Collection) Development
example/petra_power_method_LL/cxx_main.cpp
Go to the documentation of this file.
00001 //@HEADER
00002 // ************************************************************************
00003 // 
00004 //               Epetra: Linear Algebra Services Package 
00005 //                 Copyright 2011 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 <cstdio>
00044 #include <cstdlib>
00045 #include <cassert>
00046 #include <string>
00047 #include <cmath>
00048 #include <vector>
00049 #include "Epetra_Map.h"
00050 #include "Epetra_Time.h"
00051 #include "Epetra_MultiVector.h"
00052 #include "Epetra_Vector.h"
00053 #include "Epetra_CrsMatrix.h"
00054 #ifdef EPETRA_MPI
00055 #include "mpi.h"
00056 #include "Epetra_MpiComm.h"
00057 #endif
00058 #ifndef __cplusplus
00059 #define __cplusplus
00060 #endif
00061 #include "Epetra_Comm.h"
00062 #include "Epetra_SerialComm.h"
00063 #include "Epetra_Version.h"
00064 
00065 // prototype
00066 int power_method(Epetra_CrsMatrix& A, double & lambda, int niters, double tolerance,
00067      bool verbose);
00068 
00069  
00070 int main(int argc, char *argv[])
00071 {
00072   int ierr = 0, i;
00073 
00074 #ifdef EPETRA_MPI
00075 
00076   // Initialize MPI
00077 
00078   MPI_Init(&argc,&argv);
00079 
00080   Epetra_MpiComm Comm( MPI_COMM_WORLD );
00081 
00082 #else
00083 
00084   Epetra_SerialComm Comm;
00085 
00086 #endif
00087 
00088   int MyPID = Comm.MyPID();
00089   int NumProc = Comm.NumProc();
00090   bool verbose = (MyPID==0);
00091 
00092   if (verbose)
00093     std::cout << Epetra_Version() << std::endl << std::endl;
00094 
00095   std::cout << Comm << std::endl;
00096 
00097   // Get the number of local equations from the command line
00098   if (argc!=2)
00099    {
00100      if (verbose) 
00101        std::cout << "Usage: " << argv[0] << " number_of_equations" << std::endl;
00102     std::exit(1);
00103    }
00104   long long NumGlobalElements = std::atoi(argv[1]);
00105 
00106   if (NumGlobalElements < NumProc)
00107       {
00108      if (verbose)
00109        std::cout << "numGlobalBlocks = " << NumGlobalElements 
00110       << " cannot be < number of processors = " << NumProc << std::endl;
00111      std::exit(1);
00112       }
00113 
00114   // Construct a Map that puts approximately the same number of 
00115   // equations on each processor.
00116 
00117   Epetra_Map Map(NumGlobalElements, 0LL, Comm);
00118   
00119   // Get update list and number of local equations from newly created Map.
00120 
00121   int NumMyElements = Map.NumMyElements();
00122 
00123   std::vector<long long> MyGlobalElements(NumMyElements);
00124     Map.MyGlobalElements(&MyGlobalElements[0]);
00125 
00126   // Create an integer vector NumNz that is used to build the Petra Matrix.
00127   // NumNz[i] is the Number of OFF-DIAGONAL term for the ith global equation 
00128   // on this processor
00129 
00130     std::vector<int> NumNz(NumMyElements);
00131 
00132   // We are building a tridiagonal matrix where each row has (-1 2 -1)
00133   // So we need 2 off-diagonal terms (except for the first and last equation)
00134 
00135   for (i=0; i<NumMyElements; i++)
00136     if (MyGlobalElements[i]==0 || MyGlobalElements[i] == NumGlobalElements-1)
00137       NumNz[i] = 2;
00138     else
00139       NumNz[i] = 3;
00140 
00141   // Create a Epetra_Matrix
00142 
00143   Epetra_CrsMatrix A(Copy, Map, &NumNz[0]);
00144   
00145   // Add  rows one-at-a-time
00146   // Need some vectors to help
00147   // Off diagonal Values will always be -1
00148 
00149 
00150   std::vector<double> Values(2);
00151   Values[0] = -1.0; Values[1] = -1.0;
00152   std::vector<long long> Indices(2);
00153   double two = 2.0;
00154   int NumEntries;
00155   
00156   for (i=0; i<NumMyElements; i++)
00157     {
00158     if (MyGlobalElements[i]==0)
00159       {
00160   Indices[0] = 1;
00161   NumEntries = 1;
00162       }
00163     else if (MyGlobalElements[i] == NumGlobalElements-1)
00164       {
00165   Indices[0] = NumGlobalElements-2;
00166   NumEntries = 1;
00167       }
00168     else
00169       {
00170   Indices[0] = MyGlobalElements[i]-1;
00171   Indices[1] = MyGlobalElements[i]+1;
00172   NumEntries = 2;
00173       }
00174      ierr = A.InsertGlobalValues(MyGlobalElements[i], NumEntries, &Values[0], &Indices[0]);
00175      assert(ierr==0);
00176      // Put in the diagonal entry
00177      ierr = A.InsertGlobalValues(MyGlobalElements[i], 1, &two, &MyGlobalElements[i]);
00178      assert(ierr==0);
00179     }
00180    
00181   // Finish up
00182   ierr = A.FillComplete();
00183   assert(ierr==0);
00184 
00185   // Create vectors for Power method
00186 
00187 
00188   // variable needed for iteration
00189   double lambda = 0.0;
00190   int niters = (int) NumGlobalElements*10;
00191   double tolerance = 1.0e-2;
00192 
00193   // Iterate
00194   Epetra_Flops counter;
00195   A.SetFlopCounter(counter);
00196   Epetra_Time timer(Comm);
00197   ierr += power_method(A, lambda, niters, tolerance, verbose);
00198   double elapsed_time = timer.ElapsedTime();
00199   double total_flops =counter.Flops();
00200   double MFLOPs = total_flops/elapsed_time/1000000.0;
00201 
00202   if (verbose) 
00203     std::cout << "\n\nTotal MFLOPs for first solve = " << MFLOPs << std::endl<< std::endl;
00204 
00205   // Increase diagonal dominance
00206   if (verbose) 
00207     std::cout << "\nIncreasing magnitude of first diagonal term, solving again\n\n"
00208         << std::endl;
00209 
00210   if (A.MyGlobalRow(0)) {
00211     int numvals = A.NumGlobalEntries(0);
00212     std::vector<double> Rowvals(numvals);
00213     std::vector<long long> Rowinds(numvals);
00214     A.ExtractGlobalRowCopy(0, numvals, numvals, &Rowvals[0], &Rowinds[0]); // Get A[0,0]
00215     for (i=0; i<numvals; i++) if (Rowinds[i] == 0) Rowvals[i] *= 10.0;
00216 
00217     A.ReplaceGlobalValues(0, numvals, &Rowvals[0], &Rowinds[0]);
00218   }
00219  
00220   // Iterate (again)
00221   lambda = 0.0;
00222   timer.ResetStartTime();
00223   counter.ResetFlops();
00224   ierr += power_method(A, lambda, niters, tolerance, verbose);
00225   elapsed_time = timer.ElapsedTime();
00226   total_flops = counter.Flops();
00227   MFLOPs = total_flops/elapsed_time/1000000.0;
00228 
00229   if (verbose) 
00230     std::cout << "\n\nTotal MFLOPs for second solve = " << MFLOPs << std::endl<< std::endl;
00231 
00232 
00233   // Release all objects
00234 #ifdef EPETRA_MPI
00235   MPI_Finalize() ;
00236 #endif
00237 
00238 /* end main
00239 */
00240 return ierr ;
00241 }
00242 
00243 int power_method(Epetra_CrsMatrix& A, double &lambda, int niters, 
00244      double tolerance, bool verbose) {  
00245 
00246   Epetra_Vector q(A.RowMap());
00247   Epetra_Vector z(A.RowMap());
00248   Epetra_Vector resid(A.RowMap());
00249 
00250   Epetra_Flops * counter = A.GetFlopCounter();
00251   if (counter!=0) {
00252     q.SetFlopCounter(A);
00253     z.SetFlopCounter(A);
00254     resid.SetFlopCounter(A);
00255   }
00256 
00257   // Fill z with random Numbers
00258   z.Random();
00259 
00260   // variable needed for iteration
00261   double normz, residual;
00262 
00263   int ierr = 1;
00264 
00265   for (int iter = 0; iter < niters; iter++)
00266     {
00267       z.Norm2(&normz); // Compute 2-norm of z
00268       q.Scale(1.0/normz, z);
00269       A.Multiply(false, q, z); // Compute z = A*q
00270       q.Dot(z, &lambda); // Approximate maximum eigenvalue
00271       if (iter%100==0 || iter+1==niters)
00272   {
00273     resid.Update(1.0, z, -lambda, q, 0.0); // Compute A*q - lambda*q
00274     resid.Norm2(&residual);
00275     if (verbose) std::cout << "Iter = " << iter << "  Lambda = " << lambda 
00276           << "  Residual of A*q - lambda*q = " 
00277           << residual << std::endl;
00278   } 
00279       if (residual < tolerance) {
00280   ierr = 0;
00281   break;
00282       }
00283     }
00284   return(ierr);
00285 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines