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;
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 
00111   // Construct a Map that puts approximately the same Number of equations on each processor
00112 
00113   Epetra_Map Map(NumGlobalEquations, NumMyEquations, 0, Comm);
00114   
00115   // Get update list and number of local equations from newly created Map
00116   int * MyGlobalElements = new int[Map.NumMyElements()];
00117   Map.MyGlobalElements(MyGlobalElements);
00118 
00119   // Create an integer vector NumNz that is used to build the Petra Matrix.
00120   // NumNz[i] is the Number of OFF-DIAGONAL term for the ith global equation on this processor
00121 
00122   int * NumNz = new int[NumMyEquations];
00123 
00124   // We are building a tridiagonal matrix where each row has (-1 2 -1)
00125   // So we need 2 off-diagonal terms (except for the first and last equation)
00126 
00127   for (i=0; i<NumMyEquations; i++)
00128     if (MyGlobalElements[i]==0 || MyGlobalElements[i] == NumGlobalEquations-1)
00129       NumNz[i] = 1;
00130     else
00131       NumNz[i] = 2;
00132 
00133   // Create a Epetra_Matrix
00134 
00135   Epetra_CrsMatrix A(Copy, Map, NumNz);
00136   
00137   // Add  rows one-at-a-time
00138   // Need some vectors to help
00139   // Off diagonal Values will always be -1
00140 
00141 
00142   double *Values = new double[2];
00143   Values[0] = -1.0; Values[1] = -1.0;
00144   int *Indices = new int[2];
00145   double two = 2.0;
00146   int NumEntries;
00147   
00148   for (i=0; i<NumMyEquations; i++)
00149     {
00150     if (MyGlobalElements[i]==0)
00151       {
00152   Indices[0] = 1;
00153   NumEntries = 1;
00154       }
00155     else if (MyGlobalElements[i] == NumGlobalEquations-1)
00156       {
00157   Indices[0] = NumGlobalEquations-2;
00158   NumEntries = 1;
00159       }
00160     else
00161       {
00162   Indices[0] = MyGlobalElements[i]-1;
00163   Indices[1] = MyGlobalElements[i]+1;
00164   NumEntries = 2;
00165       }
00166      int ierr;
00167      ierr = A.InsertGlobalValues(MyGlobalElements[i], NumEntries, Values, Indices);
00168      IFPACK_CHK_ERR(ierr);
00169      ierr = A.InsertGlobalValues(MyGlobalElements[i], 1, &two, MyGlobalElements+i); // Put in the diagonal entry
00170      IFPACK_CHK_ERR(ierr);
00171     }
00172   
00173   // Finish up
00174   A.FillComplete();
00175 
00176   // Create vectors for Power method
00177 
00178   Epetra_Vector q(Map);
00179   Epetra_Vector z(Map);
00180   Epetra_Vector resid(Map);
00181 
00182   // variable needed for iteration
00183   double lambda = 0.0;
00184   int niters = 10000;
00185   double tolerance = 1.0e-3;
00186 
00187   // Iterate
00188   Epetra_Time timer(Comm);
00189   ierr += power_method(A, q, z, resid, &lambda, niters, tolerance, verbose);
00190   double elapsed_time = timer.ElapsedTime();
00191   double total_flops = A.Flops() + q.Flops() + z.Flops() + resid.Flops();
00192   double MFLOPs = total_flops/elapsed_time/1000000.0;
00193 
00194   if (verbose) cout << "\n\nTotal MFLOPs for first solve = " << MFLOPs << endl<< endl;
00195 
00196   // Increase diagonal dominance
00197   if (verbose) cout << "\n\nIncreasing the magnitude of first diagonal term and solving again\n\n"
00198         << endl;
00199 
00200   
00201   if (A.MyGlobalRow(0)) {
00202     int numvals = A.NumGlobalEntries(0);
00203     double * Rowvals = new double [numvals];
00204     int    * Rowinds = new int    [numvals];
00205     A.ExtractGlobalRowCopy(0, numvals, numvals, Rowvals, Rowinds); // Get A[0,0]
00206 
00207     for (i=0; i<numvals; i++) if (Rowinds[i] == 0) Rowvals[i] *= 10.0;
00208     
00209     A.ReplaceGlobalValues(0, numvals, Rowvals, Rowinds);
00210   }
00211   // Iterate (again)
00212   lambda = 0.0;
00213   timer.ResetStartTime();
00214   A.ResetFlops(); q.ResetFlops(); z.ResetFlops(); resid.ResetFlops();
00215   ierr += power_method(A, q, z, resid, &lambda, niters, tolerance, verbose);
00216   elapsed_time = timer.ElapsedTime();
00217   total_flops = A.Flops() + q.Flops() + z.Flops() + resid.Flops();
00218   MFLOPs = total_flops/elapsed_time/1000000.0;
00219 
00220   if (verbose) cout << "\n\nTotal MFLOPs for second solve = " << MFLOPs << endl<< endl;
00221 
00222 
00223   // Release all objects
00224   delete [] NumNz;
00225   delete [] Values;
00226   delete [] Indices;
00227   delete [] MyGlobalElements;
00228 
00229       
00230 
00231   if (verbose1) {
00232     // Test ostream << operator (if verbose1)
00233     // Construct a Map that puts 2 equations on each PE
00234     
00235     int NumMyElements1 = 2;
00236     int NumMyEquations1 = NumMyElements1;
00237     int NumGlobalEquations1 = NumMyEquations1*NumProc;
00238 
00239     Epetra_Map Map1(-1, NumMyElements1, 0, Comm);
00240     
00241     // Get update list and number of local equations from newly created Map
00242     int * MyGlobalElements1 = new int[Map1.NumMyElements()];
00243     Map1.MyGlobalElements(MyGlobalElements1);
00244     
00245     // Create an integer vector NumNz that is used to build the Petra Matrix.
00246     // NumNz[i] is the Number of OFF-DIAGONAL term for the ith global equation on this processor
00247     
00248     int * NumNz1 = new int[NumMyEquations1];
00249     
00250     // We are building a tridiagonal matrix where each row has (-1 2 -1)
00251     // So we need 2 off-diagonal terms (except for the first and last equation)
00252     
00253     for (i=0; i<NumMyEquations1; i++)
00254       if (MyGlobalElements1[i]==0 || MyGlobalElements1[i] == NumGlobalEquations1-1)
00255   NumNz1[i] = 1;
00256       else
00257   NumNz1[i] = 2;
00258     
00259     // Create a Epetra_Matrix
00260     
00261     Epetra_CrsMatrix A1(Copy, Map1, NumNz1);
00262     
00263     // Add  rows one-at-a-time
00264     // Need some vectors to help
00265     // Off diagonal Values will always be -1
00266     
00267     
00268     double *Values1 = new double[2];
00269     Values1[0] = -1.0; Values1[1] = -1.0;
00270     int *Indices1 = new int[2];
00271     double two1 = 2.0;
00272     int NumEntries1;
00273     
00274     for (i=0; i<NumMyEquations1; i++)
00275       {
00276   if (MyGlobalElements1[i]==0)
00277     {
00278       Indices1[0] = 1;
00279       NumEntries1 = 1;
00280     }
00281   else if (MyGlobalElements1[i] == NumGlobalEquations1-1)
00282     {
00283       Indices1[0] = NumGlobalEquations1-2;
00284       NumEntries1 = 1;
00285     }
00286   else
00287     {
00288       Indices1[0] = MyGlobalElements1[i]-1;
00289       Indices1[1] = MyGlobalElements1[i]+1;
00290       NumEntries1 = 2;
00291     }
00292         int ierr;
00293   ierr = A1.InsertGlobalValues(MyGlobalElements1[i], NumEntries1, Values1, Indices1);
00294         IFPACK_CHK_ERR(ierr);
00295   ierr = A1.InsertGlobalValues(MyGlobalElements1[i], 1, &two1, MyGlobalElements1+i); // Put in the diagonal entry
00296         IFPACK_CHK_ERR(ierr);
00297       }
00298     
00299     // Finish up
00300     A1.FillComplete();
00301     
00302     if (verbose) cout << "\n\nPrint out tridiagonal matrix, each part on each processor.\n\n" << endl;
00303     cout << A1 << endl;
00304     
00305   // Release all objects
00306   delete [] NumNz1;
00307   delete [] Values1;
00308   delete [] Indices1;
00309   delete [] MyGlobalElements1;
00310 
00311   }
00312       
00313 #ifdef EPETRA_MPI
00314   MPI_Finalize() ;
00315 #endif
00316 
00317 /* end main
00318 */
00319 return ierr ;
00320 }
00321 
00322 int power_method(Epetra_CrsMatrix& A, 
00323      Epetra_Vector& q,
00324      Epetra_Vector& z, 
00325      Epetra_Vector& resid, 
00326      double * lambda, int niters, double tolerance,
00327      bool verbose) {  
00328 
00329   // Fill z with random Numbers
00330   z.Random();
00331 
00332   // variable needed for iteration
00333   double normz, residual;
00334 
00335   int ierr = 1;
00336 
00337   for (int iter = 0; iter < niters; iter++)
00338     {
00339       z.Norm2(&normz); // Compute 2-norm of z
00340       q.Scale(1.0/normz, z);
00341       A.Multiply(false, q, z); // Compute z = A*q
00342       q.Dot(z, lambda); // Approximate maximum eigenvaluE
00343       if (iter%100==0 || iter+1==niters)
00344   {
00345     resid.Update(1.0, z, -(*lambda), q, 0.0); // Compute A*q - lambda*q
00346     resid.Norm2(&residual);
00347     if (verbose) cout << "Iter = " << iter << "  Lambda = " << *lambda 
00348            << "  Residual of A*q - lambda*q = " << residual << endl;
00349   } 
00350       if (residual < tolerance) {
00351   ierr = 0;
00352   break;
00353       }
00354     }
00355   return(ierr);
00356 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
Generated on Wed Apr 13 10:05:33 2011 for Ifpack Package Browser (Single Doxygen Collection) by  doxygen 1.6.3