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