Epetra Package Browser (Single Doxygen Collection) Development
test/Bugs_LL/Bug_5794_IndexBase_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 == 10.
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 
00022 #include <stdio.h>
00023 
00024 #include "Epetra_ConfigDefs.h"
00025 
00026 #ifdef EPETRA_MPI
00027 #include <mpi.h>
00028 #include "Epetra_MpiComm.h"
00029 #define FINALIZE MPI_Finalize()
00030 #else
00031 #include "Epetra_SerialComm.h"
00032 #define FINALIZE
00033 #endif
00034 
00035 #include "Epetra_CrsMatrix.h"
00036 #include "Epetra_Map.h"
00037 #include "Epetra_Vector.h"
00038 
00040 // Build the following Laplacian matrix using Epetra or Epetra64.
00041 //
00042 //    | 1 -1          |
00043 //    |-1  2 -1       |
00044 //    |   -1  2 -1    |
00045 //    |        ...    |
00046 //    |           -1 1|
00047 //
00049 
00050 #define MIN(a,b) ((a) < (b) ? (a) : (b))
00051 
00052 int main(int narg, char *arg[])
00053 {
00054   using std::cout;
00055 
00056 #ifdef EPETRA_MPI
00057   // Initialize MPI
00058   MPI_Init(&narg,&arg);
00059   Epetra_MpiComm comm(MPI_COMM_WORLD);
00060 #else
00061   Epetra_SerialComm comm;
00062 #endif
00063 
00064   int me = comm.MyPID();
00065   int np = comm.NumProc();
00066 
00067   ITYPE nGlobalRows = 10;
00068   if (narg > 1)
00069     nGlobalRows = (ITYPE) atol(arg[1]);
00070 
00071   bool verbose = (nGlobalRows < 20);
00072 
00073   // Linear map similar to Trilinos default,
00074   // but want to allow adding OFFSET_EPETRA64 to the indices.
00075   int nMyRows = (int) (nGlobalRows / np + (nGlobalRows % np > me));
00076   ITYPE myFirstRow = (ITYPE)(me * (nGlobalRows / np) + MIN(nGlobalRows%np, me));
00077   ITYPE *myGlobalRows = new ITYPE[nMyRows];
00078   for (int i = 0; i < nMyRows; i++)
00079     myGlobalRows[i] = (ITYPE)i + myFirstRow + OFFSET_EPETRA64;
00080   Epetra_Map *rowMap = new Epetra_Map(-1, nMyRows, &myGlobalRows[0], 0, comm);
00081   if (verbose) rowMap->Print(std::cout);
00082 
00083   // Create an integer vector nnzPerRow that is used to build the Epetra Matrix.
00084   // nnzPerRow[i] is the number of entries for the ith local equation
00085   std::vector<int> nnzPerRow(nMyRows+1, 0);
00086 
00087   // Also create lists of the nonzeros to be assigned to processors.
00088   // To save programming time and complexity, these vectors are allocated
00089   // bigger than they may actually be needed.
00090   std::vector<ITYPE> iv(3*nMyRows+1);
00091   std::vector<ITYPE> jv(3*nMyRows+1);
00092   std::vector<double> vv(3*nMyRows+1);
00093 
00094   // Generate the nonzeros for the Laplacian matrix.
00095   ITYPE nMyNonzeros = 0;
00096   for (ITYPE i = 0, myrowcnt = 0; i < nGlobalRows; i++) {
00097     if (rowMap->MyGID(i+OFFSET_EPETRA64)) {
00098       // This processor owns this row; add nonzeros.
00099       if (i > 0) {
00100         iv[nMyNonzeros] = i + OFFSET_EPETRA64;
00101         jv[nMyNonzeros] = i-1 + OFFSET_EPETRA64;
00102         vv[nMyNonzeros] = -1;
00103         if (verbose)
00104           std::cout << "(" << iv[nMyNonzeros] << "," << jv[nMyNonzeros] << ")="
00105                << vv[nMyNonzeros] << " on processor " << me
00106                << " in " << myrowcnt << std::endl;
00107         nMyNonzeros++;
00108         nnzPerRow[myrowcnt]++;
00109       }
00110 
00111       iv[nMyNonzeros] = i + OFFSET_EPETRA64;
00112       jv[nMyNonzeros] = i + OFFSET_EPETRA64;
00113       vv[nMyNonzeros] = ((i == 0 || i == nGlobalRows-1) ? 1. : 2.);
00114       if (verbose)
00115         std::cout << "(" << iv[nMyNonzeros] << "," << jv[nMyNonzeros] << ")="
00116              << vv[nMyNonzeros] << " on processor " << me
00117              << " in " << myrowcnt << std::endl;
00118       nMyNonzeros++;
00119       nnzPerRow[myrowcnt]++;
00120 
00121       if (i < nGlobalRows - 1) {
00122         iv[nMyNonzeros] = i + OFFSET_EPETRA64;
00123         jv[nMyNonzeros] = i+1 + OFFSET_EPETRA64;
00124         vv[nMyNonzeros] = -1;
00125         if (verbose)
00126           std::cout << "(" << iv[nMyNonzeros] << "," << jv[nMyNonzeros] << ")="
00127                << vv[nMyNonzeros] << " on processor " << me
00128                << " in " << myrowcnt << std::endl;
00129         nMyNonzeros++;
00130         nnzPerRow[myrowcnt]++;
00131       }
00132       myrowcnt++;
00133     }
00134   }
00135 
00136   // Create an Epetra_Matrix
00137   Epetra_CrsMatrix *A = new Epetra_CrsMatrix(Copy, *rowMap, &nnzPerRow[0], false);
00138 
00139   // Insert the nonzeros.
00140   int info;
00141   ITYPE sum = 0;
00142   for (int i=0; i < nMyRows; i++) {
00143     if (nnzPerRow[i]) {
00144       if (verbose) {
00145         std::cout << "InsertGlobalValus row " << iv[sum]
00146              << " count " << nnzPerRow[i]
00147              << " cols " << jv[sum] << " " << jv[sum+1] << " ";
00148         if (nnzPerRow[i] == 3) std::cout << jv[sum+2];
00149         std::cout << std::endl;
00150       }
00151       info = A->InsertGlobalValues(iv[sum],nnzPerRow[i],&vv[sum],&jv[sum]);
00152       assert(info==0);
00153       sum += nnzPerRow[i];
00154     }
00155   }
00156 
00157   // Finish up
00158   info = A->FillComplete();
00159   assert(info==0);
00160   if (verbose) A->Print(std::cout);
00161 
00162   // Sanity test:  Product of matrix and vector of ones should have norm == 0
00163   // and max/min/mean values of 0
00164   Epetra_Vector sanity(A->RangeMap());
00165   Epetra_Vector sanityres(A->DomainMap());
00166   sanity.PutScalar(1.);
00167   A->Multiply(false, sanity, sanityres);
00168 
00169   double jjone, jjtwo, jjmax;
00170   sanityres.Norm1(&jjone);
00171   sanityres.Norm2(&jjtwo);
00172   sanityres.NormInf(&jjmax);
00173   if (me == 0)
00174     std::cout << "SanityTest norms 1/2/inf: " << jjone << " "
00175                                          << jjtwo << " " << jjmax << std::endl;
00176 
00177   bool test_failed = (jjone != 0) || (jjtwo != 0) || (jjmax != 0);
00178 
00179   sanityres.MinValue(&jjone);
00180   sanityres.MeanValue(&jjtwo);
00181   sanityres.MaxValue(&jjmax);
00182   if (me == 0)
00183     std::cout << "SanityTest values min/max/avg: " << jjone << " "
00184                                               << jjmax << " " << jjtwo << std::endl;
00185 
00186   test_failed = test_failed || (jjone != 0) || (jjtwo != 0) || (jjmax != 0);
00187 
00188   if (me == 0) {
00189     if(test_failed)
00190       std::cout << "Bug_5794_IndexBase_LL tests FAILED" << std::endl;
00191   }
00192 
00193   delete A;
00194   delete rowMap;
00195   delete [] myGlobalRows;
00196 
00197   FINALIZE;
00198 }
00199 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines