test/CrsRiluk/cxx_main.cpp

Go to the documentation of this file.
00001 //@HEADER
00002 // ***********************************************************************
00003 // 
00004 //       Ifpack: Object-Oriented Algebraic Preconditioner Package
00005 //                 Copyright (2002) 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 "Ifpack_ConfigDefs.h"
00030 #include <stdio.h>
00031 #include <stdlib.h>
00032 #include <ctype.h>
00033 #include <assert.h>
00034 #include <string.h>
00035 #include <math.h>
00036 #include "Epetra_Comm.h"
00037 #include "Epetra_Map.h"
00038 #include "Epetra_Time.h"
00039 #include "Epetra_CrsMatrix.h"
00040 #include "Epetra_Vector.h"
00041 #include "Epetra_Object.h"
00042 #include "Ifpack_Version.h"
00043 #ifdef EPETRA_MPI
00044 #include "mpi.h"
00045 #include "Epetra_MpiComm.h"
00046 #else
00047 #include "Epetra_SerialComm.h"
00048 #endif
00049 
00050 
00051 // prototype
00052 int power_method(Epetra_CrsMatrix& A, 
00053     Epetra_Vector& q,
00054     Epetra_Vector& z, 
00055     Epetra_Vector& resid, 
00056     double * lambda, int niters, double tolerance,
00057     bool verbose);
00058 
00059 
00060 int main(int argc, char *argv[])
00061 {
00062   int ierr = 0, i, j;
00063 
00064 #ifdef EPETRA_MPI
00065 
00066   // Initialize MPI
00067 
00068   MPI_Init(&argc,&argv);
00069   int size, rank; // Number of MPI processes, My process ID
00070 
00071   MPI_Comm_size(MPI_COMM_WORLD, &size);
00072   MPI_Comm_rank(MPI_COMM_WORLD, &rank);
00073   Epetra_MpiComm Comm( MPI_COMM_WORLD );
00074 
00075 #else
00076 
00077   int size = 1; // Serial case (not using MPI)
00078   int rank = 0;
00079   Epetra_SerialComm Comm;
00080 
00081 #endif
00082 
00083   bool verbose = false;
00084 
00085   // Check if we should print results to standard out
00086   if (argc>1) if (argv[1][0]=='-' && argv[1][1]=='v') verbose = true;
00087 
00088   //  char tmp;
00089   //  if (rank==0) cout << "Press any key to continue..."<< endl;
00090   //  if (rank==0) cin >> tmp;
00091   //  Comm.Barrier();
00092 
00093   int MyPID = Comm.MyPID();
00094   int NumProc = Comm.NumProc();
00095 
00096   if (verbose && MyPID==0)
00097     cout << Ifpack_Version() << endl << endl;
00098 
00099   if (verbose) cout << "Processor "<<MyPID<<" of "<< NumProc
00100               << " is alive."<<endl;
00101 
00102   bool verbose1 = verbose;
00103 
00104   // Redefine verbose to only print on PE 0
00105   if (verbose && rank!=0) verbose = false;
00106 
00107   int NumMyEquations = 10000;
00108   int NumGlobalEquations = NumMyEquations*NumProc+EPETRA_MIN(NumProc,3);
00109   if (MyPID < 3) NumMyEquations++;
00110   int IndexBase = 0;
00111   int ElementSize = 7;
00112   bool DistributedGlobal = (NumGlobalEquations>NumMyEquations);
00113 
00114   // Construct a Map that puts approximately the same Number of equations on each processor
00115 
00116   Epetra_Map Map(NumGlobalEquations, NumMyEquations, 0, Comm);
00117   
00118   // Get update list and number of local equations from newly created Map
00119   int * MyGlobalElements = new int[Map.NumMyElements()];
00120   Map.MyGlobalElements(MyGlobalElements);
00121 
00122   // Create an integer vector NumNz that is used to build the Petra Matrix.
00123   // NumNz[i] is the Number of OFF-DIAGONAL term for the ith global equation on this processor
00124 
00125   int * NumNz = new int[NumMyEquations];
00126 
00127   // We are building a tridiagonal matrix where each row has (-1 2 -1)
00128   // So we need 2 off-diagonal terms (except for the first and last equation)
00129 
00130   for (i=0; i<NumMyEquations; i++)
00131     if (MyGlobalElements[i]==0 || MyGlobalElements[i] == NumGlobalEquations-1)
00132       NumNz[i] = 1;
00133     else
00134       NumNz[i] = 2;
00135 
00136   // Create a Epetra_Matrix
00137 
00138   Epetra_CrsMatrix A(Copy, Map, NumNz);
00139   
00140   // Add  rows one-at-a-time
00141   // Need some vectors to help
00142   // Off diagonal Values will always be -1
00143 
00144 
00145   double *Values = new double[2];
00146   Values[0] = -1.0; Values[1] = -1.0;
00147   int *Indices = new int[2];
00148   double two = 2.0;
00149   int NumEntries;
00150   
00151   for (i=0; i<NumMyEquations; i++)
00152     {
00153     if (MyGlobalElements[i]==0)
00154       {
00155   Indices[0] = 1;
00156   NumEntries = 1;
00157       }
00158     else if (MyGlobalElements[i] == NumGlobalEquations-1)
00159       {
00160   Indices[0] = NumGlobalEquations-2;
00161   NumEntries = 1;
00162       }
00163     else
00164       {
00165   Indices[0] = MyGlobalElements[i]-1;
00166   Indices[1] = MyGlobalElements[i]+1;
00167   NumEntries = 2;
00168       }
00169      int ierr;
00170      ierr = A.InsertGlobalValues(MyGlobalElements[i], NumEntries, Values, Indices);
00171      IFPACK_CHK_ERR(ierr);
00172      ierr = A.InsertGlobalValues(MyGlobalElements[i], 1, &two, MyGlobalElements+i); // Put in the diagonal entry
00173      IFPACK_CHK_ERR(ierr);
00174     }
00175   
00176   // Finish up
00177   A.FillComplete();
00178 
00179   // Create vectors for Power method
00180 
00181   Epetra_Vector q(Map);
00182   Epetra_Vector z(Map);
00183   Epetra_Vector resid(Map);
00184 
00185   // variable needed for iteration
00186   double lambda = 0.0;
00187   int niters = 10000;
00188   double tolerance = 1.0e-3;
00189 
00190   // Iterate
00191   Epetra_Time timer(Comm);
00192   ierr += power_method(A, q, z, resid, &lambda, niters, tolerance, verbose);
00193   double elapsed_time = timer.ElapsedTime();
00194   double total_flops = A.Flops() + q.Flops() + z.Flops() + resid.Flops();
00195   double MFLOPs = total_flops/elapsed_time/1000000.0;
00196 
00197   if (verbose) cout << "\n\nTotal MFLOPs for first solve = " << MFLOPs << endl<< endl;
00198 
00199   // Increase diagonal dominance
00200   if (verbose) cout << "\n\nIncreasing the magnitude of first diagonal term and solving again\n\n"
00201         << endl;
00202 
00203   
00204   if (A.MyGlobalRow(0)) {
00205     int numvals = A.NumGlobalEntries(0);
00206     double * Rowvals = new double [numvals];
00207     int    * Rowinds = new int    [numvals];
00208     A.ExtractGlobalRowCopy(0, numvals, numvals, Rowvals, Rowinds); // Get A[0,0]
00209 
00210     for (i=0; i<numvals; i++) if (Rowinds[i] == 0) Rowvals[i] *= 10.0;
00211     
00212     A.ReplaceGlobalValues(0, numvals, Rowvals, Rowinds);
00213   }
00214   // Iterate (again)
00215   lambda = 0.0;
00216   timer.ResetStartTime();
00217   A.ResetFlops(); q.ResetFlops(); z.ResetFlops(); resid.ResetFlops();
00218   ierr += power_method(A, q, z, resid, &lambda, niters, tolerance, verbose);
00219   elapsed_time = timer.ElapsedTime();
00220   total_flops = A.Flops() + q.Flops() + z.Flops() + resid.Flops();
00221   MFLOPs = total_flops/elapsed_time/1000000.0;
00222 
00223   if (verbose) cout << "\n\nTotal MFLOPs for second solve = " << MFLOPs << endl<< endl;
00224 
00225 
00226   // Release all objects
00227   delete [] NumNz;
00228   delete [] Values;
00229   delete [] Indices;
00230   delete [] MyGlobalElements;
00231 
00232       
00233 
00234   if (verbose1) {
00235     // Test ostream << operator (if verbose1)
00236     // Construct a Map that puts 2 equations on each PE
00237     
00238     int NumMyElements1 = 2;
00239     int NumMyEquations1 = NumMyElements1;
00240     int NumGlobalEquations1 = NumMyEquations1*NumProc;
00241 
00242     Epetra_Map Map1(-1, NumMyElements1, 0, Comm);
00243     
00244     // Get update list and number of local equations from newly created Map
00245     int * MyGlobalElements1 = new int[Map1.NumMyElements()];
00246     Map1.MyGlobalElements(MyGlobalElements1);
00247     
00248     // Create an integer vector NumNz that is used to build the Petra Matrix.
00249     // NumNz[i] is the Number of OFF-DIAGONAL term for the ith global equation on this processor
00250     
00251     int * NumNz1 = new int[NumMyEquations1];
00252     
00253     // We are building a tridiagonal matrix where each row has (-1 2 -1)
00254     // So we need 2 off-diagonal terms (except for the first and last equation)
00255     
00256     for (i=0; i<NumMyEquations1; i++)
00257       if (MyGlobalElements1[i]==0 || MyGlobalElements1[i] == NumGlobalEquations1-1)
00258   NumNz1[i] = 1;
00259       else
00260   NumNz1[i] = 2;
00261     
00262     // Create a Epetra_Matrix
00263     
00264     Epetra_CrsMatrix A1(Copy, Map1, NumNz1);
00265     
00266     // Add  rows one-at-a-time
00267     // Need some vectors to help
00268     // Off diagonal Values will always be -1
00269     
00270     
00271     double *Values1 = new double[2];
00272     Values1[0] = -1.0; Values1[1] = -1.0;
00273     int *Indices1 = new int[2];
00274     double two1 = 2.0;
00275     int NumEntries1;
00276     
00277     for (i=0; i<NumMyEquations1; i++)
00278       {
00279   if (MyGlobalElements1[i]==0)
00280     {
00281       Indices1[0] = 1;
00282       NumEntries1 = 1;
00283     }
00284   else if (MyGlobalElements1[i] == NumGlobalEquations1-1)
00285     {
00286       Indices1[0] = NumGlobalEquations1-2;
00287       NumEntries1 = 1;
00288     }
00289   else
00290     {
00291       Indices1[0] = MyGlobalElements1[i]-1;
00292       Indices1[1] = MyGlobalElements1[i]+1;
00293       NumEntries1 = 2;
00294     }
00295         int ierr;
00296   ierr = A1.InsertGlobalValues(MyGlobalElements1[i], NumEntries1, Values1, Indices1);
00297         IFPACK_CHK_ERR(ierr);
00298   ierr = A1.InsertGlobalValues(MyGlobalElements1[i], 1, &two1, MyGlobalElements1+i); // Put in the diagonal entry
00299         IFPACK_CHK_ERR(ierr);
00300       }
00301     
00302     // Finish up
00303     A1.FillComplete();
00304     
00305     if (verbose) cout << "\n\nPrint out tridiagonal matrix, each part on each processor.\n\n" << endl;
00306     cout << A1 << endl;
00307     
00308   // Release all objects
00309   delete [] NumNz1;
00310   delete [] Values1;
00311   delete [] Indices1;
00312   delete [] MyGlobalElements1;
00313 
00314   }
00315       
00316 #ifdef EPETRA_MPI
00317   MPI_Finalize() ;
00318 #endif
00319 
00320 /* end main
00321 */
00322 return ierr ;
00323 }
00324 
00325 int power_method(Epetra_CrsMatrix& A, 
00326      Epetra_Vector& q,
00327      Epetra_Vector& z, 
00328      Epetra_Vector& resid, 
00329      double * lambda, int niters, double tolerance,
00330      bool verbose) {  
00331 
00332   // Fill z with random Numbers
00333   z.Random();
00334 
00335   // variable needed for iteration
00336   double normz, residual;
00337 
00338   int ierr = 1;
00339 
00340   for (int iter = 0; iter < niters; iter++)
00341     {
00342       z.Norm2(&normz); // Compute 2-norm of z
00343       q.Scale(1.0/normz, z);
00344       A.Multiply(false, q, z); // Compute z = A*q
00345       q.Dot(z, lambda); // Approximate maximum eigenvaluE
00346       if (iter%100==0 || iter+1==niters)
00347   {
00348     resid.Update(1.0, z, -(*lambda), q, 0.0); // Compute A*q - lambda*q
00349     resid.Norm2(&residual);
00350     if (verbose) cout << "Iter = " << iter << "  Lambda = " << *lambda 
00351            << "  Residual of A*q - lambda*q = " << residual << endl;
00352   } 
00353       if (residual < tolerance) {
00354   ierr = 0;
00355   break;
00356       }
00357     }
00358   return(ierr);
00359 }

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