Epetra Package Browser (Single Doxygen Collection) Development
test/Bugs_LL/Bug_5797_Sort_MakeColMap_LL/cxx_main.cpp
Go to the documentation of this file.
00001 // Program for testing Epetra64 implementation.
00002 // Builds a Laplacian matrx using either Epetra or Epetra64.
00003 // Multiplies it by the vector of all ones and checks norms.
00004 //
00005 // To run:
00006 //    mpirun -np # test64.exe [n]
00007 // where n is an optional argument specifying the number of matrix rows.
00008 // Default n == 25.
00009 //
00010 // Two macro definitions below control behavior:
00011 //    ITYPE:  int --> use Epetra
00012 //            long long --> use Epetra64
00013 //    OFFSET_EPETRA64:  Add this value to each row/column index.  Resulting
00014 //                      matrix rows/columns are indexed from 
00015 //                      OFFSET_EPETRA64 to OFFSET_EPETRA64+n-1.
00016 
00017 #include <limits.h>
00018 #define ITYPE long long
00019 //#define OFFSET_EPETRA64 ((long long)(2)*(long long)INT_MAX)
00020 //#define OFFSET_EPETRA64 (INT_MAX-5)
00021 #define OFFSET_EPETRA64 0
00022 
00023 #include <stdio.h>
00024 
00025 #include "Epetra_ConfigDefs.h"
00026 
00027 #ifdef EPETRA_MPI
00028 #include <mpi.h>
00029 #include "Epetra_MpiComm.h"
00030 #define FINALIZE MPI_Finalize()
00031 #else
00032 #include "Epetra_SerialComm.h"
00033 #define FINALIZE
00034 #endif
00035  
00036 #include "Epetra_CrsMatrix.h"
00037 #include "Epetra_Map.h"
00038 #include "Epetra_Vector.h"
00039 
00041 // Build the following Laplacian matrix using Epetra or Epetra64
00042 //
00043 //    | 2 -1  ...    -1...           |
00044 //    |-1  3 -1 ...     -1...        |
00045 //    |        ...                   |
00046 //    |-1       -1  4 -1...   -1...  |
00047 //    |        ...                   |
00048 //    |         -1            -1 3 -1|
00049 //    |            -1           -1  2|
00050 //
00052 
00053 #define MIN(a,b) ((a) < (b) ? (a) : (b))
00054 
00055 int main(int narg, char *arg[])
00056 {
00057   using std::cout;
00058 
00059 #ifdef EPETRA_MPI  
00060   // Initialize MPI  
00061   MPI_Init(&narg,&arg);
00062   Epetra_MpiComm comm(MPI_COMM_WORLD);
00063 #else
00064   Epetra_SerialComm comm;
00065 #endif
00066 
00067   int me = comm.MyPID();
00068   int np = comm.NumProc();
00069 
00070   ITYPE nGlobalRows = 25;
00071   if (narg > 1) 
00072     nGlobalRows = (ITYPE) atol(arg[1]);
00073 
00074   bool verbose = (nGlobalRows < 26);
00075 
00076   // Linear map similar to Trilinos default, 
00077   // but want to allow adding OFFSET_EPETRA64 to the indices.
00078   int nMyRows = (int) (nGlobalRows / np + (nGlobalRows % np > me));
00079   ITYPE myFirstRow = (ITYPE)(me * (nGlobalRows / np) + MIN(nGlobalRows%np, me));
00080   ITYPE *myGlobalRows = new ITYPE[nMyRows];
00081   for (int i = 0; i < nMyRows; i++)
00082     myGlobalRows[i] = (ITYPE)i + myFirstRow + OFFSET_EPETRA64;
00083   Epetra_Map *rowMap = new Epetra_Map(-1, nMyRows, &myGlobalRows[0], 0, comm);
00084   if (verbose) rowMap->Print(std::cout);
00085 
00086   // Create an integer vector nnzPerRow that is used to build the Epetra Matrix.
00087   // nnzPerRow[i] is the number of entries for the ith local equation
00088   std::vector<int> nnzPerRow(nMyRows+1, 0);
00089 
00090   // Also create lists of the nonzeros to be assigned to processors.
00091   // To save programming time and complexity, these vectors are allocated 
00092   // bigger than they may actually be needed.
00093   std::vector<ITYPE> iv(5*nMyRows+1);
00094   std::vector<ITYPE> jv(5*nMyRows+1);
00095   std::vector<double> vv(5*nMyRows+1);
00096 
00097   // Generate the nonzeros for the Laplacian matrix.
00098   ITYPE nMyNonzeros = 0;
00099   for (ITYPE i = 0, myrowcnt = 0; i < nGlobalRows; i++) {
00100     if (rowMap->MyGID(i+OFFSET_EPETRA64)) { 
00101       int idegree = 0;
00102       // This processor owns this row; add nonzeros.
00103       if (i > 0) {
00104         iv[nMyNonzeros] = i + OFFSET_EPETRA64;
00105         jv[nMyNonzeros] = i-1 + OFFSET_EPETRA64;
00106         vv[nMyNonzeros] = -1;
00107         if (verbose)
00108           std::cout << "(" << iv[nMyNonzeros] << "," << jv[nMyNonzeros] << ")="
00109                << vv[nMyNonzeros] << " on processor " << me
00110                << " in " << myrowcnt << std::endl;
00111         nMyNonzeros++;
00112         nnzPerRow[myrowcnt]++;
00113         idegree++;
00114       }
00115 
00116       if (i < nGlobalRows - 1) {
00117         iv[nMyNonzeros] = i + OFFSET_EPETRA64;
00118         jv[nMyNonzeros] = i+1 + OFFSET_EPETRA64;
00119         vv[nMyNonzeros] = -1;
00120         if (verbose) 
00121           std::cout << "(" << iv[nMyNonzeros] << "," << jv[nMyNonzeros] << ")="
00122                << vv[nMyNonzeros] << " on processor " << me
00123                << " in " << myrowcnt << std::endl;
00124         nMyNonzeros++;
00125         nnzPerRow[myrowcnt]++;
00126         idegree++;
00127       }
00128 
00129       if (i+5 < nGlobalRows) {
00130         iv[nMyNonzeros] = i + OFFSET_EPETRA64;
00131         jv[nMyNonzeros] = i+5 + OFFSET_EPETRA64;
00132         vv[nMyNonzeros] = -1;
00133         if (verbose)
00134           std::cout << "(" << iv[nMyNonzeros] << "," << jv[nMyNonzeros] << ")="
00135                << vv[nMyNonzeros] << " on processor " << me
00136                << " in " << myrowcnt << std::endl;
00137         nMyNonzeros++;
00138         nnzPerRow[myrowcnt]++;
00139         idegree++;
00140       }
00141 
00142       if (i-5 >= 0) {
00143         iv[nMyNonzeros] = i + OFFSET_EPETRA64;
00144         jv[nMyNonzeros] = i-5 + OFFSET_EPETRA64;
00145         vv[nMyNonzeros] = -1;
00146         if (verbose)
00147           std::cout << "(" << iv[nMyNonzeros] << "," << jv[nMyNonzeros] << ")="
00148                << vv[nMyNonzeros] << " on processor " << me
00149                << " in " << myrowcnt << std::endl;
00150         nMyNonzeros++;
00151         nnzPerRow[myrowcnt]++;
00152         idegree++;
00153       }
00154 
00155       iv[nMyNonzeros] = i + OFFSET_EPETRA64;
00156       jv[nMyNonzeros] = i + OFFSET_EPETRA64;
00157       vv[nMyNonzeros] = idegree;
00158       if (verbose) 
00159         std::cout << "(" << iv[nMyNonzeros] << "," << jv[nMyNonzeros] << ")="
00160              << vv[nMyNonzeros] << " on processor " << me
00161              << " in " << myrowcnt << std::endl;
00162       nMyNonzeros++;
00163       nnzPerRow[myrowcnt]++;
00164 
00165       myrowcnt++;
00166     }
00167   }
00168 
00169   // Create an Epetra_Matrix
00170   // TODO:  Should be able to use StaticProfile == true, but there appears to
00171   // TODO:  be a bug in Epetra64 (or in my usage of it).  I filed bug 5791
00172   // TODO:  and will use StaticProfile == false for now.  When the bug is fixed,
00173   // TODO:  we should switch back to StaticProfile == true, as it will be more
00174   // TODO:  efficient.
00175   // #warning Once bug 5791 is fixed, change StaticProfile back to true.
00176   // -- Done.  Chetan.
00177   Epetra_CrsMatrix *A = new Epetra_CrsMatrix(Copy, *rowMap, &nnzPerRow[0], true);
00178 
00179   // Insert the nonzeros.
00180   int info;
00181   ITYPE sum = 0;
00182   for (int i=0; i < nMyRows; i++) {
00183     if (nnzPerRow[i]) {
00184       if (verbose) {
00185         std::cout << "InsertGlobalValus row " << iv[sum]
00186              << " count " << nnzPerRow[i] 
00187              << " cols " << jv[sum] << " " << jv[sum+1] << " ";
00188         if (nnzPerRow[i] == 3) std::cout << jv[sum+2];
00189         std::cout << std::endl;
00190       }
00191       info = A->InsertGlobalValues(iv[sum],nnzPerRow[i],&vv[sum],&jv[sum]);
00192       assert(info==0);
00193       sum += nnzPerRow[i];
00194     }
00195   }
00196 
00197   // Finish up
00198   info = A->FillComplete();
00199   assert(info==0);
00200   if (verbose) A->Print(std::cout);
00201 
00202   // Sanity test:  Product of matrix and vector of ones should have norm == 0
00203   // and max/min/mean values of 0
00204   Epetra_Vector sanity(A->RangeMap());
00205   Epetra_Vector sanityres(A->DomainMap());
00206   sanity.PutScalar(1.);
00207   A->Multiply(false, sanity, sanityres);
00208 
00209   double jjone, jjtwo, jjmax;
00210   sanityres.Norm1(&jjone);
00211   sanityres.Norm2(&jjtwo);
00212   sanityres.NormInf(&jjmax);
00213   if (me == 0)
00214     std::cout << "SanityTest norms 1/2/inf: " << jjone << " "
00215                                          << jjtwo << " " << jjmax << std::endl;
00216 
00217   bool test_failed = (jjone != 0) || (jjtwo != 0) || (jjmax != 0);
00218 
00219   sanityres.MinValue(&jjone);
00220   sanityres.MeanValue(&jjtwo);
00221   sanityres.MaxValue(&jjmax);
00222   if (me == 0)
00223     std::cout << "SanityTest values min/max/avg: " << jjone << " "
00224                                               << jjmax << " " << jjtwo << std::endl;
00225 
00226   test_failed = test_failed || (jjone != 0) || (jjtwo != 0) || (jjmax != 0);
00227 
00228   if (me == 0) {
00229     if(test_failed)
00230       std::cout << "Bug_5797_Sort_MakeColMap_LL tests FAILED" << std::endl;
00231   }
00232 
00233   delete A;
00234   delete rowMap;
00235   delete [] myGlobalRows;
00236 
00237   FINALIZE;
00238 }
00239 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines