test/IlukGraph/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 // FIXME: assert below are not correct if NDEBUG is defined
00030 #ifdef NDEBUG
00031 #undef NDEBUG
00032 #endif
00033 
00034 #include "Ifpack_ConfigDefs.h"
00035 #include <stdio.h>
00036 #include <stdlib.h>
00037 #include <ctype.h>
00038 #include <assert.h>
00039 #include <string.h>
00040 #include <math.h>
00041 #include "Epetra_Map.h"
00042 #include "Ifpack_Version.h"
00043 #include "Epetra_CrsGraph.h"
00044 #include "Ifpack_IlukGraph.h"
00045 #ifdef EPETRA_MPI
00046 #include "Epetra_MpiComm.h"
00047 #else
00048 #include "Epetra_SerialComm.h"
00049 #endif
00050 
00051 // Prototype
00052 
00053 int check(Epetra_CrsGraph& L, Epetra_CrsGraph& U, Ifpack_IlukGraph& LU, 
00054     int NumGlobalRows1, int NumMyRows1, int LevelFill1, bool verbose);
00055 
00056  int main(int argc, char *argv[])
00057 {
00058   int ierr = 0, i, j;
00059   int * Indices;
00060   int nx, ny;
00061 
00062 #ifdef EPETRA_MPI
00063 
00064   // Initialize MPI
00065 
00066   MPI_Init(&argc,&argv);
00067   //int size, rank; // Number of MPI processes, My process ID
00068 
00069   //MPI_Comm_size(MPI_COMM_WORLD, &size);
00070   //MPI_Comm_rank(MPI_COMM_WORLD, &rank);
00071   Epetra_MpiComm Comm( MPI_COMM_WORLD );
00072 
00073 #else
00074 
00075   //int size = 1; // Serial case (not using MPI)
00076   //int rank = 0;
00077 
00078   Epetra_SerialComm Comm;
00079 #endif
00080 
00081   bool verbose = false;
00082 
00083   int nextarg = 1;
00084   // Check if we should print results to standard out
00085   if (argc>1) if (argv[1][0]=='-' && argv[1][1]=='v') {
00086     verbose = true;
00087     nextarg++;
00088   }
00089 
00090   // char tmp;
00091   // if (rank==0) cout << "Press any key to continue..."<< endl;
00092   // if (rank==0) cin >> tmp;
00093   // Comm.Barrier();
00094 
00095   int MyPID = Comm.MyPID();
00096   int NumProc = Comm.NumProc();
00097 
00098   if (verbose && MyPID==0)
00099     cout << Ifpack_Version() << endl << endl;
00100 
00101   if (verbose) cout << Comm <<endl;
00102 
00103   //int sqrtNumProc = (int) ceil(sqrt((double) NumProc));
00104 
00105   bool verbose1 = verbose;
00106   verbose = verbose && (MyPID==0);
00107 
00108   if (verbose1 && argc != 4) {
00109     nx = 10;
00110     ny = 12*NumProc;
00111     cout << "Setting nx = " << nx << ", ny = " << ny << endl;
00112   }
00113     else if (!verbose1 && argc != 3) {
00114     nx = 10;
00115     ny = 12*NumProc;
00116     }
00117   else {
00118     nx = atoi(argv[nextarg++]);
00119     if (nx<3) {cout << "nx = " << nx << ": Must be greater than 2 for meaningful graph." << endl; exit(1);}
00120     ny = atoi(argv[nextarg++]);
00121     if (ny<3) {cout << "ny = " << ny << ": Must be greater than 2 for meaningful graph." << endl; exit(1);}
00122   }
00123   
00124   int NumGlobalPoints = nx*ny;
00125   int IndexBase = 0;
00126 
00127   if (verbose)
00128     cout << "\n\n*****Building 5 point matrix, Level 1 and 2 filled matrices for" << endl
00129    << "  nx = " << nx << ",  ny = " << ny << endl<< endl;
00130 
00131 
00132   // Create a 5 point stencil graph, level 1 fill of it and level 2 fill of it
00133 
00134   Epetra_Map Map(NumGlobalPoints, IndexBase, Comm);
00135 
00136   int NumMyPoints = Map.NumMyPoints();
00137 
00138   Epetra_CrsGraph A(Copy, Map, 5);
00139   Epetra_CrsGraph L0(Copy, Map, Map, 2);
00140   Epetra_CrsGraph U0(Copy, Map, Map, 2);
00141   Epetra_CrsGraph L1(Copy, Map, Map, 3);
00142   Epetra_CrsGraph U1(Copy, Map, Map, 3);
00143   Epetra_CrsGraph L2(Copy, Map, Map, 4);
00144   Epetra_CrsGraph U2(Copy, Map, Map, 4);
00145   
00146   // Add  rows one-at-a-time
00147 
00148   Indices = new int[4]; // Work space
00149   
00150   for (j=0; j<ny; j++) {
00151     for (i=0; i<nx; i++) {
00152       int Row = i+j*nx;
00153       if (Map.MyGID(Row)) { // Only work on rows I own
00154   
00155   //**** Work on lower triangle of all three matrices ****
00156   
00157   // Define entries (i-1,j), (i,j-1)
00158   
00159   int k = 0;
00160   if (i>0)    Indices[k++] = i-1 + j   *nx;
00161   if (j>0)    Indices[k++] = i   +(j-1)*nx;
00162   
00163   // Define lower triangular terms of original matrix and L(0)
00164   assert(A.InsertGlobalIndices(Row, k, Indices)>=0);
00165   assert(L0.InsertGlobalIndices(Row, k, Indices)>=0);
00166   
00167   // Define entry (i+1,j-1)
00168   if ((i<nx-1) && (j>0   )) Indices[k++] = i+1 +(j-1)*nx;
00169   
00170   
00171   // Define lower triangle of level(1) fill matrix
00172   assert(L1.InsertGlobalIndices(Row, k, Indices)>=0);
00173   
00174   // Define entry (i+2, j-1)
00175   
00176   if ((i<nx-2) && (j>0   )) Indices[k++] = i+2 +(j-1)*nx;
00177   
00178   // Define lower triangle of level(2) fill matrix
00179   assert(L2.InsertGlobalIndices(Row, k, Indices)>=0);
00180   
00181   // Define main diagonal of original matrix
00182   assert(A.InsertGlobalIndices(Row, 1, &Row)>=0);
00183   
00184   k = 0; // Reset index counter
00185   
00186   //**** Work on upper triangle of all three matrices ****
00187   
00188   // Define entries (i+1,j), ( i,j+1)
00189   
00190   if (i<nx-1) Indices[k++] = i+1 + j   *nx;
00191   if (j<ny-1) Indices[k++] = i   +(j+1)*nx;
00192   
00193   // Define upper  triangular terms of original matrix and L(0)
00194   assert(A.InsertGlobalIndices(Row, k, Indices)>=0);
00195   assert(U0.InsertGlobalIndices(Row, k, Indices)>=0);
00196   
00197   // Define entry (i-1,j+1)
00198   
00199   if ((i>0   ) && (j<ny-1)) Indices[k++] = i-1 +(j+1)*nx;
00200   
00201   // Define upper triangle of level(1) fill matrix
00202   assert(U1.InsertGlobalIndices(Row, k, Indices)>=0);
00203   
00204   // Define entry (i-2, j+1)
00205   
00206   if ((i>1   ) && (j<ny-1)) Indices[k++] = i-2 +(j+1)*nx;
00207   
00208   // Define upper triangle of level(2) fill matrix
00209   assert(U2.InsertGlobalIndices(Row, k, Indices)>=0);
00210       }
00211     }
00212   }
00213 
00214   delete [] Indices;
00215 
00216   // Finish up
00217   if (verbose) cout << "\n\nCompleting A" << endl<< endl;
00218   assert(A.FillComplete()==0);
00219   if (verbose) cout << "\n\nCompleting L0" << endl<< endl;
00220   assert(L0.FillComplete()==0);
00221   if (verbose) cout << "\n\nCompleting U0" << endl<< endl;
00222   assert(U0.FillComplete()==0);
00223   if (verbose) cout << "\n\nCompleting L1" << endl<< endl;
00224   assert(L1.FillComplete()==0);
00225   if (verbose) cout << "\n\nCompleting U1" << endl<< endl;
00226   assert(U1.FillComplete()==0);
00227   if (verbose) cout << "\n\nCompleting L2" << endl<< endl;
00228   assert(L2.FillComplete()==0);
00229   if (verbose) cout << "\n\nCompleting U2" << endl<< endl;
00230   assert(U2.FillComplete()==0);
00231 
00232   if (verbose) cout << "\n\n*****Testing ILU(0) constructor on A" << endl<< endl;
00233 
00234   Ifpack_IlukGraph ILU0(A, 0, 0);
00235   assert(ILU0.ConstructFilledGraph()==0);
00236 
00237   assert(check(L0, U0, ILU0, NumGlobalPoints, NumMyPoints, 0, verbose)==0);
00238 
00239   if (verbose) cout << "\n\n*****Testing ILU(1) constructor on A" << endl<< endl;
00240 
00241   Ifpack_IlukGraph ILU1(A, 1, 0);
00242   assert(ILU1.ConstructFilledGraph()==0);
00243 
00244   assert(check(L1, U1, ILU1, NumGlobalPoints, NumMyPoints, 1, verbose)==0);
00245 
00246   if (verbose) cout << "\n\n*****Testing ILU(2) constructor on A" << endl<< endl;
00247 
00248   Ifpack_IlukGraph ILU2(A, 2, 0);
00249   assert(ILU2.ConstructFilledGraph()==0);
00250 
00251   assert(check(L2, U2, ILU2, NumGlobalPoints, NumMyPoints, 2, verbose)==0);
00252 
00253   if (verbose) cout << "\n\n*****Testing copy constructor" << endl<< endl;
00254 
00255   Ifpack_IlukGraph ILUC(ILU2);
00256   
00257   assert(check(L2, U2, ILUC, NumGlobalPoints, NumMyPoints, 2, verbose)==0);
00258 
00259   if (verbose) cout << "\n\n*****Testing copy constructor" << endl<< endl;
00260 
00261   Ifpack_IlukGraph * OverlapGraph;
00262   for (int overlap = 1; overlap < 4; overlap++) {
00263     if (verbose) cout << "\n\n*********************************************" << endl;
00264     if (verbose) cout << "\n\nConstruct Level 1 fill with Overlap = " << overlap << ".\n\n" << endl;
00265     
00266     OverlapGraph = new Ifpack_IlukGraph(A, 1, overlap);
00267     assert(OverlapGraph->ConstructFilledGraph()==0);
00268     
00269     
00270     if (verbose) {
00271       cout << "Number of Global Rows     = " << OverlapGraph->NumGlobalRows() << endl;
00272       cout << "Number of Global Nonzeros = " << OverlapGraph->NumGlobalNonzeros() << endl;
00273       cout << "Number of Local Rows     = " << OverlapGraph->NumMyRows() << endl;
00274       cout << "Number of Local Nonzeros = " << OverlapGraph->NumMyNonzeros() << endl;
00275     }
00276     delete OverlapGraph;
00277   }
00278     
00279   if (verbose1) {
00280     // Test ostream << operator (if verbose1)
00281     // Construct a Map that puts 6 equations on each PE
00282     
00283     int NumElements1 = 6;
00284     int NumPoints1 = NumElements1;
00285 
00286     // Create an integer vector NumNz that is used to build the Petra Matrix.
00287     // NumNz[i] is the Number of terms for the ith global equation on this processor
00288     
00289     int * NumNz1 = new int[NumPoints1];
00290     
00291     // We are building a tridiagonal matrix where each row has (-1 2 -1)
00292     // So we need 2 off-diagonal terms (except for the first and last equation)
00293     
00294     for (i=0; i<NumPoints1; i++)
00295       if (i==0 || i == NumPoints1-1)
00296   NumNz1[i] = 2;
00297       else
00298   NumNz1[i] = 3;
00299     
00300     // Create a Epetra_Matrix
00301     
00302     Epetra_Map Map1(NumPoints1, NumPoints1, 1, Comm);
00303     Epetra_CrsGraph A1(Copy, Map1, NumNz1);
00304     
00305     // Add  rows one-at-a-time
00306     // Need some vectors to help
00307     // Off diagonal Values will always be -1
00308     
00309     
00310     int *Indices1 = new int[2];
00311     int NumEntries1;
00312     
00313     for (i=0; i<NumPoints1; i++)
00314       {
00315   if (i==0)
00316     {
00317       Indices1[0] = 2;
00318       NumEntries1 = 1;
00319     }
00320   else if (i == NumPoints1-1)
00321     {
00322       Indices1[0] = NumPoints1-1;
00323       NumEntries1 = 1;
00324     }
00325   else
00326     {
00327       Indices1[0] = i;
00328       Indices1[1] = i+2;
00329       NumEntries1 = 2;
00330     }
00331   assert(A1.InsertGlobalIndices(i+1, NumEntries1, Indices1)==0);
00332   int ip1 = i+1;
00333   assert(A1.InsertGlobalIndices(ip1, 1, &ip1)==0); // Put in the diagonal entry
00334       }
00335     
00336     // Finish up
00337     assert(A1.FillComplete()==0);
00338     
00339     if (verbose) cout << "\n\nPrint out tridiagonal matrix with IndexBase = 1.\n\n" << endl;
00340     cout << A1 << endl;
00341 
00342     if (verbose) cout << "\n\nConstruct Level 1 fill with IndexBase = 1.\n\n" << endl;
00343 
00344     Ifpack_IlukGraph ILU11(A1, 1, 0);
00345     assert(ILU11.ConstructFilledGraph()==0);
00346 
00347     if (verbose) cout << "\n\nPrint out Level 1 ILU graph of tridiagonal matrix with IndexBase = 1.\n\n" << endl;
00348     if (verbose1) cout << ILU11 << endl;
00349 
00350     
00351   // Release all objects
00352   delete [] NumNz1;
00353   delete [] Indices1;
00354 
00355   }
00356 
00357 #ifdef EPETRA_MPI
00358   MPI_Finalize() ;
00359 #endif
00360 
00361 
00362 /* end main
00363 */
00364 return ierr ;
00365 }
00366 
00367 int check(Epetra_CrsGraph& L, Epetra_CrsGraph& U, Ifpack_IlukGraph& LU, 
00368     int NumGlobalRows1, int NumMyRows1, int LevelFill1, bool verbose) {  
00369 
00370   int i, j;
00371   int NumIndices, * Indices;
00372   int NumIndices1, * Indices1;
00373 
00374   bool debug = true;
00375 
00376   Epetra_CrsGraph& L1 = LU.L_Graph();
00377   Epetra_CrsGraph& U1 = LU.U_Graph();
00378 
00379   // Test entries and count nonzeros
00380 
00381   int Nout = 0;
00382 
00383   for (i=0; i<LU.NumMyRows(); i++) {
00384 
00385     assert(L.ExtractMyRowView(i, NumIndices, Indices)==0);
00386     assert(L1.ExtractMyRowView(i, NumIndices1, Indices1)==0);
00387     assert(NumIndices==NumIndices1);
00388     for (j=0; j<NumIndices1; j++) {
00389       if (debug &&(Indices[j]!=Indices1[j])) {
00390   int MyPID = L.RowMap().Comm().MyPID();
00391   cout << "Proc " << MyPID
00392        << " Local Row = " << i
00393        << "  L.Indices["<< j <<"]  = " << Indices[j]
00394        << " L1.Indices["<< j <<"] = " << Indices1[j] << endl;
00395       }
00396       assert(Indices[j]==Indices1[j]);
00397     }
00398     Nout += (NumIndices-NumIndices1);
00399 
00400     assert(U.ExtractMyRowView(i, NumIndices, Indices)==0);
00401     assert(U1.ExtractMyRowView(i, NumIndices1, Indices1)==0);
00402     assert(NumIndices==NumIndices1);
00403     for (j=0; j<NumIndices1; j++)  {
00404       if (debug &&(Indices[j]!=Indices1[j])) {
00405   int MyPID = L.RowMap().Comm().MyPID();
00406   cout << "Proc " << MyPID
00407        << " Local Row = " << i
00408        << "  U.Indices["<< j <<"]  = " << Indices[j]
00409        << " U1.Indices["<< j <<"] = " << Indices1[j] << endl;
00410       }
00411       assert(Indices[j]==Indices1[j]);
00412     }
00413     Nout += (NumIndices-NumIndices1);
00414   }
00415 
00416   // Test query functions
00417 
00418   int NumGlobalRows = LU.NumGlobalRows();
00419   if (verbose) cout << "\n\nNumber of Global Rows = " << NumGlobalRows << endl<< endl;
00420 
00421   assert(NumGlobalRows==NumGlobalRows1);
00422 
00423   int NumGlobalNonzeros = LU.NumGlobalNonzeros();
00424   if (verbose) cout << "\n\nNumber of Global Nonzero entries = " 
00425         << NumGlobalNonzeros << endl<< endl;
00426 
00427   int NoutG = 0;
00428 
00429   L.RowMap().Comm().SumAll(&Nout, &NoutG, 1);
00430 
00431   assert(NumGlobalNonzeros==L.NumGlobalNonzeros()+U.NumGlobalNonzeros()-NoutG);
00432 
00433   int NumMyRows = LU.NumMyRows();
00434   if (verbose) cout << "\n\nNumber of Rows = " << NumMyRows << endl<< endl;
00435 
00436   assert(NumMyRows==NumMyRows1);
00437 
00438   int NumMyNonzeros = LU.NumMyNonzeros();
00439   if (verbose) cout << "\n\nNumber of Nonzero entries = " << NumMyNonzeros << endl<< endl;
00440 
00441   assert(NumMyNonzeros==L.NumMyNonzeros()+U.NumMyNonzeros()-Nout);
00442 
00443   if (verbose) cout << "\n\nLU check OK" << endl<< endl;
00444 
00445   return(0);
00446 }

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