Epetra Package Browser (Single Doxygen Collection) Development
test/RowMatrix/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 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_JadMatrix.h"
00047 #include "Epetra_Vector.h"
00048 #include "Epetra_Flops.h"
00049 #ifdef EPETRA_MPI
00050 #include "Epetra_MpiComm.h"
00051 #include "mpi.h"
00052 #else
00053 #include "Epetra_SerialComm.h"
00054 #endif
00055 #include "../epetra_test_err.h"
00056 #include "Epetra_Version.h"
00057 #include <vector>
00058 #include <algorithm>
00059 #include <string>
00060 
00061 // prototypes
00062 
00063 int checkValues( double x, double y, string message = "", bool verbose = false) { 
00064   if (fabs((x-y)/x) > 0.01) {
00065     return(1); 
00066     if (verbose) cout << "********** " << message << " check failed.********** " << endl;
00067   }
00068   else {
00069     if (verbose) cout << message << " check OK." << endl;    
00070     return(0);
00071   }
00072 }
00073 
00074 int checkMultiVectors( Epetra_MultiVector & X, Epetra_MultiVector & Y, string message = "", bool verbose = false) {
00075   int numVectors = X.NumVectors();
00076   int length = Y.MyLength();
00077   int badvalue = 0;
00078   int globalbadvalue = 0;
00079   for (int j=0; j<numVectors; j++)
00080     for (int i=0; i< length; i++) 
00081       if (checkValues(X[j][i], Y[j][i])==1) badvalue = 1;
00082   X.Map().Comm().MaxAll(&badvalue, &globalbadvalue, 1);
00083 
00084   if (verbose) {
00085     if (globalbadvalue==0) cout << message << " check OK." << endl;
00086     else cout << "********* " << message << " check failed.********** " << endl;
00087   }
00088   return(globalbadvalue);
00089 }
00090 
00091 int check(Epetra_RowMatrix & A, Epetra_RowMatrix & B, bool verbose);
00092 
00093 int powerMethodTests(Epetra_RowMatrix & A, Epetra_RowMatrix & JadA, Epetra_Map & Map, 
00094          Epetra_Vector & q, Epetra_Vector & z, Epetra_Vector & resid, bool verbose);
00095 
00096 int power_method(bool TransA, Epetra_RowMatrix& A, 
00097      Epetra_Vector& q,
00098      Epetra_Vector& z0, 
00099      Epetra_Vector& resid, 
00100      double * lambda, int niters, double tolerance,
00101      bool verbose);
00102 
00103 int main(int argc, char *argv[])
00104 {
00105   int ierr = 0, i, forierr = 0;
00106 #ifdef EPETRA_MPI
00107 
00108   // Initialize MPI
00109 
00110   MPI_Init(&argc,&argv);
00111   int rank; // My process ID
00112 
00113   MPI_Comm_rank(MPI_COMM_WORLD, &rank);
00114   Epetra_MpiComm Comm( MPI_COMM_WORLD );
00115 
00116 #else
00117 
00118   int rank = 0;
00119   Epetra_SerialComm Comm;
00120 
00121 #endif
00122 
00123   bool verbose = false;
00124 
00125   // Check if we should print results to standard out
00126   if (argc>1) if (argv[1][0]=='-' && argv[1][1]=='v') verbose = true;
00127 
00128   int verbose_int = verbose ? 1 : 0;
00129   Comm.Broadcast(&verbose_int, 1, 0);
00130   verbose = verbose_int==1 ? true : false;
00131 
00132 
00133   //  char tmp;
00134   //  if (rank==0) cout << "Press any key to continue..."<< endl;
00135   //  if (rank==0) cin >> tmp;
00136   //  Comm.Barrier();
00137 
00138   Comm.SetTracebackMode(0); // This should shut down any error traceback reporting
00139   int MyPID = Comm.MyPID();
00140   int NumProc = Comm.NumProc();
00141 
00142   if(verbose && MyPID==0)
00143     cout << Epetra_Version() << endl << endl;
00144 
00145   if (verbose) cout << "Processor "<<MyPID<<" of "<< NumProc
00146         << " is alive."<<endl;
00147 
00148   // Redefine verbose to only print on PE 0
00149   if(verbose && rank!=0) 
00150     verbose = false;
00151 
00152   int NumMyEquations = 10000;
00153   int NumGlobalEquations = (NumMyEquations * NumProc) + EPETRA_MIN(NumProc,3);
00154   if(MyPID < 3) 
00155     NumMyEquations++;
00156 
00157   // Construct a Map that puts approximately the same Number of equations on each processor
00158 
00159   Epetra_Map Map(NumGlobalEquations, NumMyEquations, 0, Comm);
00160   
00161   // Get update list and number of local equations from newly created Map
00162   vector<int> MyGlobalElements(Map.NumMyElements());
00163   Map.MyGlobalElements(&MyGlobalElements[0]);
00164 
00165   // Create an integer vector NumNz that is used to build the Petra Matrix.
00166   // NumNz[i] is the Number of OFF-DIAGONAL term for the ith global equation on this processor
00167 
00168   vector<int> NumNz(NumMyEquations);
00169 
00170   // We are building a tridiagonal matrix where each row has (-1 2 -1)
00171   // So we need 2 off-diagonal terms (except for the first and last equation)
00172 
00173   for(i = 0; i < NumMyEquations; i++)
00174     if((MyGlobalElements[i] == 0) || (MyGlobalElements[i] == NumGlobalEquations - 1))
00175       NumNz[i] = 1;
00176     else
00177       NumNz[i] = 2;
00178 
00179   // Create a Epetra_Matrix
00180 
00181   Epetra_CrsMatrix A(Copy, Map, &NumNz[0]);
00182   EPETRA_TEST_ERR(A.IndicesAreGlobal(),ierr);
00183   EPETRA_TEST_ERR(A.IndicesAreLocal(),ierr);
00184   
00185   // Add  rows one-at-a-time
00186   // Need some vectors to help
00187   // Off diagonal Values will always be -1
00188 
00189 
00190   vector<double> Values(2);
00191   Values[0] = -1.0; 
00192   Values[1] = -1.0;
00193   vector<int> Indices(2);
00194   double two = 2.0;
00195   int NumEntries;
00196 
00197   forierr = 0;
00198   for(i = 0; i < NumMyEquations; i++) {
00199     if(MyGlobalElements[i] == 0) {
00200       Indices[0] = 1;
00201       NumEntries = 1;
00202     }
00203     else if (MyGlobalElements[i] == NumGlobalEquations-1) {
00204       Indices[0] = NumGlobalEquations-2;
00205       NumEntries = 1;
00206     }
00207     else {
00208       Indices[0] = MyGlobalElements[i]-1;
00209       Indices[1] = MyGlobalElements[i]+1;
00210       NumEntries = 2;
00211     }
00212     forierr += !(A.InsertGlobalValues(MyGlobalElements[i], NumEntries, &Values[0], &Indices[0])==0);
00213     forierr += !(A.InsertGlobalValues(MyGlobalElements[i], 1, &two, &MyGlobalElements[i])>0); // Put in the diagonal entry
00214   }
00215   EPETRA_TEST_ERR(forierr,ierr);
00216 
00217   // Finish up
00218   A.FillComplete();
00219   A.OptimizeStorage();
00220 
00221   Epetra_JadMatrix JadA(A);
00222   Epetra_JadMatrix JadA1(A);
00223   Epetra_JadMatrix JadA2(A);
00224 
00225   // Create vectors for Power method
00226 
00227   Epetra_Vector q(Map);
00228   Epetra_Vector z(Map); z.Random();
00229   Epetra_Vector resid(Map);
00230 
00231   Epetra_Flops flopcounter;
00232   A.SetFlopCounter(flopcounter);
00233   q.SetFlopCounter(A);
00234   z.SetFlopCounter(A);
00235   resid.SetFlopCounter(A);
00236   JadA.SetFlopCounter(A);
00237   JadA1.SetFlopCounter(A);
00238   JadA2.SetFlopCounter(A);
00239   
00240 
00241   if (verbose) cout << "=======================================" << endl
00242         << "Testing Jad using CrsMatrix as input..." << endl
00243         << "=======================================" << endl;
00244 
00245   A.ResetFlops();
00246   powerMethodTests(A, JadA, Map, q, z, resid, verbose);
00247 
00248   // Increase diagonal dominance
00249 
00250   if (verbose) cout << "\n\nIncreasing the magnitude of first diagonal term and solving again\n\n"
00251         << endl;
00252 
00253   
00254   if (A.MyGlobalRow(0)) {
00255     int numvals = A.NumGlobalEntries(0);
00256     vector<double> Rowvals(numvals);
00257     vector<int> Rowinds(numvals);
00258     A.ExtractGlobalRowCopy(0, numvals, numvals, &Rowvals[0], &Rowinds[0]); // Get A[0,0]
00259 
00260     for (i=0; i<numvals; i++) if (Rowinds[i] == 0) Rowvals[i] *= 10.0;
00261     
00262     A.ReplaceGlobalValues(0, numvals, &Rowvals[0], &Rowinds[0]);
00263   }
00264   JadA.UpdateValues(A);
00265   A.ResetFlops();
00266   powerMethodTests(A, JadA, Map, q, z, resid, verbose);
00267 
00268   if (verbose) cout << "================================================================" << endl
00269               << "Testing Jad using Jad matrix as input matrix for construction..." << endl
00270               << "================================================================" << endl;
00271   JadA1.ResetFlops();
00272   powerMethodTests(JadA1, JadA2, Map, q, z, resid, verbose);
00273 
00274 #ifdef EPETRA_MPI
00275   MPI_Finalize() ;
00276 #endif
00277 
00278 return ierr ;
00279 }
00280 
00281 int powerMethodTests(Epetra_RowMatrix & A, Epetra_RowMatrix & JadA, Epetra_Map & Map, 
00282          Epetra_Vector & q, Epetra_Vector & z, Epetra_Vector & resid, bool verbose) {
00283 
00284   // variable needed for iteration
00285   double lambda = 0.0;
00286   // int niters = 10000;
00287   int niters = 300;
00288   double tolerance = 1.0e-2;
00289   int ierr = 0;
00290 
00292   
00293   // Iterate
00294 
00295   Epetra_Time timer(Map.Comm());
00296   
00297   double startTime = timer.ElapsedTime();
00298   EPETRA_TEST_ERR(power_method(false, A, q, z, resid, &lambda, niters, tolerance, verbose),ierr);
00299   double elapsed_time = timer.ElapsedTime() - startTime;
00300   double total_flops = q.Flops();
00301   double MFLOPs = total_flops/elapsed_time/1000000.0;
00302   double lambdaref = lambda;
00303   double flopsref = total_flops;
00304 
00305   if (verbose) 
00306     cout << "\n\nTotal MFLOPs for reference first solve = " << MFLOPs << endl
00307       <<     "Total FLOPS                            = " <<total_flops <<endl<<endl;
00308 
00309   lambda = 0.0;
00310   startTime = timer.ElapsedTime();
00311   EPETRA_TEST_ERR(power_method(false, JadA, q, z, resid, &lambda, niters, tolerance, verbose),ierr);
00312   elapsed_time = timer.ElapsedTime() - startTime;
00313   total_flops = q.Flops();
00314   MFLOPs = total_flops/elapsed_time/1000000.0;
00315 
00316   if (verbose) 
00317     cout << "\n\nTotal MFLOPs for candidate first solve = " << MFLOPs << endl
00318       <<     "Total FLOPS                            = " <<total_flops <<endl<<endl;
00319 
00320   EPETRA_TEST_ERR(checkValues(lambda,lambdaref," No-transpose Power Method result", verbose),ierr);
00321   EPETRA_TEST_ERR(checkValues(total_flops,flopsref," No-transpose Power Method flop count", verbose),ierr);
00322 
00324   
00325   // Solve transpose problem
00326 
00327   if (verbose) cout << "\n\nUsing transpose of matrix and solving again (should give same result).\n\n"
00328         << endl;
00329   // Iterate
00330   lambda = 0.0;
00331   startTime = timer.ElapsedTime();
00332   EPETRA_TEST_ERR(power_method(true, A, q, z, resid, &lambda, niters, tolerance, verbose),ierr);
00333   elapsed_time = timer.ElapsedTime() - startTime;
00334   total_flops = q.Flops();
00335   MFLOPs = total_flops/elapsed_time/1000000.0;
00336   lambdaref = lambda;
00337   flopsref = total_flops;
00338 
00339   if (verbose) 
00340    cout << "\n\nTotal MFLOPs for reference transpose solve = " << MFLOPs << endl
00341      <<     "Total FLOPS                                = " <<total_flops <<endl<<endl;
00342 
00343   lambda = 0.0;
00344   startTime = timer.ElapsedTime();
00345   EPETRA_TEST_ERR(power_method(true, JadA, q, z, resid, &lambda, niters, tolerance, verbose),ierr);
00346   elapsed_time = timer.ElapsedTime() - startTime;
00347   total_flops = q.Flops();
00348   MFLOPs = total_flops/elapsed_time/1000000.0;
00349 
00350   if (verbose) 
00351     cout << "\n\nTotal MFLOPs for candidate transpose solve = " << MFLOPs << endl
00352       <<     "Total FLOPS                                = " <<total_flops <<endl<<endl;
00353 
00354   EPETRA_TEST_ERR(checkValues(lambda,lambdaref,"Transpose Power Method result", verbose),ierr);
00355   EPETRA_TEST_ERR(checkValues(total_flops,flopsref,"Transpose Power Method flop count", verbose),ierr);
00356 
00357   EPETRA_TEST_ERR(check(A, JadA, verbose),ierr);
00358 
00359   return(0);
00360 }
00361 int power_method(bool TransA, Epetra_RowMatrix& A, Epetra_Vector& q, Epetra_Vector& z0, 
00362      Epetra_Vector& resid, double* lambda, int niters, double tolerance, bool verbose) 
00363 {  
00364   
00365   // Fill z with random Numbers
00366   Epetra_Vector z(z0);
00367   
00368   // variable needed for iteration
00369   double normz, residual;
00370 
00371   int ierr = 1;
00372   
00373   for(int iter = 0; iter < niters; iter++) {
00374     z.Norm2(&normz); // Compute 2-norm of z
00375     q.Scale(1.0/normz, z);
00376     A.Multiply(TransA, q, z); // Compute z = A*q // SEGFAULT HAPPENS HERE
00377     q.Dot(z, lambda); // Approximate maximum eigenvaluE
00378     if(iter%100==0 || iter+1==niters) {
00379       resid.Update(1.0, z, -(*lambda), q, 0.0); // Compute A*q - lambda*q
00380       resid.Norm2(&residual);
00381       if(verbose) cout << "Iter = " << iter << "  Lambda = " << *lambda 
00382                        << "  Residual of A*q - lambda*q = " << residual << endl;
00383     }
00384     if(residual < tolerance) {
00385       ierr = 0;
00386       break;
00387     }
00388   }
00389   return(ierr);
00390 }
00391 
00392 int check(Epetra_RowMatrix& A, Epetra_RowMatrix & B, bool verbose)  {  
00393 
00394   int ierr = 0;
00395   EPETRA_TEST_ERR(!A.Comm().NumProc()==B.Comm().NumProc(),ierr);
00396   EPETRA_TEST_ERR(!A.Comm().MyPID()==B.Comm().MyPID(),ierr);
00397   EPETRA_TEST_ERR(!A.Filled()==B.Filled(),ierr);
00398   EPETRA_TEST_ERR(!A.HasNormInf()==B.HasNormInf(),ierr);
00399   EPETRA_TEST_ERR(!A.LowerTriangular()==B.LowerTriangular(),ierr);
00400   EPETRA_TEST_ERR(!A.Map().SameAs(B.Map()),ierr);
00401   EPETRA_TEST_ERR(!A.MaxNumEntries()==B.MaxNumEntries(),ierr);
00402   EPETRA_TEST_ERR(!A.NumGlobalCols()==B.NumGlobalCols(),ierr);
00403   EPETRA_TEST_ERR(!A.NumGlobalDiagonals()==B.NumGlobalDiagonals(),ierr);
00404   EPETRA_TEST_ERR(!A.NumGlobalNonzeros()==B.NumGlobalNonzeros(),ierr);
00405   EPETRA_TEST_ERR(!A.NumGlobalRows()==B.NumGlobalRows(),ierr);
00406   EPETRA_TEST_ERR(!A.NumMyCols()==B.NumMyCols(),ierr);
00407   EPETRA_TEST_ERR(!A.NumMyDiagonals()==B.NumMyDiagonals(),ierr);
00408   EPETRA_TEST_ERR(!A.NumMyNonzeros()==B.NumMyNonzeros(),ierr);
00409   for (int i=0; i<A.NumMyRows(); i++) {
00410     int nA, nB;
00411     A.NumMyRowEntries(i,nA); B.NumMyRowEntries(i,nB);
00412     EPETRA_TEST_ERR(!nA==nB,ierr);
00413   }
00414   EPETRA_TEST_ERR(!A.NumMyRows()==B.NumMyRows(),ierr);
00415   EPETRA_TEST_ERR(!A.OperatorDomainMap().SameAs(B.OperatorDomainMap()),ierr);
00416   EPETRA_TEST_ERR(!A.OperatorRangeMap().SameAs(B.OperatorRangeMap()),ierr);
00417   EPETRA_TEST_ERR(!A.RowMatrixColMap().SameAs(B.RowMatrixColMap()),ierr);
00418   EPETRA_TEST_ERR(!A.RowMatrixRowMap().SameAs(B.RowMatrixRowMap()),ierr);
00419   EPETRA_TEST_ERR(!A.UpperTriangular()==B.UpperTriangular(),ierr);
00420   EPETRA_TEST_ERR(!A.UseTranspose()==B.UseTranspose(),ierr);
00421 
00422   int NumVectors = 5;
00423   { // No transpose case
00424     Epetra_MultiVector X(A.OperatorDomainMap(), NumVectors);
00425     Epetra_MultiVector YA1(A.OperatorRangeMap(), NumVectors);
00426     Epetra_MultiVector YA2(YA1);
00427     Epetra_MultiVector YB1(YA1);
00428     Epetra_MultiVector YB2(YA1);
00429     X.Random();
00430     
00431     bool transA = false;
00432     A.SetUseTranspose(transA);
00433     B.SetUseTranspose(transA);
00434     A.Apply(X,YA1);
00435     A.Multiply(transA, X, YA2);
00436     EPETRA_TEST_ERR(checkMultiVectors(YA1,YA2,"A Multiply and A Apply", verbose),ierr);
00437     B.Apply(X,YB1);
00438     EPETRA_TEST_ERR(checkMultiVectors(YA1,YB1,"A Multiply and B Multiply", verbose),ierr);
00439     B.Multiply(transA, X, YB2);
00440     EPETRA_TEST_ERR(checkMultiVectors(YA1,YB2,"A Multiply and B Apply", verbose), ierr);
00441     
00442   }
00443   {// transpose case
00444     Epetra_MultiVector X(A.OperatorRangeMap(), NumVectors);
00445     Epetra_MultiVector YA1(A.OperatorDomainMap(), NumVectors);
00446     Epetra_MultiVector YA2(YA1);
00447     Epetra_MultiVector YB1(YA1);
00448     Epetra_MultiVector YB2(YA1);
00449     X.Random();
00450     
00451     bool transA = true;
00452     A.SetUseTranspose(transA);
00453     B.SetUseTranspose(transA);
00454     A.Apply(X,YA1);
00455     A.Multiply(transA, X, YA2);
00456     EPETRA_TEST_ERR(checkMultiVectors(YA1,YA2, "A Multiply and A Apply (transpose)", verbose),ierr);
00457     B.Apply(X,YB1);
00458     EPETRA_TEST_ERR(checkMultiVectors(YA1,YB1, "A Multiply and B Multiply (transpose)", verbose),ierr);
00459     B.Multiply(transA, X,YB2);
00460     EPETRA_TEST_ERR(checkMultiVectors(YA1,YB2, "A Multiply and B Apply (transpose)", verbose),ierr);
00461     
00462   }
00463 
00464   Epetra_Vector diagA(A.RowMatrixRowMap());
00465   EPETRA_TEST_ERR(A.ExtractDiagonalCopy(diagA),ierr);
00466   Epetra_Vector diagB(B.RowMatrixRowMap());
00467   EPETRA_TEST_ERR(B.ExtractDiagonalCopy(diagB),ierr);
00468   EPETRA_TEST_ERR(checkMultiVectors(diagA,diagB, "ExtractDiagonalCopy", verbose),ierr);
00469 
00470   Epetra_Vector rowA(A.RowMatrixRowMap());
00471   EPETRA_TEST_ERR(A.InvRowSums(rowA),ierr);
00472   Epetra_Vector rowB(B.RowMatrixRowMap());
00473   EPETRA_TEST_ERR(B.InvRowSums(rowB),ierr)
00474   EPETRA_TEST_ERR(checkMultiVectors(rowA,rowB, "InvRowSums", verbose),ierr);
00475 
00476   Epetra_Vector colA(A.RowMatrixColMap());
00477   EPETRA_TEST_ERR(A.InvColSums(colA),ierr);
00478   Epetra_Vector colB(B.RowMatrixColMap());
00479   EPETRA_TEST_ERR(B.InvColSums(colB),ierr);
00480   EPETRA_TEST_ERR(checkMultiVectors(colA,colB, "InvColSums", verbose),ierr);
00481 
00482   EPETRA_TEST_ERR(checkValues(A.NormInf(), B.NormInf(), "NormInf before scaling", verbose), ierr);
00483   EPETRA_TEST_ERR(checkValues(A.NormOne(), B.NormOne(), "NormOne before scaling", verbose),ierr);
00484 
00485   EPETRA_TEST_ERR(A.RightScale(colA),ierr);  
00486   EPETRA_TEST_ERR(B.RightScale(colB),ierr);
00487 
00488 
00489   EPETRA_TEST_ERR(A.LeftScale(rowA),ierr);
00490   EPETRA_TEST_ERR(B.LeftScale(rowB),ierr);
00491 
00492 
00493   EPETRA_TEST_ERR(checkValues(A.NormInf(), B.NormInf(), "NormInf after scaling", verbose), ierr);
00494   EPETRA_TEST_ERR(checkValues(A.NormOne(), B.NormOne(), "NormOne after scaling", verbose),ierr);
00495 
00496   vector<double> valuesA(A.MaxNumEntries());
00497   vector<int> indicesA(A.MaxNumEntries());  
00498   vector<double> valuesB(B.MaxNumEntries());
00499   vector<int> indicesB(B.MaxNumEntries());
00500   return(0);
00501   for (int i=0; i<A.NumMyRows(); i++) {
00502     int nA, nB;
00503     EPETRA_TEST_ERR(A.ExtractMyRowCopy(i, A.MaxNumEntries(), nA, &valuesA[0], &indicesA[0]),ierr); 
00504     EPETRA_TEST_ERR(B.ExtractMyRowCopy(i, B.MaxNumEntries(), nB, &valuesB[0], &indicesB[0]),ierr);
00505     EPETRA_TEST_ERR(!nA==nB,ierr);
00506     for (int j=0; j<nA; j++) {
00507       double curVal = valuesA[j];
00508       int curIndex = indicesA[j];
00509       bool notfound = true;
00510       int jj = 0;
00511       while (notfound && jj< nB) {
00512   if (!checkValues(curVal, valuesB[jj])) notfound = false;
00513   jj++;
00514       }
00515       EPETRA_TEST_ERR(notfound, ierr);
00516       vector<int>::iterator p = find(indicesB.begin(),indicesB.end(),curIndex);  // find curIndex in indicesB
00517       EPETRA_TEST_ERR(p==indicesB.end(), ierr);
00518     }
00519 
00520   }
00521   if (verbose) cout << "RowMatrix Methods check OK" << endl;
00522 
00523   return (ierr);
00524 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines