Ifpack Package Browser (Single Doxygen Collection) Development
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 #include "Teuchos_RefCountPtr.hpp"
00046 #ifdef EPETRA_MPI
00047 #include "Epetra_MpiComm.h"
00048 #else
00049 #include "Epetra_SerialComm.h"
00050 #endif
00051 
00052 // Prototype
00053 
00054 int check(Epetra_CrsGraph& L, Epetra_CrsGraph& U, Ifpack_IlukGraph& LU, 
00055     int NumGlobalRows1, int NumMyRows1, int LevelFill1, bool verbose);
00056 
00057  int main(int argc, char *argv[])
00058 {
00059   int ierr = 0, i, j;
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   std::vector<int> Indices(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])>=0);
00165   assert(L0.InsertGlobalIndices(Row, k, &Indices[0])>=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])>=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])>=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])>=0);
00195   assert(U0.InsertGlobalIndices(Row, k, &Indices[0])>=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])>=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])>=0);
00210       }
00211     }
00212   }
00213 
00214   // Finish up
00215   if (verbose) cout << "\n\nCompleting A" << endl<< endl;
00216   assert(A.FillComplete()==0);
00217   if (verbose) cout << "\n\nCompleting L0" << endl<< endl;
00218   assert(L0.FillComplete()==0);
00219   if (verbose) cout << "\n\nCompleting U0" << endl<< endl;
00220   assert(U0.FillComplete()==0);
00221   if (verbose) cout << "\n\nCompleting L1" << endl<< endl;
00222   assert(L1.FillComplete()==0);
00223   if (verbose) cout << "\n\nCompleting U1" << endl<< endl;
00224   assert(U1.FillComplete()==0);
00225   if (verbose) cout << "\n\nCompleting L2" << endl<< endl;
00226   assert(L2.FillComplete()==0);
00227   if (verbose) cout << "\n\nCompleting U2" << endl<< endl;
00228   assert(U2.FillComplete()==0);
00229 
00230   if (verbose) cout << "\n\n*****Testing ILU(0) constructor on A" << endl<< endl;
00231 
00232   Ifpack_IlukGraph ILU0(A, 0, 0);
00233   assert(ILU0.ConstructFilledGraph()==0);
00234 
00235   assert(check(L0, U0, ILU0, NumGlobalPoints, NumMyPoints, 0, verbose)==0);
00236 
00237   if (verbose) cout << "\n\n*****Testing ILU(1) constructor on A" << endl<< endl;
00238 
00239   Ifpack_IlukGraph ILU1(A, 1, 0);
00240   assert(ILU1.ConstructFilledGraph()==0);
00241 
00242   assert(check(L1, U1, ILU1, NumGlobalPoints, NumMyPoints, 1, verbose)==0);
00243 
00244   if (verbose) cout << "\n\n*****Testing ILU(2) constructor on A" << endl<< endl;
00245 
00246   Ifpack_IlukGraph ILU2(A, 2, 0);
00247   assert(ILU2.ConstructFilledGraph()==0);
00248 
00249   assert(check(L2, U2, ILU2, NumGlobalPoints, NumMyPoints, 2, verbose)==0);
00250 
00251   if (verbose) cout << "\n\n*****Testing copy constructor" << endl<< endl;
00252 
00253   Ifpack_IlukGraph ILUC(ILU2);
00254   
00255   assert(check(L2, U2, ILUC, NumGlobalPoints, NumMyPoints, 2, verbose)==0);
00256 
00257   if (verbose) cout << "\n\n*****Testing copy constructor" << endl<< endl;
00258 
00259   Teuchos::RefCountPtr<Ifpack_IlukGraph> OverlapGraph;
00260   for (int overlap = 1; overlap < 4; overlap++) {
00261     if (verbose) cout << "\n\n*********************************************" << endl;
00262     if (verbose) cout << "\n\nConstruct Level 1 fill with Overlap = " << overlap << ".\n\n" << endl;
00263     
00264     OverlapGraph = Teuchos::rcp( new Ifpack_IlukGraph(A, 1, overlap) );
00265     assert(OverlapGraph->ConstructFilledGraph()==0);
00266     
00267     if (verbose) {
00268       cout << "Number of Global Rows     = " << OverlapGraph->NumGlobalRows() << endl;
00269       cout << "Number of Global Nonzeros = " << OverlapGraph->NumGlobalNonzeros() << endl;
00270       cout << "Number of Local Rows     = " << OverlapGraph->NumMyRows() << endl;
00271       cout << "Number of Local Nonzeros = " << OverlapGraph->NumMyNonzeros() << endl;
00272     }
00273   }
00274     
00275   if (verbose1) {
00276     // Test ostream << operator (if verbose1)
00277     // Construct a Map that puts 6 equations on each PE
00278     
00279     int NumElements1 = 6;
00280     int NumPoints1 = NumElements1;
00281 
00282     // Create an integer vector NumNz that is used to build the Petra Matrix.
00283     // NumNz[i] is the Number of terms for the ith global equation on this processor
00284     
00285     std::vector<int> NumNz1(NumPoints1);
00286     
00287     // We are building a tridiagonal matrix where each row has (-1 2 -1)
00288     // So we need 2 off-diagonal terms (except for the first and last equation)
00289     
00290     for (i=0; i<NumPoints1; i++)
00291       if (i==0 || i == NumPoints1-1)
00292   NumNz1[i] = 2;
00293       else
00294   NumNz1[i] = 3;
00295     
00296     // Create a Epetra_Matrix
00297     
00298     Epetra_Map Map1(NumPoints1, NumPoints1, 1, Comm);
00299     Epetra_CrsGraph A1(Copy, Map1, &NumNz1[0]);
00300     
00301     // Add  rows one-at-a-time
00302     // Need some vectors to help
00303     // Off diagonal Values will always be -1
00304     
00305     
00306     std::vector<int> Indices1(2);
00307     int NumEntries1;
00308     
00309     for (i=0; i<NumPoints1; i++)
00310       {
00311   if (i==0)
00312     {
00313       Indices1[0] = 2;
00314       NumEntries1 = 1;
00315     }
00316   else if (i == NumPoints1-1)
00317     {
00318       Indices1[0] = NumPoints1-1;
00319       NumEntries1 = 1;
00320     }
00321   else
00322     {
00323       Indices1[0] = i;
00324       Indices1[1] = i+2;
00325       NumEntries1 = 2;
00326     }
00327   assert(A1.InsertGlobalIndices(i+1, NumEntries1, &Indices1[0])==0);
00328   int ip1 = i+1;
00329   assert(A1.InsertGlobalIndices(ip1, 1, &ip1)==0); // Put in the diagonal entry
00330       }
00331     
00332     // Finish up
00333     assert(A1.FillComplete()==0);
00334     
00335     if (verbose) cout << "\n\nPrint out tridiagonal matrix with IndexBase = 1.\n\n" << endl;
00336     cout << A1 << endl;
00337 
00338     if (verbose) cout << "\n\nConstruct Level 1 fill with IndexBase = 1.\n\n" << endl;
00339 
00340     Ifpack_IlukGraph ILU11(A1, 1, 0);
00341     assert(ILU11.ConstructFilledGraph()==0);
00342 
00343     if (verbose) cout << "\n\nPrint out Level 1 ILU graph of tridiagonal matrix with IndexBase = 1.\n\n" << endl;
00344     if (verbose1) cout << ILU11 << endl;
00345 
00346   }
00347 
00348 #ifdef EPETRA_MPI
00349   MPI_Finalize() ;
00350 #endif
00351 
00352 
00353 /* end main
00354 */
00355 return ierr ;
00356 }
00357 
00358 int check(Epetra_CrsGraph& L, Epetra_CrsGraph& U, Ifpack_IlukGraph& LU, 
00359     int NumGlobalRows1, int NumMyRows1, int LevelFill1, bool verbose) {  
00360 
00361   int i, j;
00362   int NumIndices, * Indices;
00363   int NumIndices1, * Indices1;
00364 
00365   bool debug = true;
00366 
00367   Epetra_CrsGraph& L1 = LU.L_Graph();
00368   Epetra_CrsGraph& U1 = LU.U_Graph();
00369 
00370   // Test entries and count nonzeros
00371 
00372   int Nout = 0;
00373 
00374   for (i=0; i<LU.NumMyRows(); i++) {
00375 
00376     assert(L.ExtractMyRowView(i, NumIndices, Indices)==0);
00377     assert(L1.ExtractMyRowView(i, NumIndices1, Indices1)==0);
00378     assert(NumIndices==NumIndices1);
00379     for (j=0; j<NumIndices1; j++) {
00380       if (debug &&(Indices[j]!=Indices1[j])) {
00381   int MyPID = L.RowMap().Comm().MyPID();
00382   cout << "Proc " << MyPID
00383        << " Local Row = " << i
00384        << "  L.Indices["<< j <<"]  = " << Indices[j]
00385        << " L1.Indices["<< j <<"] = " << Indices1[j] << endl;
00386       }
00387       assert(Indices[j]==Indices1[j]);
00388     }
00389     Nout += (NumIndices-NumIndices1);
00390 
00391     assert(U.ExtractMyRowView(i, NumIndices, Indices)==0);
00392     assert(U1.ExtractMyRowView(i, NumIndices1, Indices1)==0);
00393     assert(NumIndices==NumIndices1);
00394     for (j=0; j<NumIndices1; j++)  {
00395       if (debug &&(Indices[j]!=Indices1[j])) {
00396   int MyPID = L.RowMap().Comm().MyPID();
00397   cout << "Proc " << MyPID
00398        << " Local Row = " << i
00399        << "  U.Indices["<< j <<"]  = " << Indices[j]
00400        << " U1.Indices["<< j <<"] = " << Indices1[j] << endl;
00401       }
00402       assert(Indices[j]==Indices1[j]);
00403     }
00404     Nout += (NumIndices-NumIndices1);
00405   }
00406 
00407   // Test query functions
00408 
00409   int NumGlobalRows = LU.NumGlobalRows();
00410   if (verbose) cout << "\n\nNumber of Global Rows = " << NumGlobalRows << endl<< endl;
00411 
00412   assert(NumGlobalRows==NumGlobalRows1);
00413 
00414   int NumGlobalNonzeros = LU.NumGlobalNonzeros();
00415   if (verbose) cout << "\n\nNumber of Global Nonzero entries = " 
00416         << NumGlobalNonzeros << endl<< endl;
00417 
00418   int NoutG = 0;
00419 
00420   L.RowMap().Comm().SumAll(&Nout, &NoutG, 1);
00421 
00422   assert(NumGlobalNonzeros==L.NumGlobalNonzeros()+U.NumGlobalNonzeros()-NoutG);
00423 
00424   int NumMyRows = LU.NumMyRows();
00425   if (verbose) cout << "\n\nNumber of Rows = " << NumMyRows << endl<< endl;
00426 
00427   assert(NumMyRows==NumMyRows1);
00428 
00429   int NumMyNonzeros = LU.NumMyNonzeros();
00430   if (verbose) cout << "\n\nNumber of Nonzero entries = " << NumMyNonzeros << endl<< endl;
00431 
00432   assert(NumMyNonzeros==L.NumMyNonzeros()+U.NumMyNonzeros()-Nout);
00433 
00434   if (verbose) cout << "\n\nLU check OK" << endl<< endl;
00435 
00436   return(0);
00437 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines