example/petra_power_method/cxx_main.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 terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //  
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //  
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00025 // 
00026 // ************************************************************************
00027 //@HEADER
00028 
00029 #include <cstdio>
00030 #include <cstdlib>
00031 #include <cassert>
00032 #include <string>
00033 #include <cmath>
00034 #include <vector>
00035 #include "Epetra_Map.h"
00036 #include "Epetra_Time.h"
00037 #include "Epetra_MultiVector.h"
00038 #include "Epetra_Vector.h"
00039 #include "Epetra_CrsMatrix.h"
00040 #ifdef EPETRA_MPI
00041 #include "mpi.h"
00042 #include "Epetra_MpiComm.h"
00043 #endif
00044 #ifndef __cplusplus
00045 #define __cplusplus
00046 #endif
00047 #include "Epetra_Comm.h"
00048 #include "Epetra_SerialComm.h"
00049 #include "Epetra_Version.h"
00050 
00051 // prototype
00052 int power_method(Epetra_CrsMatrix& A, double & lambda, int niters, double tolerance,
00053      bool verbose);
00054 
00055  
00056 int main(int argc, char *argv[])
00057 {
00058   int ierr = 0, i;
00059 
00060 #ifdef EPETRA_MPI
00061 
00062   // Initialize MPI
00063 
00064   MPI_Init(&argc,&argv);
00065 
00066   Epetra_MpiComm Comm( MPI_COMM_WORLD );
00067 
00068 #else
00069 
00070   Epetra_SerialComm Comm;
00071 
00072 #endif
00073 
00074   int MyPID = Comm.MyPID();
00075   int NumProc = Comm.NumProc();
00076   bool verbose = (MyPID==0);
00077 
00078   if (verbose)
00079     cout << Epetra_Version() << endl << endl;
00080 
00081   cout << Comm << endl;
00082 
00083   // Get the number of local equations from the command line
00084   if (argc!=2)
00085    {
00086      if (verbose) 
00087        cout << "Usage: " << argv[0] << " number_of_equations" << endl;
00088     std::exit(1);
00089    }
00090   int NumGlobalElements = std::atoi(argv[1]);
00091 
00092   if (NumGlobalElements < NumProc)
00093       {
00094      if (verbose)
00095        cout << "numGlobalBlocks = " << NumGlobalElements 
00096       << " cannot be < number of processors = " << NumProc << endl;
00097      std::exit(1);
00098       }
00099 
00100   // Construct a Map that puts approximately the same number of 
00101   // equations on each processor.
00102 
00103   Epetra_Map Map(NumGlobalElements, 0, Comm);
00104   
00105   // Get update list and number of local equations from newly created Map.
00106 
00107   int NumMyElements = Map.NumMyElements();
00108 
00109   std::vector<int> MyGlobalElements(NumMyElements);
00110     Map.MyGlobalElements(&MyGlobalElements[0]);
00111 
00112   // Create an integer vector NumNz that is used to build the Petra Matrix.
00113   // NumNz[i] is the Number of OFF-DIAGONAL term for the ith global equation 
00114   // on this processor
00115 
00116     std::vector<int> NumNz(NumMyElements);
00117 
00118   // We are building a tridiagonal matrix where each row has (-1 2 -1)
00119   // So we need 2 off-diagonal terms (except for the first and last equation)
00120 
00121   for (i=0; i<NumMyElements; i++)
00122     if (MyGlobalElements[i]==0 || MyGlobalElements[i] == NumGlobalElements-1)
00123       NumNz[i] = 2;
00124     else
00125       NumNz[i] = 3;
00126 
00127   // Create a Epetra_Matrix
00128 
00129   Epetra_CrsMatrix A(Copy, Map, &NumNz[0]);
00130   
00131   // Add  rows one-at-a-time
00132   // Need some vectors to help
00133   // Off diagonal Values will always be -1
00134 
00135 
00136   std::vector<double> Values(2);
00137   Values[0] = -1.0; Values[1] = -1.0;
00138   std::vector<int> Indices(2);
00139   double two = 2.0;
00140   int NumEntries;
00141   
00142   for (i=0; i<NumMyElements; i++)
00143     {
00144     if (MyGlobalElements[i]==0)
00145       {
00146   Indices[0] = 1;
00147   NumEntries = 1;
00148       }
00149     else if (MyGlobalElements[i] == NumGlobalElements-1)
00150       {
00151   Indices[0] = NumGlobalElements-2;
00152   NumEntries = 1;
00153       }
00154     else
00155       {
00156   Indices[0] = MyGlobalElements[i]-1;
00157   Indices[1] = MyGlobalElements[i]+1;
00158   NumEntries = 2;
00159       }
00160      ierr = A.InsertGlobalValues(MyGlobalElements[i], NumEntries, &Values[0], &Indices[0]);
00161      assert(ierr==0);
00162      // Put in the diagonal entry
00163      ierr = A.InsertGlobalValues(MyGlobalElements[i], 1, &two, &MyGlobalElements[i]);
00164      assert(ierr==0);
00165     }
00166    
00167   // Finish up
00168   ierr = A.FillComplete();
00169   assert(ierr==0);
00170 
00171   // Create vectors for Power method
00172 
00173 
00174   // variable needed for iteration
00175   double lambda = 0.0;
00176   int niters = NumGlobalElements*10;
00177   double tolerance = 1.0e-2;
00178 
00179   // Iterate
00180   Epetra_Flops counter;
00181   A.SetFlopCounter(counter);
00182   Epetra_Time timer(Comm);
00183   ierr += power_method(A, lambda, niters, tolerance, verbose);
00184   double elapsed_time = timer.ElapsedTime();
00185   double total_flops =counter.Flops();
00186   double MFLOPs = total_flops/elapsed_time/1000000.0;
00187 
00188   if (verbose) 
00189     cout << "\n\nTotal MFLOPs for first solve = " << MFLOPs << endl<< endl;
00190 
00191   // Increase diagonal dominance
00192   if (verbose) 
00193     cout << "\nIncreasing magnitude of first diagonal term, solving again\n\n"
00194         << endl;
00195 
00196   if (A.MyGlobalRow(0)) {
00197     int numvals = A.NumGlobalEntries(0);
00198     std::vector<double> Rowvals(numvals);
00199     std::vector<int> Rowinds(numvals);
00200     A.ExtractGlobalRowCopy(0, numvals, numvals, &Rowvals[0], &Rowinds[0]); // Get A[0,0]
00201     for (i=0; i<numvals; i++) if (Rowinds[i] == 0) Rowvals[i] *= 10.0;
00202 
00203     A.ReplaceGlobalValues(0, numvals, &Rowvals[0], &Rowinds[0]);
00204   }
00205  
00206   // Iterate (again)
00207   lambda = 0.0;
00208   timer.ResetStartTime();
00209   counter.ResetFlops();
00210   ierr += power_method(A, lambda, niters, tolerance, verbose);
00211   elapsed_time = timer.ElapsedTime();
00212   total_flops = counter.Flops();
00213   MFLOPs = total_flops/elapsed_time/1000000.0;
00214 
00215   if (verbose) 
00216     cout << "\n\nTotal MFLOPs for second solve = " << MFLOPs << endl<< endl;
00217 
00218 
00219   // Release all objects
00220 #ifdef EPETRA_MPI
00221   MPI_Finalize() ;
00222 #endif
00223 
00224 /* end main
00225 */
00226 return ierr ;
00227 }
00228 
00229 int power_method(Epetra_CrsMatrix& A, double &lambda, int niters, 
00230      double tolerance, bool verbose) {  
00231 
00232   Epetra_Vector q(A.RowMap());
00233   Epetra_Vector z(A.RowMap());
00234   Epetra_Vector resid(A.RowMap());
00235 
00236   Epetra_Flops * counter = A.GetFlopCounter();
00237   if (counter!=0) {
00238     q.SetFlopCounter(A);
00239     z.SetFlopCounter(A);
00240     resid.SetFlopCounter(A);
00241   }
00242 
00243   // Fill z with random Numbers
00244   z.Random();
00245 
00246   // variable needed for iteration
00247   double normz, residual;
00248 
00249   int ierr = 1;
00250 
00251   for (int iter = 0; iter < niters; iter++)
00252     {
00253       z.Norm2(&normz); // Compute 2-norm of z
00254       q.Scale(1.0/normz, z);
00255       A.Multiply(false, q, z); // Compute z = A*q
00256       q.Dot(z, &lambda); // Approximate maximum eigenvalue
00257       if (iter%100==0 || iter+1==niters)
00258   {
00259     resid.Update(1.0, z, -lambda, q, 0.0); // Compute A*q - lambda*q
00260     resid.Norm2(&residual);
00261     if (verbose) cout << "Iter = " << iter << "  Lambda = " << lambda 
00262           << "  Residual of A*q - lambda*q = " 
00263           << residual << endl;
00264   } 
00265       if (residual < tolerance) {
00266   ierr = 0;
00267   break;
00268       }
00269     }
00270   return(ierr);
00271 }

Generated on Thu Sep 18 12:37:56 2008 for Epetra Package Browser (Single Doxygen Collection) by doxygen 1.3.9.1