Epetra Package Browser (Single Doxygen Collection) Development
test/EpetraBenchmarkTest/cxx_main.cpp
Go to the documentation of this file.
00001 //@HEADER
00002 // ************************************************************************
00003 //
00004 //               Epetra: Linear Algebra Services Package
00005 //                 Copyright 2011 Sandia Corporation
00006 //
00007 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00008 // the U.S. Government retains certain rights in this software.
00009 //
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00038 //
00039 // ************************************************************************
00040 //@HEADER
00041 
00042 #include <iostream>
00043 #include <cstring>
00044 #include "Epetra_Map.h"
00045 #include "Epetra_LocalMap.h"
00046 #include "Epetra_BlockMap.h"
00047 #include "Epetra_Time.h"
00048 #include "Epetra_CrsMatrix.h"
00049 #include "Epetra_VbrMatrix.h"
00050 #include "Epetra_Vector.h"
00051 #include "Epetra_IntVector.h"
00052 #include "Epetra_MultiVector.h"
00053 #include "Epetra_IntSerialDenseVector.h"
00054 #include "Epetra_SerialDenseVector.h"
00055 #ifdef EPETRA_MPI
00056 #include "Epetra_MpiComm.h"
00057 #include "mpi.h"
00058 #else
00059 #include "Epetra_SerialComm.h"
00060 #endif
00061 #include "Epetra_Version.h"
00062 
00063 // prototypes
00064 
00065 void GenerateCrsProblem(int numNodesX, int numNodesY, int numProcsX, int numProcsY, int numPoints,
00066       int * xoff, int * yoff,
00067       const Epetra_Comm  &comm,
00068       Epetra_Map *& map,
00069       Epetra_CrsMatrix *& A,
00070       Epetra_Vector *& b,
00071       Epetra_Vector *& bt,
00072       Epetra_Vector *&xexact, bool StaticProfile, bool MakeLocalOnly);
00073 
00074 void GenerateCrsProblem(int numNodesX, int numNodesY, int numProcsX, int numProcsY, int numPoints,
00075       int * xoff, int * yoff, int nrhs,
00076       const Epetra_Comm  &comm,
00077       Epetra_Map *& map,
00078       Epetra_CrsMatrix *& A,
00079       Epetra_MultiVector *& b,
00080       Epetra_MultiVector *& bt,
00081       Epetra_MultiVector *&xexact, bool StaticProfile, bool MakeLocalOnly);
00082 
00083 
00084 void GenerateMyGlobalElements(int numNodesX, int numNodesY, int numProcsX, int numProcs,
00085             int myPID, int * & myGlobalElements);
00086 
00087 void runMatrixTests(Epetra_CrsMatrix * A,  Epetra_MultiVector * b, Epetra_MultiVector * bt,
00088         Epetra_MultiVector * xexact, bool StaticProfile);
00089 
00090 int main(int argc, char *argv[])
00091 {
00092 
00093   const string EpetraBenchmarkTest = "Epetra Benchmark Test Version 0.2 08/30/2007"; // Change this number and date when changes are made to functionality
00094   int ierr = 0;
00095   double elapsed_time;
00096   double total_flops;
00097   double MFLOPs;
00098   double global_dimension;
00099   double global_nonzero_count;
00100 
00101 
00102 #ifdef EPETRA_MPI
00103 
00104   // Initialize MPI
00105   MPI_Init(&argc,&argv);
00106   Epetra_MpiComm comm( MPI_COMM_WORLD );
00107 #else
00108   Epetra_SerialComm comm;
00109 #endif
00110 
00111   // Check if we should print only selected timing info to standard out
00112   bool printflops = true, printmflops = true, printtime = true;
00113   if (argc>5) if (argv[5][0]=='-' && argv[5][1]=='f') {printmflops = false; printtime = false;}
00114   bool mflopsonly = false;
00115   if (argc>5) if (argv[5][0]=='-' && argv[5][1]=='m') {printflops = false; printtime = false;}
00116   bool timeonly = false;
00117   if (argc>5) if (argv[5][0]=='-' && argv[5][1]=='t') {printflops = false; printmflops = false;}
00118 
00119   // Check if we should print header to standard out
00120   bool makeheader = false;
00121   if (argc>6) if (argv[6][0]=='-' && argv[6][1]=='h') makeheader = true;
00122 
00123   if(argc < 5) {
00124     cerr << "Usage: " << argv[0]
00125          << " NumNodesX NumNodesY NumProcX NumProcY [-a|-f|-m|-t [-h]]" << endl
00126          << "where:" << endl
00127          << "NumNodesX         - Number of mesh nodes in X direction per processor" << endl
00128          << "NumNodesY         - Number of mesh nodes in Y direction per processor" << endl
00129          << "NumProcX          - Number of processors to use in X direction" << endl
00130          << "NumProcY          - Number of processors to use in Y direction" << endl
00131    << "-a|-f|-m|-t       - Type of information to print: a=all, f=FLOPS, m=MFLOP/s, t=time (sec).  Default is -a."
00132          << "-h                - (Optional) Printer output table header if -h present (typically done on first call)." << endl
00133          << " NOTES: NumProcX*NumProcY must equal the number of processors used to run the problem." << endl << endl
00134    << " Serial example:" << endl << endl
00135          << argv[0] << " 200 300 1 1 -m" << endl
00136    << " Run this program on 1 processor using a 200 by 300 grid, printing only MFLOP/s information without a header."<< endl <<endl
00137    << " MPI example:" << endl << endl
00138          << "mpirun -np 32 " << argv[0] << " 250 200 4 8 -a -h" << endl
00139    << " Run this program on 32 processors putting a 250 by 200 subgrid on each processor using 4 processors "<< endl
00140    << " in the X direction and 8 in the Y direction.  Total grid size is 1000 points in X and 1600 in Y for a total of 1.6M equations."<< endl
00141    << " Print all information. Print header." << endl
00142          << endl;
00143     return(1);
00144 
00145   }
00146     //char tmp;
00147     //if (comm.MyPID()==0) cout << "Press any key to continue..."<< endl;
00148     //if (comm.MyPID()==0) cin >> tmp;
00149     //comm.Barrier();
00150 
00151   if (makeheader && comm.MyPID()==0)
00152     cout << EpetraBenchmarkTest << endl
00153    << "Using " << Epetra_Version() << endl << endl;
00154   if (makeheader) cout << comm <<endl;
00155 
00156 
00157   // Redefine makeheader to only print on PE 0
00158 
00159   if (makeheader && comm.MyPID()!=0) makeheader = false;
00160 
00161   int numNodesX = atoi(argv[1]);
00162   int numNodesY = atoi(argv[2]);
00163   int numProcsX = atoi(argv[3]);
00164   int numProcsY = atoi(argv[4]);
00165   int numPoints = 25;
00166 
00167   if (makeheader) {
00168     cout << " Number of local nodes in X direction  = " << numNodesX << endl
00169    << " Number of local nodes in Y direction  = " << numNodesY << endl
00170    << " Number of global nodes in X direction = " << numNodesX*numProcsX << endl
00171    << " Number of global nodes in Y direction = " << numNodesY*numProcsY << endl
00172    << " Number of local nonzero entries       = " << numNodesX*numNodesY*numPoints << endl
00173    << " Number of global nonzero entries      = " << numNodesX*numNodesY*numPoints*numProcsX*numProcsY << endl
00174    << " Number of Processors in X direction   = " << numProcsX << endl
00175    << " Number of Processors in Y direction   = " << numProcsY << endl
00176    << " Number of Points in stencil           = " << numPoints << endl << endl;
00177     cout << " Timing the following:" <<endl
00178    << " SpMV - Sparse matrix vector product using Epetra_CrsMatrix class" << endl
00179    << " SpMM2- Sparse matrix times 2-column multivector using Epetra_CrsMatrix class" << endl
00180    << " SpMM4- Sparse matrix times 4-column multivector using Epetra_CrsMatrix class" << endl
00181    << " SpMM8- Sparse matrix times 8-column multivector using Epetra_CrsMatrix class" << endl
00182    << " 2-norm of an Epetra_MultiVector" << endl
00183    << " Dot-product of 2 Epetra_MultiVectors" << endl
00184    << " AXPY of 2 Epetra_MultiVectors" << endl << endl;
00185   }
00186 
00187   if (numProcsX*numProcsY!=comm.NumProc()) {
00188     cerr << "Number of processors = " << comm.NumProc() << endl
00189    << " is not the product of " << numProcsX << " and " << numProcsY << endl << endl;
00190     return(1);
00191   }
00192 
00193   if (numNodesX*numNodesY<=0) {
00194     cerr << "Product of number of nodes is <= zero" << endl << endl;
00195     return(1);
00196   }
00197 
00198   Epetra_IntSerialDenseVector Xoff, XLoff, XUoff;
00199   Epetra_IntSerialDenseVector Yoff, YLoff, YUoff;
00200   // Generate a 25-point 2D Finite Difference matrix
00201   Xoff.Size(25);
00202   Yoff.Size(25);
00203   int xi = 0, yi = 0;
00204   int xo = -2, yo = -2;
00205   Xoff[xi++] = xo++;  Xoff[xi++] = xo++; Xoff[xi++] = xo++; Xoff[xi++] = xo++; Xoff[xi++] = xo++;
00206   Yoff[yi++] = yo  ;  Yoff[yi++] = yo  ; Yoff[yi++] = yo  ; Yoff[yi++] = yo  ; Yoff[yi++] = yo  ;
00207   xo = -2, yo++;
00208   Xoff[xi++] = xo++;  Xoff[xi++] = xo++; Xoff[xi++] = xo++; Xoff[xi++] = xo++; Xoff[xi++] = xo++;
00209   Yoff[yi++] = yo  ;  Yoff[yi++] = yo  ; Yoff[yi++] = yo  ; Yoff[yi++] = yo  ; Yoff[yi++] = yo  ;
00210   xo = -2, yo++;
00211   Xoff[xi++] = xo++;  Xoff[xi++] = xo++; Xoff[xi++] = xo++; Xoff[xi++] = xo++; Xoff[xi++] = xo++;
00212   Yoff[yi++] = yo  ;  Yoff[yi++] = yo  ; Yoff[yi++] = yo  ; Yoff[yi++] = yo  ; Yoff[yi++] = yo  ;
00213   xo = -2, yo++;
00214   Xoff[xi++] = xo++;  Xoff[xi++] = xo++; Xoff[xi++] = xo++; Xoff[xi++] = xo++; Xoff[xi++] = xo++;
00215   Yoff[yi++] = yo  ;  Yoff[yi++] = yo  ; Yoff[yi++] = yo  ; Yoff[yi++] = yo  ; Yoff[yi++] = yo  ;
00216   xo = -2, yo++;
00217   Xoff[xi++] = xo++;  Xoff[xi++] = xo++; Xoff[xi++] = xo++; Xoff[xi++] = xo++; Xoff[xi++] = xo++;
00218   Yoff[yi++] = yo  ;  Yoff[yi++] = yo  ; Yoff[yi++] = yo  ; Yoff[yi++] = yo  ; Yoff[yi++] = yo  ;
00219 
00220 
00221 
00222   Epetra_Map * map;
00223   Epetra_CrsMatrix * A;
00224   Epetra_MultiVector * b;
00225   Epetra_MultiVector * bt;
00226   Epetra_MultiVector * xexact;
00227   Epetra_SerialDenseVector resvec(0);
00228   global_dimension = (double) (numNodesX*numNodesY); // local grid dimension
00229   global_dimension *= (double) (numProcsX*numProcsY); // times number of processors
00230   global_nonzero_count = global_dimension * (double) numPoints; // number of nonzeros (w/o accounting for boundaries)
00231 
00232   //Timings
00233   Epetra_Time timer(comm);
00234   const int maxtimings = 7;
00235   const int maxtrials = 10;
00236   double results[maxtimings][3]; // timing results stored here
00237 
00238   int nrhs = 1;
00239   int ntimings = 0;
00240   for (int k=0; k<4; k++) {
00241     if (k>0) nrhs = nrhs*2;
00242 
00243     GenerateCrsProblem(numNodesX, numNodesY, numProcsX, numProcsY, numPoints,
00244            Xoff.Values(), Yoff.Values(), nrhs, comm,
00245            map, A, b, bt, xexact, true, false);
00246 
00247 
00248     Epetra_MultiVector z(*b);
00249     Epetra_MultiVector r(*b);
00250     Epetra_SerialDenseVector resvec(b->NumVectors());
00251 
00252     //Timings
00253     Epetra_Time timer(A->Comm());
00254     A->OptimizeStorage();
00255 
00256     timer.ResetStartTime();
00257 
00258     //maxtrials matvecs
00259     for( int i = 0; i < maxtrials; ++i )
00260       A->Multiply(false, *xexact, z); // Compute z = A*xexact or z = A'*xexact
00261 
00262     double elapsed_time = timer.ElapsedTime();
00263     double total_flops = 2.0*global_nonzero_count *((double) maxtrials);
00264 
00265     // Compute residual
00266     r.Update(-1.0, z, 1.0, *b, 0.0); // r = b - z
00267 
00268     r.Norm2(resvec.Values());
00269     double diff = resvec.NormInf();
00270     if (diff>1.0e-8 && comm.MyPID()==0) cerr << "Warning: Residual vector unusually large = " << diff << endl;
00271 
00272     double MFLOPs = total_flops/elapsed_time/1000000.0;
00273 
00274     results[ntimings][0] = total_flops;
00275     results[ntimings][1] = elapsed_time;
00276     results[ntimings++][2] = MFLOPs;
00277 
00278     delete A;
00279     delete b;
00280     delete bt;
00281     delete xexact;
00282   } // end of k loop
00283 
00284 
00285   // *************** Vector timings *********************
00286 
00287   if (ntimings+3>maxtimings) cerr << "Variable maxtimings = " << maxtimings << " must be at least = " << ntimings+3 << endl;
00288 
00289   Epetra_MultiVector q(*map, nrhs);
00290   Epetra_MultiVector z(q);
00291   Epetra_MultiVector r(q);
00292 
00293   delete map;
00294 
00295   resvec.Resize(nrhs);
00296 
00297 
00298   timer.ResetStartTime();
00299 
00300   //maxtrials norms
00301   for( int i = 0; i < maxtrials; ++i )
00302     q.Norm2( resvec.Values() );
00303 
00304   elapsed_time = timer.ElapsedTime();
00305   total_flops = 2.0*global_dimension *((double) maxtrials);
00306   MFLOPs = total_flops/elapsed_time/1000000.0;
00307   results[ntimings][0] = total_flops;
00308   results[ntimings][1] = elapsed_time;
00309   results[ntimings++][2] = MFLOPs;
00310 
00311 
00312 
00313   timer.ResetStartTime();
00314 
00315   //maxtrials dot's
00316   for( int i = 0; i < maxtrials; ++i )
00317     q.Dot(z, resvec.Values());
00318 
00319   elapsed_time = timer.ElapsedTime();
00320   total_flops = 2.0*global_dimension *((double) maxtrials);
00321   MFLOPs = total_flops/elapsed_time/1000000.0;
00322   results[ntimings][0] = total_flops;
00323   results[ntimings][1] = elapsed_time;
00324   results[ntimings++][2] = MFLOPs;
00325 
00326   timer.ResetStartTime();
00327 
00328   //maxtrials updates
00329   for( int i = 0; i < maxtrials; ++i )
00330     q.Update(1.0, z, 1.0, r, 0.0);
00331 
00332   elapsed_time = timer.ElapsedTime();
00333   total_flops = 2.0*global_dimension *((double) maxtrials);
00334   MFLOPs = total_flops/elapsed_time/1000000.0;
00335   results[ntimings][0] = total_flops;
00336   results[ntimings][1] = elapsed_time;
00337   results[ntimings++][2] = MFLOPs;
00338 
00339   if (makeheader)
00340     cout << "Metric_\t\t_Procs_\t__SpMV__\t_SpMM2__\t_SpMM4__\t_SpMM8__\t__NORM___\t__DOT___\t__AXPY__" << endl;
00341   if (comm.MyPID()==0) {
00342     cout.setf(std::ios::scientific);
00343     cout.precision(2);
00344     for (int j=0; j<3; j++) {
00345       bool doloop = false;
00346       if (j==0 && printflops) { cout << "FLOPS\t"; doloop = true;}
00347       else if (j==1 && printtime) {cout << "Time(s)\t"; doloop = true;}
00348       else if (j==2 && printmflops) {cout << "MFLOP/s\t"; doloop = true;}
00349       if (doloop) {
00350   cout << "\t" << comm.NumProc();
00351   for (int i=0; i<maxtimings; i++)
00352     cout << "\t" << results[i][j];
00353   cout << endl;
00354       }
00355     }
00356   }
00357 
00358 #ifdef EPETRA_MPI
00359   MPI_Finalize() ;
00360 #endif
00361 
00362 return ierr ;
00363 }
00364 
00365 // Constructs a 2D PDE finite difference matrix using the list of x and y offsets.
00366 //
00367 // nx      (In) - number of grid points in x direction
00368 // ny      (In) - number of grid points in y direction
00369 //   The total number of equations will be nx*ny ordered such that the x direction changes
00370 //   most rapidly:
00371 //      First equation is at point (0,0)
00372 //      Second at                  (1,0)
00373 //       ...
00374 //      nx equation at             (nx-1,0)
00375 //      nx+1st equation at         (0,1)
00376 
00377 // numPoints (In) - number of points in finite difference stencil
00378 // xoff    (In) - stencil offsets in x direction (of length numPoints)
00379 // yoff    (In) - stencil offsets in y direction (of length numPoints)
00380 //   A standard 5-point finite difference stencil would be described as:
00381 //     numPoints = 5
00382 //     xoff = [-1, 1, 0,  0, 0]
00383 //     yoff = [ 0, 0, 0, -1, 1]
00384 
00385 // nrhs - Number of rhs to generate. (First interface produces vectors, so nrhs is not needed
00386 
00387 // comm    (In) - an Epetra_Comm object describing the parallel machine (numProcs and my proc ID)
00388 // map    (Out) - Epetra_Map describing distribution of matrix and vectors/multivectors
00389 // A      (Out) - Epetra_CrsMatrix constructed for nx by ny grid using prescribed stencil
00390 //                Off-diagonal values are random between 0 and 1.  If diagonal is part of stencil,
00391 //                diagonal will be slightly diag dominant.
00392 // b      (Out) - Generated RHS.  Values satisfy b = A*xexact
00393 // bt     (Out) - Generated RHS.  Values satisfy b = A'*xexact
00394 // xexact (Out) - Generated exact solution to Ax = b and b' = A'xexact
00395 
00396 // Note: Caller of this function is responsible for deleting all output objects.
00397 
00398 void GenerateCrsProblem(int numNodesX, int numNodesY, int numProcsX, int numProcsY, int numPoints,
00399       int * xoff, int * yoff,
00400       const Epetra_Comm  &comm,
00401       Epetra_Map *& map,
00402       Epetra_CrsMatrix *& A,
00403       Epetra_Vector *& b,
00404       Epetra_Vector *& bt,
00405       Epetra_Vector *&xexact, bool StaticProfile, bool MakeLocalOnly) {
00406 
00407   Epetra_MultiVector * b1, * bt1, * xexact1;
00408   
00409   GenerateCrsProblem(numNodesX, numNodesY, numProcsX, numProcsY, numPoints,
00410          xoff, yoff, 1, comm,
00411          map, A, b1, bt1, xexact1, StaticProfile, MakeLocalOnly);
00412 
00413   b = dynamic_cast<Epetra_Vector *>(b1);
00414   bt = dynamic_cast<Epetra_Vector *>(bt1);
00415   xexact = dynamic_cast<Epetra_Vector *>(xexact1);
00416 
00417   return;
00418 }
00419 
00420 void GenerateCrsProblem(int numNodesX, int numNodesY, int numProcsX, int numProcsY, int numPoints,
00421       int * xoff, int * yoff, int nrhs,
00422       const Epetra_Comm  &comm,
00423       Epetra_Map *& map,
00424       Epetra_CrsMatrix *& A,
00425       Epetra_MultiVector *& b,
00426       Epetra_MultiVector *& bt,
00427       Epetra_MultiVector *&xexact, bool StaticProfile, bool MakeLocalOnly) {
00428 
00429   Epetra_Time timer(comm);
00430   // Determine my global IDs
00431   int * myGlobalElements;
00432   GenerateMyGlobalElements(numNodesX, numNodesY, numProcsX, numProcsY, comm.MyPID(), myGlobalElements);
00433 
00434   int numMyEquations = numNodesX*numNodesY;
00435 
00436   map = new Epetra_Map(-1, numMyEquations, myGlobalElements, 0, comm); // Create map with 2D block partitioning.
00437   delete [] myGlobalElements;
00438 
00439   int numGlobalEquations = map->NumGlobalElements();
00440 
00441   int profile = 0; if (StaticProfile) profile = numPoints;
00442 
00443 #ifdef EPETRA_HAVE_STATICPROFILE
00444 
00445   if (MakeLocalOnly)
00446     A = new Epetra_CrsMatrix(Copy, *map, *map, profile, StaticProfile); // Construct matrix with rowmap=colmap
00447   else
00448     A = new Epetra_CrsMatrix(Copy, *map, profile, StaticProfile); // Construct matrix
00449 
00450 #else
00451 
00452   if (MakeLocalOnly)
00453     A = new Epetra_CrsMatrix(Copy, *map, *map, profile); // Construct matrix with rowmap=colmap
00454   else
00455     A = new Epetra_CrsMatrix(Copy, *map, profile); // Construct matrix
00456 
00457 #endif
00458 
00459   int * indices = new int[numPoints];
00460   double * values = new double[numPoints];
00461 
00462   double dnumPoints = (double) numPoints;
00463   int nx = numNodesX*numProcsX;
00464 
00465   for (int i=0; i<numMyEquations; i++) {
00466 
00467     int rowID = map->GID(i);
00468     int numIndices = 0;
00469 
00470     for (int j=0; j<numPoints; j++) {
00471       int colID = rowID + xoff[j] + nx*yoff[j]; // Compute column ID based on stencil offsets
00472       if (colID>-1 && colID<numGlobalEquations) {
00473   indices[numIndices] = colID;
00474   double value = - ((double) rand())/ ((double) RAND_MAX);
00475   if (colID==rowID)
00476     values[numIndices++] = dnumPoints - value; // Make diagonal dominant
00477   else
00478     values[numIndices++] = value;
00479       }
00480     }
00481     //cout << "Building row " << rowID << endl;
00482     A->InsertGlobalValues(rowID, numIndices, values, indices);
00483   }
00484 
00485   delete [] indices;
00486   delete [] values;
00487   double insertTime = timer.ElapsedTime();
00488   timer.ResetStartTime();
00489   A->FillComplete();
00490   double fillCompleteTime = timer.ElapsedTime();
00491 
00492   if (nrhs<=1) {
00493     b = new Epetra_Vector(*map);
00494     bt = new Epetra_Vector(*map);
00495     xexact = new Epetra_Vector(*map);
00496   }
00497   else {
00498     b = new Epetra_MultiVector(*map, nrhs);
00499     bt = new Epetra_MultiVector(*map, nrhs);
00500     xexact = new Epetra_MultiVector(*map, nrhs);
00501   }
00502 
00503   xexact->Random(); // Fill xexact with random values
00504 
00505   A->Multiply(false, *xexact, *b);
00506   A->Multiply(true, *xexact, *bt);
00507 
00508   return;
00509 }
00510 
00511 
00512 
00513 
00514 void GenerateMyGlobalElements(int numNodesX, int numNodesY, int numProcsX, int numProcs,
00515             int myPID, int * & myGlobalElements) {
00516 
00517   myGlobalElements = new int[numNodesX*numNodesY];
00518   int myProcX = myPID%numProcsX;
00519   int myProcY = myPID/numProcsX;
00520   int curGID = myProcY*(numProcsX*numNodesX)*numNodesY+myProcX*numNodesX;
00521   for (int j=0; j<numNodesY; j++) {
00522     for (int i=0; i<numNodesX; i++) {
00523       myGlobalElements[j*numNodesX+i] = curGID+i;
00524     }
00525     curGID+=numNodesX*numProcsX;
00526   }
00527   //for (int i=0; i<numNodesX*numNodesY; i++) cout << "MYPID " << myPID <<" GID "<< myGlobalElements[i] << endl;
00528 
00529   return;
00530 }
00531 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines