test/CrsMatrix/cxx_main.cpp

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

Generated on Thu Sep 18 12:37:56 2008 for Epetra Package Browser (Single Doxygen Collection) by doxygen 1.3.9.1