Epetra Package Browser (Single Doxygen Collection) Development
test/CrsMatrix/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 
00043 #include "Epetra_Map.h"
00044 #include "Epetra_Time.h"
00045 #include "Epetra_CrsMatrix.h"
00046 #include "Epetra_Vector.h"
00047 #include "Epetra_Flops.h"
00048 #ifdef EPETRA_MPI
00049 #include "Epetra_MpiComm.h"
00050 #include "mpi.h"
00051 #else
00052 #include "Epetra_SerialComm.h"
00053 #endif
00054 #include "../epetra_test_err.h"
00055 #include "Epetra_Version.h"
00056 
00057 // prototypes
00058 
00059 int check(Epetra_CrsMatrix& A, int NumMyRows1, int NumGlobalRows1, int NumMyNonzeros1,
00060     int NumGlobalNonzeros1, int * MyGlobalElements, bool verbose);
00061 
00062 int power_method(bool TransA, Epetra_CrsMatrix& A,
00063      Epetra_Vector& q,
00064      Epetra_Vector& z,
00065      Epetra_Vector& resid,
00066      double * lambda, int niters, double tolerance,
00067      bool verbose);
00068 
00069 int check_graph_sharing(Epetra_Comm& Comm);
00070 
00071 int main(int argc, char *argv[])
00072 {
00073   int ierr = 0, forierr = 0;
00074   bool debug = false;
00075 
00076 #ifdef EPETRA_MPI
00077 
00078   // Initialize MPI
00079 
00080   MPI_Init(&argc,&argv);
00081   int rank; // My process ID
00082 
00083   MPI_Comm_rank(MPI_COMM_WORLD, &rank);
00084   Epetra_MpiComm Comm( MPI_COMM_WORLD );
00085 
00086 #else
00087 
00088   int rank = 0;
00089   Epetra_SerialComm Comm;
00090 
00091 #endif
00092 
00093   bool verbose = false;
00094 
00095   // Check if we should print results to standard out
00096   if (argc>1) if (argv[1][0]=='-' && argv[1][1]=='v') verbose = true;
00097 
00098   int verbose_int = verbose ? 1 : 0;
00099   Comm.Broadcast(&verbose_int, 1, 0);
00100   verbose = verbose_int==1 ? true : false;
00101 
00102 
00103   //  char tmp;
00104   //  if (rank==0) cout << "Press any key to continue..."<< std::endl;
00105   //  if (rank==0) cin >> tmp;
00106   //  Comm.Barrier();
00107 
00108   Comm.SetTracebackMode(0); // This should shut down any error traceback reporting
00109   int MyPID = Comm.MyPID();
00110   int NumProc = Comm.NumProc();
00111 
00112   if(verbose && MyPID==0)
00113     cout << Epetra_Version() << std::endl << std::endl;
00114 
00115   if (verbose) cout << "Processor "<<MyPID<<" of "<< NumProc
00116         << " is alive."<<endl;
00117 
00118   bool verbose1 = verbose;
00119 
00120   // Redefine verbose to only print on PE 0
00121   if(verbose && rank!=0)
00122     verbose = false;
00123 
00124   int NumMyEquations = 10000;
00125   int NumGlobalEquations = (NumMyEquations * NumProc) + EPETRA_MIN(NumProc,3);
00126   if(MyPID < 3)
00127     NumMyEquations++;
00128 
00129   // Construct a Map that puts approximately the same Number of equations on each processor
00130 
00131   Epetra_Map Map(NumGlobalEquations, NumMyEquations, 0, Comm);
00132 
00133   // Get update list and number of local equations from newly created Map
00134   int* MyGlobalElements = new int[Map.NumMyElements()];
00135   Map.MyGlobalElements(MyGlobalElements);
00136 
00137   // Create an integer vector NumNz that is used to build the Petra Matrix.
00138   // NumNz[i] is the Number of OFF-DIAGONAL term for the ith global equation on this processor
00139 
00140   int* NumNz = new int[NumMyEquations];
00141 
00142   // We are building a tridiagonal matrix where each row has (-1 2 -1)
00143   // So we need 2 off-diagonal terms (except for the first and last equation)
00144 
00145   for (int i = 0; i < NumMyEquations; i++)
00146     if((MyGlobalElements[i] == 0) || (MyGlobalElements[i] == NumGlobalEquations - 1))
00147       NumNz[i] = 1;
00148     else
00149       NumNz[i] = 2;
00150 
00151   // Create a Epetra_Matrix
00152 
00153   Epetra_CrsMatrix A(Copy, Map, NumNz);
00154   EPETRA_TEST_ERR(A.IndicesAreGlobal(),ierr);
00155   EPETRA_TEST_ERR(A.IndicesAreLocal(),ierr);
00156 
00157   // Add  rows one-at-a-time
00158   // Need some vectors to help
00159   // Off diagonal Values will always be -1
00160 
00161 
00162   double* Values = new double[2];
00163   Values[0] = -1.0;
00164   Values[1] = -1.0;
00165   int* Indices = new int[2];
00166   double two = 2.0;
00167   int NumEntries;
00168 
00169   forierr = 0;
00170   for (int i = 0; i < NumMyEquations; i++) {
00171     if(MyGlobalElements[i] == 0) {
00172       Indices[0] = 1;
00173       NumEntries = 1;
00174     }
00175     else if (MyGlobalElements[i] == NumGlobalEquations-1) {
00176       Indices[0] = NumGlobalEquations-2;
00177       NumEntries = 1;
00178     }
00179     else {
00180       Indices[0] = MyGlobalElements[i]-1;
00181       Indices[1] = MyGlobalElements[i]+1;
00182       NumEntries = 2;
00183     }
00184     forierr += !(A.InsertGlobalValues(MyGlobalElements[i], NumEntries, Values, Indices)==0);
00185     forierr += !(A.InsertGlobalValues(MyGlobalElements[i], 1, &two, MyGlobalElements+i)>0); // Put in the diagonal entry
00186   }
00187   EPETRA_TEST_ERR(forierr,ierr);
00188 
00189   int * indexOffsetTmp;
00190   int * indicesTmp;
00191   double * valuesTmp;
00192   // Finish up
00193   EPETRA_TEST_ERR(!(A.IndicesAreGlobal()),ierr);
00194   EPETRA_TEST_ERR(!(A.ExtractCrsDataPointers(indexOffsetTmp, indicesTmp, valuesTmp)==-1),ierr);  // Should fail
00195   EPETRA_TEST_ERR(!(A.FillComplete(false)==0),ierr);
00196   EPETRA_TEST_ERR(!(A.ExtractCrsDataPointers(indexOffsetTmp, indicesTmp, valuesTmp)==-1),ierr);  // Should fail
00197   EPETRA_TEST_ERR(!(A.IndicesAreLocal()),ierr);
00198   EPETRA_TEST_ERR(A.StorageOptimized(),ierr);
00199   A.OptimizeStorage();
00200   EPETRA_TEST_ERR(!(A.StorageOptimized()),ierr);
00201   EPETRA_TEST_ERR(!(A.ExtractCrsDataPointers(indexOffsetTmp, indicesTmp, valuesTmp)==0),ierr);  // Should succeed
00202   const Epetra_CrsGraph & GofA = A.Graph();
00203   EPETRA_TEST_ERR((indicesTmp!=GofA[0] || valuesTmp!=A[0]),ierr); // Extra check to see if operator[] is consistent
00204   EPETRA_TEST_ERR(A.UpperTriangular(),ierr);
00205   EPETRA_TEST_ERR(A.LowerTriangular(),ierr);
00206   
00207   int NumMyNonzeros = 3 * NumMyEquations;
00208   if(A.LRID(0) >= 0)
00209     NumMyNonzeros--; // If I own first global row, then there is one less nonzero
00210   if(A.LRID(NumGlobalEquations-1) >= 0)
00211     NumMyNonzeros--; // If I own last global row, then there is one less nonzero
00212   EPETRA_TEST_ERR(check(A, NumMyEquations, NumGlobalEquations, NumMyNonzeros, 3*NumGlobalEquations-2,
00213          MyGlobalElements, verbose),ierr);
00214   forierr = 0;
00215   for (int i = 0; i < NumMyEquations; i++)
00216     forierr += !(A.NumGlobalEntries(MyGlobalElements[i])==NumNz[i]+1);
00217   EPETRA_TEST_ERR(forierr,ierr);
00218   forierr = 0;
00219   for (int i = 0; i < NumMyEquations; i++)
00220     forierr += !(A.NumMyEntries(i)==NumNz[i]+1);
00221   EPETRA_TEST_ERR(forierr,ierr);
00222 
00223   if (verbose) cout << "\n\nNumEntries function check OK" << std::endl<< std::endl;
00224 
00225   EPETRA_TEST_ERR(check_graph_sharing(Comm),ierr);
00226 
00227   // Create vectors for Power method
00228 
00229   Epetra_Vector q(Map);
00230   Epetra_Vector z(Map);
00231   Epetra_Vector resid(Map);
00232 
00233   // variable needed for iteration
00234   double lambda = 0.0;
00235   // int niters = 10000;
00236   int niters = 200;
00237   double tolerance = 1.0e-1;
00238 
00240   
00241   // Iterate
00242 
00243   Epetra_Flops flopcounter;
00244   A.SetFlopCounter(flopcounter);
00245   q.SetFlopCounter(A);
00246   z.SetFlopCounter(A);
00247   resid.SetFlopCounter(A);
00248   
00249 
00250   Epetra_Time timer(Comm);
00251   EPETRA_TEST_ERR(power_method(false, A, q, z, resid, &lambda, niters, tolerance, verbose),ierr);
00252   double elapsed_time = timer.ElapsedTime();
00253   double total_flops = A.Flops() + q.Flops() + z.Flops() + resid.Flops();
00254   double MFLOPs = total_flops/elapsed_time/1000000.0;
00255 
00256   if (verbose) cout << "\n\nTotal MFLOPs for first solve = " << MFLOPs << std::endl<< std::endl;
00257 
00259   
00260   // Solve transpose problem
00261 
00262   if (verbose) cout << "\n\nUsing transpose of matrix and solving again (should give same result).\n\n"
00263         << std::endl;
00264   // Iterate
00265   lambda = 0.0;
00266   flopcounter.ResetFlops();
00267   timer.ResetStartTime();
00268   EPETRA_TEST_ERR(power_method(true, A, q, z, resid, &lambda, niters, tolerance, verbose),ierr);
00269   elapsed_time = timer.ElapsedTime();
00270   total_flops = A.Flops() + q.Flops() + z.Flops() + resid.Flops();
00271   MFLOPs = total_flops/elapsed_time/1000000.0;
00272 
00273   if (verbose) cout << "\n\nTotal MFLOPs for transpose solve = " << MFLOPs << std::endl<< endl;
00274 
00276 
00277   // Increase diagonal dominance
00278 
00279   if (verbose) cout << "\n\nIncreasing the magnitude of first diagonal term and solving again\n\n"
00280         << endl;
00281 
00282 
00283   if (A.MyGlobalRow(0)) {
00284     int numvals = A.NumGlobalEntries(0);
00285     double * Rowvals = new double [numvals];
00286     int    * Rowinds = new int    [numvals];
00287     A.ExtractGlobalRowCopy(0, numvals, numvals, Rowvals, Rowinds); // Get A[0,0]
00288 
00289     for (int i=0; i<numvals; i++) if (Rowinds[i] == 0) Rowvals[i] *= 10.0;
00290 
00291     A.ReplaceGlobalValues(0, numvals, Rowvals, Rowinds);
00292     delete [] Rowvals;
00293     delete [] Rowinds;
00294   }
00295   // Iterate (again)
00296   lambda = 0.0;
00297   flopcounter.ResetFlops();
00298   timer.ResetStartTime();
00299   EPETRA_TEST_ERR(power_method(false, A, q, z, resid, &lambda, niters, tolerance, verbose),ierr);
00300   elapsed_time = timer.ElapsedTime();
00301   total_flops = A.Flops() + q.Flops() + z.Flops() + resid.Flops();
00302   MFLOPs = total_flops/elapsed_time/1000000.0;
00303 
00304   if (verbose) cout << "\n\nTotal MFLOPs for second solve = " << MFLOPs << endl<< endl;
00305 
00307 
00308   // Solve transpose problem
00309 
00310   if (verbose) cout << "\n\nUsing transpose of matrix and solving again (should give same result).\n\n"
00311         << endl;
00312 
00313   // Iterate (again)
00314   lambda = 0.0;
00315   flopcounter.ResetFlops();
00316   timer.ResetStartTime();
00317   EPETRA_TEST_ERR(power_method(true, A, q, z, resid, &lambda, niters, tolerance, verbose),ierr);
00318   elapsed_time = timer.ElapsedTime();
00319   total_flops = A.Flops() + q.Flops() + z.Flops() + resid.Flops();
00320   MFLOPs = total_flops/elapsed_time/1000000.0;
00321 
00322 
00323   if (verbose) cout << "\n\nTotal MFLOPs for tranpose of second solve = " << MFLOPs << endl<< endl;
00324 
00325   if (verbose) cout << "\n\n*****Testing constant entry constructor" << endl<< endl;
00326 
00327   Epetra_CrsMatrix AA(Copy, Map, 5);
00328 
00329   if (debug) Comm.Barrier();
00330 
00331   double dble_one = 1.0;
00332   for (int i=0; i< NumMyEquations; i++) AA.InsertGlobalValues(MyGlobalElements[i], 1, &dble_one, MyGlobalElements+i);
00333 
00334   // Note:  All processors will call the following Insert routines, but only the processor
00335   //        that owns it will actually do anything
00336 
00337   int One = 1;
00338   if (AA.MyGlobalRow(0)) {
00339     EPETRA_TEST_ERR(!(AA.InsertGlobalValues(0, 0, &dble_one, &One)==0),ierr);
00340   }
00341   else EPETRA_TEST_ERR(!(AA.InsertGlobalValues(0, 1, &dble_one, &One)==-1),ierr);
00342   EPETRA_TEST_ERR(!(AA.FillComplete(false)==0),ierr);
00343   EPETRA_TEST_ERR(AA.StorageOptimized(),ierr);
00344   EPETRA_TEST_ERR(!(AA.UpperTriangular()),ierr);
00345   EPETRA_TEST_ERR(!(AA.LowerTriangular()),ierr);
00346 
00347   if (debug) Comm.Barrier();
00348   EPETRA_TEST_ERR(check(AA, NumMyEquations, NumGlobalEquations, NumMyEquations, NumGlobalEquations,
00349          MyGlobalElements, verbose),ierr);
00350 
00351   if (debug) Comm.Barrier();
00352 
00353   forierr = 0;
00354   for (int i=0; i<NumMyEquations; i++) forierr += !(AA.NumGlobalEntries(MyGlobalElements[i])==1);
00355   EPETRA_TEST_ERR(forierr,ierr);
00356 
00357   if (verbose) cout << "\n\nNumEntries function check OK" << endl<< endl;
00358 
00359   if (debug) Comm.Barrier();
00360 
00361   if (verbose) cout << "\n\n*****Testing copy constructor" << endl<< endl;
00362 
00363   Epetra_CrsMatrix B(AA);
00364   EPETRA_TEST_ERR(check(B, NumMyEquations, NumGlobalEquations, NumMyEquations, NumGlobalEquations,
00365          MyGlobalElements, verbose),ierr);
00366 
00367   forierr = 0;
00368   for (int i=0; i<NumMyEquations; i++) forierr += !(B.NumGlobalEntries(MyGlobalElements[i])==1);
00369   EPETRA_TEST_ERR(forierr,ierr);
00370 
00371   if (verbose) cout << "\n\nNumEntries function check OK" << endl<< endl;
00372 
00373   if (debug) Comm.Barrier();
00374 
00375   if (verbose) cout << "\n\n*****Testing local view constructor" << endl<< endl;
00376 
00377   Epetra_CrsMatrix BV(View, AA.RowMap(), AA.ColMap(), 0);
00378 
00379   forierr = 0;
00380   int* Inds;
00381   double* Vals;
00382   for (int i = 0; i < NumMyEquations; i++) {
00383     forierr += !(AA.ExtractMyRowView(i, NumEntries, Vals, Inds)==0);
00384     forierr += !(BV.InsertMyValues(i, NumEntries, Vals, Inds)==0);
00385   }
00386   BV.FillComplete(false);
00387   EPETRA_TEST_ERR(check(BV, NumMyEquations, NumGlobalEquations, NumMyEquations, NumGlobalEquations,
00388                         MyGlobalElements, verbose),ierr);
00389 
00390   forierr = 0;
00391   for (int i=0; i<NumMyEquations; i++) forierr += !(BV.NumGlobalEntries(MyGlobalElements[i])==1);
00392   EPETRA_TEST_ERR(forierr,ierr);
00393 
00394   if (verbose) cout << "\n\nNumEntries function check OK" << endl<< endl;
00395 
00396   if (debug) Comm.Barrier();
00397   if (verbose) cout << "\n\n*****Testing post construction modifications" << endl<< endl;
00398 
00399   EPETRA_TEST_ERR(!(B.InsertGlobalValues(0, 1, &dble_one, &One)==-2),ierr);
00400 
00401 
00402   // Release all objects
00403   delete [] NumNz;
00404   delete [] Values;
00405   delete [] Indices;
00406   delete [] MyGlobalElements;
00407       
00408 
00409   if (verbose1) {
00410     // Test ostream << operator (if verbose1)
00411     // Construct a Map that puts 2 equations on each PE
00412 
00413     int NumMyElements1 = 2;
00414     int NumMyEquations1 = NumMyElements1;
00415     int NumGlobalEquations1 = NumMyEquations1*NumProc;
00416 
00417     Epetra_Map Map1(-1, NumMyElements1, 0, Comm);
00418 
00419     // Get update list and number of local equations from newly created Map
00420     int * MyGlobalElements1 = new int[Map1.NumMyElements()];
00421     Map1.MyGlobalElements(MyGlobalElements1);
00422 
00423     // Create an integer vector NumNz that is used to build the Petra Matrix.
00424     // NumNz[i] is the Number of OFF-DIAGONAL term for the ith global equation on this processor
00425 
00426     int * NumNz1 = new int[NumMyEquations1];
00427 
00428     // We are building a tridiagonal matrix where each row has (-1 2 -1)
00429     // So we need 2 off-diagonal terms (except for the first and last equation)
00430 
00431     for (int i=0; i<NumMyEquations1; i++)
00432       if (MyGlobalElements1[i]==0 || MyGlobalElements1[i] == NumGlobalEquations1-1)
00433   NumNz1[i] = 1;
00434       else
00435   NumNz1[i] = 2;
00436 
00437     // Create a Epetra_Matrix
00438 
00439     Epetra_CrsMatrix A1(Copy, Map1, NumNz1);
00440 
00441     // Add  rows one-at-a-time
00442     // Need some vectors to help
00443     // Off diagonal Values will always be -1
00444 
00445 
00446     double *Values1 = new double[2];
00447     Values1[0] = -1.0; Values1[1] = -1.0;
00448     int *Indices1 = new int[2];
00449     double two1 = 2.0;
00450     int NumEntries1;
00451 
00452     forierr = 0;
00453     for (int i=0; i<NumMyEquations1; i++)
00454       {
00455   if (MyGlobalElements1[i]==0)
00456     {
00457       Indices1[0] = 1;
00458       NumEntries1 = 1;
00459     }
00460   else if (MyGlobalElements1[i] == NumGlobalEquations1-1)
00461     {
00462       Indices1[0] = NumGlobalEquations1-2;
00463       NumEntries1 = 1;
00464     }
00465   else
00466     {
00467       Indices1[0] = MyGlobalElements1[i]-1;
00468       Indices1[1] = MyGlobalElements1[i]+1;
00469       NumEntries1 = 2;
00470     }
00471   forierr += !(A1.InsertGlobalValues(MyGlobalElements1[i], NumEntries1, Values1, Indices1)==0);
00472   forierr += !(A1.InsertGlobalValues(MyGlobalElements1[i], 1, &two1, MyGlobalElements1+i)>0); // Put in the diagonal entry
00473       }
00474     EPETRA_TEST_ERR(forierr,ierr);
00475     delete [] Indices1;
00476     delete [] Values1;
00477 
00478     // Finish up
00479     EPETRA_TEST_ERR(!(A1.FillComplete(false)==0),ierr);
00480 
00481     // Test diagonal extraction function
00482 
00483     Epetra_Vector checkDiag(Map1);
00484     EPETRA_TEST_ERR(!(A1.ExtractDiagonalCopy(checkDiag)==0),ierr);
00485 
00486     forierr = 0;
00487     for (int i=0; i<NumMyEquations1; i++) forierr += !(checkDiag[i]==two1);
00488     EPETRA_TEST_ERR(forierr,ierr);
00489 
00490     // Test diagonal replacement method
00491 
00492     forierr = 0;
00493     for (int i=0; i<NumMyEquations1; i++) checkDiag[i]=two1*two1;
00494     EPETRA_TEST_ERR(forierr,ierr);
00495 
00496     EPETRA_TEST_ERR(!(A1.ReplaceDiagonalValues(checkDiag)==0),ierr);
00497 
00498     Epetra_Vector checkDiag1(Map1);
00499     EPETRA_TEST_ERR(!(A1.ExtractDiagonalCopy(checkDiag1)==0),ierr);
00500 
00501     forierr = 0;
00502     for (int i=0; i<NumMyEquations1; i++) forierr += !(checkDiag[i]==checkDiag1[i]);
00503     EPETRA_TEST_ERR(forierr,ierr);
00504 
00505     if (verbose) cout << "\n\nDiagonal extraction and replacement OK.\n\n" << endl;
00506 
00507     double orignorm = A1.NormOne();
00508     EPETRA_TEST_ERR(!(A1.Scale(4.0)==0),ierr);
00509     EPETRA_TEST_ERR(!(A1.NormOne()!=orignorm),ierr);
00510 
00511     if (verbose) cout << "\n\nMatrix scale OK.\n\n" << endl;
00512 
00513     if (verbose) cout << "\n\nPrint out tridiagonal matrix, each part on each processor.\n\n" << endl;
00514     cout << A1 << endl;
00515 
00516 
00517   // Release all objects
00518   delete [] NumNz1;
00519   delete [] MyGlobalElements1;
00520 
00521   }
00522 
00523   if (verbose) cout << "\n\n*****Testing LeftScale and RightScale" << endl << endl;
00524 
00525   int NumMyElements2 = 7;
00526   int NumMyRows2 = 1;//This value should not be changed without editing the
00527     // code below.
00528   Epetra_Map RowMap(-1,NumMyRows2,0,Comm);
00529   Epetra_Map ColMap(NumMyElements2,NumMyElements2,0,Comm);
00530   // The DomainMap needs to be different from the ColMap for the test to
00531   // be meaningful.
00532   Epetra_Map DomainMap(NumMyElements2,0,Comm);
00533   int NumMyRangeElements2 = 0;
00534   // We need to distribute the elements differently for the range map also.
00535   if (MyPID % 2 == 0)
00536     NumMyRangeElements2 = NumMyRows2*2; //put elements on even number procs
00537   if (NumProc % 2 == 1 && MyPID == NumProc-1)
00538     NumMyRangeElements2 = NumMyRows2; //If number of procs is odd, put
00539       // the last NumMyElements2 elements on the last proc
00540   Epetra_Map RangeMap(-1,NumMyRangeElements2,0,Comm);
00541   Epetra_CrsMatrix A2(Copy,RowMap,ColMap,NumMyElements2);
00542   double * Values2 = new double[NumMyElements2];
00543   int * Indices2 = new int[NumMyElements2];
00544 
00545   for (int i=0; i<NumMyElements2; i++) {
00546     Values2[i] = i+MyPID;
00547     Indices2[i]=i;
00548   }
00549 
00550   A2.InsertMyValues(0,NumMyElements2,Values2,Indices2);
00551   A2.FillComplete(DomainMap,RangeMap,false);
00552   Epetra_CrsMatrix A2copy(A2);
00553 
00554   double * RowLeftScaleValues = new double[NumMyRows2];
00555   double * ColRightScaleValues = new double[NumMyElements2];
00556   int RowLoopLength = RowMap.MaxMyGID()-RowMap.MinMyGID()+1;
00557   for (int i=0; i<RowLoopLength; i++)
00558     RowLeftScaleValues[i] = (i + RowMap.MinMyGID() ) % 2 + 1;
00559   // For the column map, all procs own all elements
00560   for (int  i=0; i<NumMyElements2;i++)
00561     ColRightScaleValues[i] = i % 2 + 1;
00562 
00563   int RangeLoopLength = RangeMap.MaxMyGID()-RangeMap.MinMyGID()+1;
00564   double * RangeLeftScaleValues = new double[RangeLoopLength];
00565   int DomainLoopLength = DomainMap.MaxMyGID()-DomainMap.MinMyGID()+1;
00566    double * DomainRightScaleValues = new double[DomainLoopLength];
00567   for (int i=0; i<RangeLoopLength; i++)
00568     RangeLeftScaleValues[i] = 1.0/((i + RangeMap.MinMyGID() ) % 2 + 1);
00569   for (int  i=0; i<DomainLoopLength;i++)
00570     DomainRightScaleValues[i] = 1.0/((i + DomainMap.MinMyGID() ) % 2 + 1);
00571 
00572   Epetra_Vector xRow(View,RowMap,RowLeftScaleValues);
00573   Epetra_Vector xCol(View,ColMap,ColRightScaleValues);
00574   Epetra_Vector xRange(View,RangeMap,RangeLeftScaleValues);
00575   Epetra_Vector xDomain(View,DomainMap,DomainRightScaleValues);
00576 
00577   double A2infNorm = A2.NormInf();
00578   double A2oneNorm = A2.NormOne();
00579 
00580   if (verbose1) cout << A2;
00581   EPETRA_TEST_ERR(A2.LeftScale(xRow),ierr);
00582   double A2infNorm1 = A2.NormInf();
00583   double A2oneNorm1 = A2.NormOne();
00584   bool ScalingBroke = false;
00585   if (A2infNorm1>2*A2infNorm||A2infNorm1<A2infNorm) {
00586     EPETRA_TEST_ERR(-31,ierr);
00587     ScalingBroke = true;
00588   }
00589   if (A2oneNorm1>2*A2oneNorm||A2oneNorm1<A2oneNorm) {
00590 
00591     EPETRA_TEST_ERR(-32,ierr);
00592     ScalingBroke = true;
00593   }
00594   if (verbose1) cout << A2;
00595   EPETRA_TEST_ERR(A2.RightScale(xCol),ierr);
00596   double A2infNorm2 = A2.NormInf();
00597   double A2oneNorm2 = A2.NormOne();
00598   if (A2infNorm2>=2*A2infNorm1||A2infNorm2<=A2infNorm1) {
00599     EPETRA_TEST_ERR(-33,ierr);
00600     ScalingBroke = true;
00601   }
00602   if (A2oneNorm2>2*A2oneNorm1||A2oneNorm2<=A2oneNorm1) {
00603     EPETRA_TEST_ERR(-34,ierr);
00604     ScalingBroke = true;
00605   }
00606   if (verbose1) cout << A2;
00607   EPETRA_TEST_ERR(A2.RightScale(xDomain),ierr);
00608   double A2infNorm3 = A2.NormInf();
00609   double A2oneNorm3 = A2.NormOne();
00610   // The last two scaling ops cancel each other out
00611   if (A2infNorm3!=A2infNorm1) {
00612     EPETRA_TEST_ERR(-35,ierr)
00613     ScalingBroke = true;
00614   }
00615   if (A2oneNorm3!=A2oneNorm1) {
00616     EPETRA_TEST_ERR(-36,ierr)
00617     ScalingBroke = true;
00618   }
00619   if (verbose1) cout << A2;
00620   EPETRA_TEST_ERR(A2.LeftScale(xRange),ierr);
00621   double A2infNorm4 = A2.NormInf();
00622   double A2oneNorm4 = A2.NormOne();
00623   // The 4 scaling ops all cancel out
00624   if (A2infNorm4!=A2infNorm) {
00625     EPETRA_TEST_ERR(-37,ierr)
00626     ScalingBroke = true;
00627   }
00628   if (A2oneNorm4!=A2oneNorm) {
00629     EPETRA_TEST_ERR(-38,ierr)
00630     ScalingBroke = true;
00631   }
00632 
00633   //
00634   //  Now try changing the values underneath and make sure that
00635   //  telling one process about the change causes NormInf() and
00636   //  NormOne() to recompute the norm on all processes.
00637   //
00638 
00639   double *values;
00640   int num_my_rows = A2.NumMyRows() ;
00641   int num_entries;
00642 
00643   for ( int  i=0 ; i< num_my_rows; i++ ) {
00644     EPETRA_TEST_ERR( A2.ExtractMyRowView( i, num_entries, values ), ierr );
00645     for ( int j = 0 ; j <num_entries; j++ ) {
00646       values[j] *= 2.0;
00647     }
00648   }
00649 
00650 
00651   if ( MyPID == 0 )
00652     A2.SumIntoGlobalValues( 0, 0, 0, 0 ) ;
00653 
00654   double A2infNorm5 = A2.NormInf();
00655   double A2oneNorm5 = A2.NormOne();
00656 
00657   if (A2infNorm5!=2.0 * A2infNorm4) {
00658     EPETRA_TEST_ERR(-39,ierr)
00659     ScalingBroke = true;
00660   }
00661   if (A2oneNorm5!= 2.0 * A2oneNorm4) {
00662     EPETRA_TEST_ERR(-40,ierr)
00663     ScalingBroke = true;
00664   }
00665 
00666   //
00667   //  Restore the values underneath
00668   //
00669   for ( int  i=0 ; i< num_my_rows; i++ ) {
00670     EPETRA_TEST_ERR( A2.ExtractMyRowView( i, num_entries, values ), ierr );
00671     for ( int j = 0 ; j <num_entries; j++ ) {
00672       values[j] /= 2.0;
00673     }
00674   }
00675 
00676   if (verbose1) cout << A2;
00677 
00678   if (ScalingBroke) {
00679     if (verbose) cout << endl << "LeftScale and RightScale tests FAILED" << endl << endl;
00680   }
00681   else {
00682     if (verbose) cout << endl << "LeftScale and RightScale tests PASSED" << endl << endl;
00683   }
00684 
00685   Comm.Barrier();
00686 
00687   if (verbose) cout << "\n\n*****Testing InvRowMaxs and InvColMaxs" << endl << endl;
00688 
00689   if (verbose1) cout << A2 << endl;
00690   EPETRA_TEST_ERR(A2.InvRowMaxs(xRow),ierr);
00691   EPETRA_TEST_ERR(A2.InvRowMaxs(xRange),ierr);
00692   if (verbose1) cout << xRow << endl << xRange << endl;
00693 
00694   if (verbose) cout << "\n\n*****Testing InvRowSums and InvColSums" << endl << endl;
00695   bool InvSumsBroke = false;
00696 // Works!
00697   EPETRA_TEST_ERR(A2.InvRowSums(xRow),ierr);
00698   if (verbose1) cout << xRow;
00699   EPETRA_TEST_ERR(A2.LeftScale(xRow),ierr);
00700   float A2infNormFloat = A2.NormInf();
00701   if (verbose1) cout << A2 << endl;
00702   if (fabs(1.0-A2infNormFloat) > 1.e-5) {
00703     EPETRA_TEST_ERR(-41,ierr);
00704     InvSumsBroke = true;
00705   }
00706 
00707   // Works
00708   int expectedcode = 1;
00709   if (Comm.NumProc()>1) expectedcode = 0;
00710   EPETRA_TEST_ERR(!(A2.InvColSums(xDomain)==expectedcode),ierr); // This matrix has a single row, the first column has a zero, so a warning is issued.
00711   if (verbose1) cout << xDomain << endl;
00712   EPETRA_TEST_ERR(A2.RightScale(xDomain),ierr);
00713   float A2oneNormFloat2 = A2.NormOne();
00714   if (verbose1) cout << A2;
00715   if (fabs(1.0-A2oneNormFloat2)>1.e-5) {
00716     EPETRA_TEST_ERR(-42,ierr)
00717     InvSumsBroke = true;
00718   }
00719 
00720 // Works!
00721   EPETRA_TEST_ERR(A2.InvRowSums(xRange),ierr);
00722 
00723   if (verbose1) cout << xRange;
00724   EPETRA_TEST_ERR(A2.LeftScale(xRange),ierr);
00725   float A2infNormFloat2 = A2.NormInf(); // We use a float so that rounding error
00726   // will not prevent the sum from being 1.0.
00727   if (verbose1) cout << A2;
00728   if (fabs(1.0-A2infNormFloat2)>1.e-5) {
00729     cout << "InfNorm should be = 1, but InfNorm = " << A2infNormFloat2 << endl;
00730     EPETRA_TEST_ERR(-43,ierr);
00731     InvSumsBroke = true;
00732   }
00733 
00734   // Doesn't work - may not need this test because column ownership is not unique
00735   /*  EPETRA_TEST_ERR(A2.InvColSums(xCol),ierr);
00736 cout << xCol;
00737   EPETRA_TEST_ERR(A2.RightScale(xCol),ierr);
00738   float A2oneNormFloat = A2.NormOne();
00739 cout << A2;
00740   if (fabs(1.0-A2oneNormFloat)>1.e-5) {
00741     EPETRA_TEST_ERR(-44,ierr);
00742     InvSumsBroke = true;
00743   }
00744   */
00745   delete [] ColRightScaleValues;
00746   delete [] DomainRightScaleValues;
00747   if (verbose) cout << "Begin partial sum testing." << endl;
00748   // Test with a matrix that has partial sums for a subset of the rows
00749   // on multiple processors. (Except for the serial case, of course.)
00750   int NumMyRows3 = 2; // Changing this requires further changes below
00751   int * myGlobalElements = new int[NumMyRows3];
00752   for (int i=0; i<NumMyRows3; i++) myGlobalElements[i] = MyPID+i;
00753   Epetra_Map RowMap3(NumProc*2, NumMyRows3, myGlobalElements, 0, Comm);
00754   int NumMyElements3 = 5;
00755   Epetra_CrsMatrix A3(Copy, RowMap3, NumMyElements3);
00756   double * Values3 = new double[NumMyElements3];
00757   int * Indices3 = new int[NumMyElements3];
00758   for (int i=0; i < NumMyElements3; i++) {
00759     Values3[i] = (int) (MyPID + (i+1));
00760     Indices3[i]=i;
00761   }
00762   for (int i=0; i<NumMyRows3; i++) {
00763     A3.InsertGlobalValues(myGlobalElements[i],NumMyElements3,Values3,Indices3);
00764   }
00765   Epetra_Map RangeMap3(NumProc+1, 0, Comm);
00766   Epetra_Map DomainMap3(NumMyElements3, 0, Comm);
00767   EPETRA_TEST_ERR(A3.FillComplete(DomainMap3, RangeMap3,false),ierr);
00768   if (verbose1) cout << A3;
00769   Epetra_Vector xRange3(RangeMap3,false);
00770   Epetra_Vector xDomain3(DomainMap3,false);
00771 
00772   EPETRA_TEST_ERR(A3.InvRowSums(xRange3),ierr);
00773 
00774   if (verbose1) cout << xRange3;
00775   EPETRA_TEST_ERR(A3.LeftScale(xRange3),ierr);
00776   float A3infNormFloat = A3.NormInf();
00777   if (verbose1) cout << A3;
00778   if (1.0!=A3infNormFloat) {
00779     cout << "InfNorm should be = 1, but InfNorm = " << A3infNormFloat <<endl;
00780     EPETRA_TEST_ERR(-61,ierr);
00781     InvSumsBroke = true;
00782   }
00783   // we want to take the transpose of our matrix and fill in different values.
00784   int NumMyColumns3 = NumMyRows3;
00785   Epetra_Map ColMap3cm(RowMap3);
00786   Epetra_Map RowMap3cm(A3.ColMap());
00787 
00788   Epetra_CrsMatrix A3cm(Copy,RowMap3cm,ColMap3cm,NumProc+1);
00789   double *Values3cm = new double[NumMyColumns3];
00790   int * Indices3cm = new int[NumMyColumns3];
00791   for (int i=0; i<NumMyColumns3; i++) {
00792     Values3cm[i] = MyPID + i + 1;
00793     Indices3cm[i]= i + MyPID;
00794   }
00795   for (int ii=0; ii<NumMyElements3; ii++) {
00796     A3cm.InsertGlobalValues(ii, NumMyColumns3, Values3cm, Indices3cm);
00797   }
00798 
00799   // The DomainMap and the RangeMap from the last test will work fine for
00800   // the RangeMap and DomainMap, respectively, but I will make copies to
00801   // avaoid confusion when passing what looks like a DomainMap where we
00802   // need a RangeMap and vice vera.
00803   Epetra_Map RangeMap3cm(DomainMap3);
00804   Epetra_Map DomainMap3cm(RangeMap3);
00805   EPETRA_TEST_ERR(A3cm.FillComplete(DomainMap3cm,RangeMap3cm),ierr);
00806   if (verbose1) cout << A3cm << endl;
00807 
00808   // Again, we can copy objects from the last example.
00809   //Epetra_Vector xRange3cm(xDomain3); //Don't use at this time
00810   Epetra_Vector xDomain3cm(DomainMap3cm,false);
00811 
00812   EPETRA_TEST_ERR(A3cm.InvColSums(xDomain3cm),ierr);
00813 
00814   if (verbose1) cout << xDomain3cm << endl;
00815 
00816   EPETRA_TEST_ERR(A3cm.RightScale(xDomain3cm),ierr);
00817   float A3cmOneNormFloat  = A3cm.NormOne();
00818   if (verbose1) cout << A3cm << endl;
00819   if (1.0!=A3cmOneNormFloat) {
00820     cout << "OneNorm should be = 1, but OneNorm = " << A3cmOneNormFloat << endl;
00821     EPETRA_TEST_ERR(-62,ierr);
00822     InvSumsBroke = true;
00823   }
00824 
00825   if (verbose) cout << "End partial sum testing" << endl;
00826   if (verbose) cout << "Begin replicated testing" << endl;
00827 
00828   // We will now view the shared row as a repliated row, rather than one
00829   // that has partial sums of its entries on mulitple processors.
00830   // We will reuse much of the data used for the partial sum tesitng.
00831   Epetra_Vector xRow3(RowMap3,false);
00832   Epetra_CrsMatrix A4(Copy, RowMap3, NumMyElements3);
00833   for (int ii=0; ii < NumMyElements3; ii++) {
00834     Values3[ii] = (int)((ii*.6)+1.0);
00835   }
00836   for (int ii=0; ii<NumMyRows3; ii++) {
00837     A4.InsertGlobalValues(myGlobalElements[ii],NumMyElements3,Values3,Indices3);
00838   }
00839   EPETRA_TEST_ERR(A4.FillComplete(DomainMap3, RangeMap3,false),ierr);
00840   if (verbose1) cout << A4 << endl;
00841   // The next two lines should be expanded into a verifiable test.
00842   EPETRA_TEST_ERR(A4.InvRowMaxs(xRow3),ierr);
00843   EPETRA_TEST_ERR(A4.InvRowMaxs(xRange3),ierr);
00844   if (verbose1) cout << xRow3 << xRange3;
00845 
00846   EPETRA_TEST_ERR(A4.InvRowSums(xRow3),ierr);
00847   if (verbose1) cout << xRow3;
00848   EPETRA_TEST_ERR(A4.LeftScale(xRow3),ierr);
00849   float A4infNormFloat = A4.NormInf();
00850   if (verbose1) cout << A4;
00851   if (2.0!=A4infNormFloat && NumProc != 1) {
00852     if (verbose1) cout << "InfNorm should be = 2 (because one column is replicated on two processors and NormOne() does not handle replication), but InfNorm = " << A4infNormFloat <<endl;
00853     EPETRA_TEST_ERR(-63,ierr);
00854     InvSumsBroke = true;
00855   }
00856   else if (1.0!=A4infNormFloat && NumProc == 1) {
00857     if (verbose1) cout << "InfNorm should be = 1, but InfNorm = " << A4infNormFloat <<endl;
00858     EPETRA_TEST_ERR(-63,ierr);
00859     InvSumsBroke = true;
00860   }
00861 
00862   Epetra_Vector xCol3cm(ColMap3cm,false);
00863   Epetra_CrsMatrix A4cm(Copy, RowMap3cm, ColMap3cm, NumProc+1);
00864   //Use values from A3cm
00865   for (int ii=0; ii<NumMyElements3; ii++) {
00866     A4cm.InsertGlobalValues(ii,NumMyColumns3,Values3cm,Indices3cm);
00867   }
00868   EPETRA_TEST_ERR(A4cm.FillComplete(DomainMap3cm, RangeMap3cm,false),ierr);
00869   if (verbose1) cout << A4cm << endl;
00870   // The next two lines should be expanded into a verifiable test.
00871   EPETRA_TEST_ERR(A4cm.InvColMaxs(xCol3cm),ierr);
00872   EPETRA_TEST_ERR(A4cm.InvColMaxs(xDomain3cm),ierr);
00873   if (verbose1) cout << xCol3cm << xDomain3cm;
00874 
00875   EPETRA_TEST_ERR(A4cm.InvColSums(xCol3cm),ierr);
00876 
00877   if (verbose1) cout << xCol3cm << endl;
00878   EPETRA_TEST_ERR(A4cm.RightScale(xCol3cm),ierr);
00879   float A4cmOneNormFloat = A4cm.NormOne();
00880   if (verbose1) cout << A4cm << endl;
00881   if (2.0!=A4cmOneNormFloat && NumProc != 1) {
00882     if (verbose1) cout << "OneNorm should be = 2 (because one column is replicated on two processors and NormOne() does not handle replication), but OneNorm = " << A4cmOneNormFloat << endl;
00883     EPETRA_TEST_ERR(-64,ierr);
00884     InvSumsBroke = true;
00885   }
00886   else if (1.0!=A4cmOneNormFloat && NumProc == 1) {
00887     if (verbose1) cout << "OneNorm should be = 1, but OneNorm = " << A4infNormFloat <<endl;
00888     EPETRA_TEST_ERR(-64,ierr);
00889     InvSumsBroke = true;
00890   }
00891 
00892   if (verbose) cout << "End replicated testing" << endl;
00893 
00894   if (InvSumsBroke) {
00895     if (verbose) cout << endl << "InvRowSums tests FAILED" << endl << endl;
00896   }
00897   else
00898     if (verbose) cout << endl << "InvRowSums tests PASSED" << endl << endl;
00899 
00900   A3cm.PutScalar(2.0);
00901   int nnz_A3cm = A3cm.Graph().NumGlobalNonzeros();
00902   double check_frobnorm = sqrt(nnz_A3cm*4.0);
00903   double frobnorm = A3cm.NormFrobenius();
00904 
00905   bool frobnorm_test_failed = false;
00906   if (fabs(check_frobnorm-frobnorm) > 5.e-5) {
00907     frobnorm_test_failed = true;
00908   }
00909 
00910   if (frobnorm_test_failed) {
00911     if (verbose) std::cout << "Frobenius-norm test FAILED."<<std::endl;
00912     EPETRA_TEST_ERR(-65, ierr);
00913   }
00914 
00915   // Subcommunicator test - only processor 1 has unknowns
00916   {
00917     int rv=0;
00918     int NumMyRows = (MyPID==0) ? NumGlobalEquations : 0;
00919     Epetra_Map Map1(-1,NumMyRows, 0, Comm);
00920 
00921     Epetra_CrsMatrix *A1 = new Epetra_CrsMatrix(Copy, Map1,0);
00922     double value = 1.0;
00923     for(int i=0; i<NumMyRows; i++) {
00924       int GID = Map1.GID(i);
00925       EPETRA_TEST_ERR(A1->InsertGlobalValues(GID, 1,&value,&GID),ierr);
00926     }
00927     EPETRA_TEST_ERR(A1->FillComplete(),ierr);
00928 
00929     Epetra_BlockMap *Map2 = Map1.RemoveEmptyProcesses();
00930     rv=A1->RemoveEmptyProcessesInPlace(Map2);
00931     if(rv!=0) {
00932       if (verbose) std::cout << "Subcommunicator test FAILED."<<std::endl;
00933       EPETRA_TEST_ERR(-66, ierr);
00934     }
00935 
00936     delete Map2;
00937   }
00938 
00939   delete [] Values2;
00940   delete [] Indices2;
00941   delete [] myGlobalElements;
00942   delete [] Values3;
00943   delete [] Indices3;
00944   delete [] Values3cm;
00945   delete [] Indices3cm;
00946   delete [] RangeLeftScaleValues;
00947   delete [] RowLeftScaleValues;
00948 #ifdef EPETRA_MPI
00949   MPI_Finalize() ;
00950 #endif
00951 
00952 /* end main
00953 */
00954 return ierr ;
00955 }
00956 
00957 int power_method(bool TransA, Epetra_CrsMatrix& A, Epetra_Vector& q, Epetra_Vector& z,
00958                  Epetra_Vector& resid, double* lambda, int niters, double tolerance, bool verbose)
00959 {
00960   
00961   // Fill z with random Numbers
00962   z.Random();
00963   
00964   // variable needed for iteration
00965   double normz, residual;
00966 
00967   int ierr = 1;
00968   
00969   for(int iter = 0; iter < niters; iter++) {
00970     z.Norm2(&normz); // Compute 2-norm of z
00971     q.Scale(1.0/normz, z);
00972     A.Multiply(TransA, q, z); // Compute z = A*q // SEGFAULT HAPPENS HERE
00973     q.Dot(z, lambda); // Approximate maximum eigenvaluE
00974     if(iter%100==0 || iter+1==niters) {
00975       resid.Update(1.0, z, -(*lambda), q, 0.0); // Compute A*q - lambda*q
00976       resid.Norm2(&residual);
00977       if(verbose) cout << "Iter = " << iter << "  Lambda = " << *lambda
00978                        << "  Residual of A*q - lambda*q = " << residual << endl;
00979     }
00980     if(residual < tolerance) {
00981       ierr = 0;
00982       break;
00983     }
00984   }
00985   return(ierr);
00986 }
00987 
00988 int check(Epetra_CrsMatrix& A, int NumMyRows1, int NumGlobalRows1, int NumMyNonzeros1,
00989           int NumGlobalNonzeros1, int* MyGlobalElements, bool verbose)
00990 {
00991   (void)MyGlobalElements;
00992   int ierr = 0, forierr = 0;
00993   int NumGlobalIndices;
00994   int NumMyIndices;
00995   int* MyViewIndices = 0;
00996   int* GlobalViewIndices = 0;
00997   double* MyViewValues = 0;
00998   double* GlobalViewValues = 0;
00999   int MaxNumIndices = A.Graph().MaxNumIndices();
01000   int* MyCopyIndices = new int[MaxNumIndices];
01001   int* GlobalCopyIndices = new int[MaxNumIndices];
01002   double* MyCopyValues = new double[MaxNumIndices];
01003   double* GlobalCopyValues = new double[MaxNumIndices];
01004 
01005   // Test query functions
01006 
01007   int NumMyRows = A.NumMyRows();
01008   if (verbose) cout << "\n\nNumber of local Rows = " << NumMyRows << endl<< endl;
01009 
01010   EPETRA_TEST_ERR(!(NumMyRows==NumMyRows1),ierr);
01011 
01012   int NumMyNonzeros = A.NumMyNonzeros();
01013   if (verbose) cout << "\n\nNumber of local Nonzero entries = " << NumMyNonzeros << endl<< endl;
01014 
01015   EPETRA_TEST_ERR(!(NumMyNonzeros==NumMyNonzeros1),ierr);
01016 
01017   int NumGlobalRows = A.NumGlobalRows();
01018   if (verbose) cout << "\n\nNumber of global Rows = " << NumGlobalRows << endl<< endl;
01019 
01020   EPETRA_TEST_ERR(!(NumGlobalRows==NumGlobalRows1),ierr);
01021 
01022   int NumGlobalNonzeros = A.NumGlobalNonzeros();
01023   if (verbose) cout << "\n\nNumber of global Nonzero entries = " << NumGlobalNonzeros << endl<< endl;
01024 
01025   EPETRA_TEST_ERR(!(NumGlobalNonzeros==NumGlobalNonzeros1),ierr);
01026 
01027   // GlobalRowView should be illegal (since we have local indices)
01028 
01029   EPETRA_TEST_ERR(!(A.ExtractGlobalRowView(A.RowMap().MaxMyGID(), NumGlobalIndices, GlobalViewValues, GlobalViewIndices)==-2),ierr);
01030 
01031   // Other binary tests
01032 
01033   EPETRA_TEST_ERR(A.NoDiagonal(),ierr);
01034   EPETRA_TEST_ERR(!(A.Filled()),ierr);
01035   EPETRA_TEST_ERR(!(A.MyGRID(A.RowMap().MaxMyGID())),ierr);
01036   EPETRA_TEST_ERR(!(A.MyGRID(A.RowMap().MinMyGID())),ierr);
01037   EPETRA_TEST_ERR(A.MyGRID(1+A.RowMap().MaxMyGID()),ierr);
01038   EPETRA_TEST_ERR(A.MyGRID(-1+A.RowMap().MinMyGID()),ierr);
01039   EPETRA_TEST_ERR(!(A.MyLRID(0)),ierr);
01040   EPETRA_TEST_ERR(!(A.MyLRID(NumMyRows-1)),ierr);
01041   EPETRA_TEST_ERR(A.MyLRID(-1),ierr);
01042   EPETRA_TEST_ERR(A.MyLRID(NumMyRows),ierr);
01043 
01044   forierr = 0;
01045   for (int i = 0; i < NumMyRows; i++) {
01046     int Row = A.GRID(i);
01047     A.ExtractGlobalRowCopy(Row, MaxNumIndices, NumGlobalIndices, GlobalCopyValues, GlobalCopyIndices);
01048     A.ExtractMyRowView(i, NumMyIndices, MyViewValues, MyViewIndices); // this is where the problem comes from
01049     forierr += !(NumGlobalIndices == NumMyIndices);
01050     for(int j = 1; j < NumMyIndices; j++) {
01051       forierr += !(MyViewIndices[j-1] < MyViewIndices[j]); // this is where the test fails
01052     }
01053     for(int j = 0; j < NumGlobalIndices; j++) {
01054       forierr += !(GlobalCopyIndices[j] == A.GCID(MyViewIndices[j]));
01055       forierr += !(A.LCID(GlobalCopyIndices[j]) == MyViewIndices[j]);
01056       forierr += !(GlobalCopyValues[j] == MyViewValues[j]);
01057     }
01058   }
01059   EPETRA_TEST_ERR(forierr,ierr);
01060 
01061   forierr = 0;
01062   for (int i = 0; i < NumMyRows; i++) {
01063     int Row = A.GRID(i);
01064     A.ExtractGlobalRowCopy(Row, MaxNumIndices, NumGlobalIndices, GlobalCopyValues, GlobalCopyIndices);
01065     A.ExtractMyRowCopy(i, MaxNumIndices, NumMyIndices, MyCopyValues, MyCopyIndices);
01066     forierr += !(NumGlobalIndices == NumMyIndices);
01067     for (int j = 1; j < NumMyIndices; j++)
01068       forierr += !(MyCopyIndices[j-1] < MyCopyIndices[j]);
01069     for (int j = 0; j < NumGlobalIndices; j++) {
01070       forierr += !(GlobalCopyIndices[j] == A.GCID(MyCopyIndices[j]));
01071       forierr += !(A.LCID(GlobalCopyIndices[j]) == MyCopyIndices[j]);
01072       forierr += !(GlobalCopyValues[j] == MyCopyValues[j]);
01073     }
01074 
01075   }
01076   EPETRA_TEST_ERR(forierr,ierr);
01077 
01078   delete [] MyCopyIndices;
01079   delete [] GlobalCopyIndices;
01080   delete [] MyCopyValues;
01081   delete [] GlobalCopyValues;
01082 
01083   if (verbose) cout << "\n\nRows sorted check OK" << endl<< endl;
01084 
01085   return (ierr);
01086 }
01087 
01088 int check_graph_sharing(Epetra_Comm& Comm)
01089 {
01090   int numLocalElems = 5;
01091   int localProc = Comm.MyPID();
01092   int firstElem = localProc*numLocalElems;
01093   int err;
01094   Epetra_Map map(-1, numLocalElems, 0, Comm);
01095 
01096   Epetra_CrsMatrix* A = new Epetra_CrsMatrix(Copy, map, 1);
01097 
01098   for (int i=0; i<numLocalElems; ++i) {
01099     int row = firstElem+i;
01100     int col = row;
01101     double val = 1.0;
01102 
01103     err = A->InsertGlobalValues(row, 1, &val, &col);
01104     if (err != 0) {
01105       cerr << "A->InsertGlobalValues("<<row<<") returned err="<<err<<endl;
01106       return(err);
01107     }
01108   }
01109 
01110   A->FillComplete(false);
01111 
01112   Epetra_CrsMatrix B(Copy, A->Graph());
01113 
01114   delete A;
01115 
01116   for (int i=0; i<numLocalElems; ++i) {
01117     int row = firstElem+i;
01118     int col = row;
01119     double val = 1.0;
01120 
01121     err = B.ReplaceGlobalValues(row, 1, &val, &col);
01122     if (err != 0) {
01123       cerr << "B.InsertGlobalValues("<<row<<") returned err="<<err<<endl;
01124       return(err);
01125     }
01126   }
01127 
01128   return(0);
01129 }
01130 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines