EpetraExt Package Browser (Single Doxygen Collection) Development
test/hypre/cxx_main.cpp
Go to the documentation of this file.
00001 //@HEADER
00002 // ***********************************************************************
00003 //
00004 //     EpetraExt: Epetra Extended - Linear Algebra Services Package
00005 //                 Copyright (2011) Sandia Corporation
00006 //
00007 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00008 // the U.S. Government retains certain rights in this software.
00009 //
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00038 //
00039 // ***********************************************************************
00040 //@HEADER
00041 
00042 #include "Epetra_Map.h"
00043 #include "Epetra_Time.h"
00044 #include "Epetra_CrsMatrix.h"
00045 #include "EpetraExt_HypreIJMatrix.h"
00046 #include "HYPRE_IJ_mv.h"
00047 #include "hypre_Helpers.hpp"
00048 #include "Epetra_Vector.h"
00049 #include "Epetra_Flops.h"
00050 #ifdef EPETRA_MPI
00051 #include "Epetra_MpiComm.h"
00052 #include "mpi.h"
00053 #else
00054 #include "Epetra_SerialComm.h"
00055 #endif
00056 #include "../epetra_test_err.h"
00057 #include "Epetra_Version.h"
00058 #include <vector>
00059 #include <algorithm>
00060 #include <string>
00061 #include <iostream>
00062 #include <fstream>
00063 
00064 // prototypes
00065 
00066 int checkValues( double x, double y, string message = "", bool verbose = false) { 
00067   if (fabs((x-y)/x) > 0.01) {
00068     return(1); 
00069     if (verbose) cout << "********** " << message << " check failed.********** " << endl;
00070   }
00071   else {
00072     if (verbose) cout << message << " check OK." << endl;    
00073     return(0);
00074   }
00075 }
00076 
00077 int checkMultiVectors( Epetra_MultiVector & X, Epetra_MultiVector & Y, string message = "", bool verbose = false) {
00078   int numVectors = X.NumVectors();
00079   int length = Y.MyLength();
00080   int badvalue = 0;
00081   int globalbadvalue = 0;
00082   for (int j=0; j<numVectors; j++)
00083     for (int i=0; i< length; i++) 
00084       if (checkValues(X[j][i], Y[j][i])==1) badvalue = 1;
00085   X.Map().Comm().MaxAll(&badvalue, &globalbadvalue, 1);
00086 
00087   if (verbose) {
00088     if (globalbadvalue==0) cout << message << " check OK." << endl;
00089     else cout << "********* " << message << " check failed.********** " << endl;
00090   }
00091   return(globalbadvalue);
00092 }
00093 
00094 int check(Epetra_RowMatrix & A, Epetra_RowMatrix & B, bool verbose);
00095 
00096 int powerMethodTests(Epetra_RowMatrix & A, Epetra_RowMatrix & JadA, Epetra_Map & Map, 
00097          Epetra_Vector & q, Epetra_Vector & z, Epetra_Vector & resid, bool verbose);
00098 
00099 int power_method(bool TransA, Epetra_RowMatrix& A, 
00100      Epetra_Vector& q,
00101      Epetra_Vector& z0, 
00102      Epetra_Vector& resid, 
00103      double * lambda, int niters, double tolerance,
00104      bool verbose);
00105 
00106 int main(int argc, char *argv[])
00107 {
00108   int ierr = 0, i, forierr = 0;
00109 #ifdef EPETRA_MPI
00110 
00111   // Initialize MPI
00112 
00113   MPI_Init(&argc,&argv);
00114   int rank; // My process ID
00115 
00116   MPI_Comm_rank(MPI_COMM_WORLD, &rank);
00117   Epetra_MpiComm Comm( MPI_COMM_WORLD );
00118 
00119 #else
00120 
00121   int rank = 0;
00122   Epetra_SerialComm Comm;
00123 
00124 #endif
00125 
00126   bool verbose = false;
00127 
00128   // Check if we should print results to standard out
00129   if (argc>1) if (argv[1][0]=='-' && argv[1][1]=='v') verbose = true;
00130 
00131   int verbose_int = verbose ? 1 : 0;
00132   Comm.Broadcast(&verbose_int, 1, 0);
00133   verbose = verbose_int==1 ? true : false;
00134 
00135 
00136   //  char tmp;
00137   //  if (rank==0) cout << "Press any key to continue..."<< endl;
00138   //  if (rank==0) cin >> tmp;
00139   //  Comm.Barrier();
00140 
00141   Comm.SetTracebackMode(0); // This should shut down any error traceback reporting
00142   int MyPID = Comm.MyPID();
00143   int NumProc = Comm.NumProc();
00144 
00145   if(verbose && MyPID==0)
00146     cout << Epetra_Version() << endl << endl;
00147 
00148   if (verbose) cout << "Processor "<<MyPID<<" of "<< NumProc
00149         << " is alive."<<endl;
00150 
00151   // Redefine verbose to only print on PE 0
00152   if(verbose && rank!=0) 
00153     verbose = false;
00154 
00155   int NumMyEquations = 10000;
00156   int NumGlobalEquations = (NumMyEquations * NumProc) + EPETRA_MIN(NumProc,3);
00157   if(MyPID < 3) 
00158     NumMyEquations++;
00159 
00160   // Construct a Map that puts approximately the same Number of equations on each processor
00161 
00162   Epetra_Map Map(NumGlobalEquations, NumMyEquations, 0, Comm);
00163   
00164   // Get update list and number of local equations from newly created Map
00165   vector<int> MyGlobalElements(Map.NumMyElements());
00166   Map.MyGlobalElements(&MyGlobalElements[0]);
00167 
00168   // Create an integer vector NumNz that is used to build the Petra Matrix.
00169   // NumNz[i] is the Number of OFF-DIAGONAL term for the ith global equation on this processor
00170 
00171   vector<int> NumNz(NumMyEquations);
00172 
00173   // We are building a tridiagonal matrix where each row has (-1 2 -1)
00174   // So we need 2 off-diagonal terms (except for the first and last equation)
00175 
00176   for(i = 0; i < NumMyEquations; i++)
00177     if((MyGlobalElements[i] == 0) || (MyGlobalElements[i] == NumGlobalEquations - 1))
00178       NumNz[i] = 1;
00179     else
00180       NumNz[i] = 2;
00181 
00182   // Create a Epetra_Matrix
00183 
00184   Epetra_CrsMatrix A(Copy, Map, &NumNz[0]);
00185   EPETRA_TEST_ERR(A.IndicesAreGlobal(),ierr);
00186   EPETRA_TEST_ERR(A.IndicesAreLocal(),ierr);
00187   
00188   // Add  rows one-at-a-time
00189   // Need some vectors to help
00190   // Off diagonal Values will always be -1
00191 
00192 
00193   vector<double> Values(2);
00194   Values[0] = -1.0; 
00195   Values[1] = -1.0;
00196   vector<int> Indices(2);
00197   double two = 2.0;
00198   int NumEntries;
00199 
00200   forierr = 0;
00201   for(i = 0; i < NumMyEquations; i++) {
00202     if(MyGlobalElements[i] == 0) {
00203       Indices[0] = 1;
00204       NumEntries = 1;
00205     }
00206     else if (MyGlobalElements[i] == NumGlobalEquations-1) {
00207       Indices[0] = NumGlobalEquations-2;
00208       NumEntries = 1;
00209     }
00210     else {
00211       Indices[0] = MyGlobalElements[i]-1;
00212       Indices[1] = MyGlobalElements[i]+1;
00213       NumEntries = 2;
00214     }
00215     forierr += !(A.InsertGlobalValues(MyGlobalElements[i], NumEntries, &Values[0], &Indices[0])==0);
00216     forierr += !(A.InsertGlobalValues(MyGlobalElements[i], 1, &two, &MyGlobalElements[i])>0); // Put in the diagonal entry
00217   }
00218   EPETRA_TEST_ERR(forierr,ierr);
00219 
00220   // Finish up
00221   A.FillComplete();
00222   A.OptimizeStorage();
00223 
00224   HYPRE_IJMatrix Matrix;
00225   int ilower = Map.MinMyGID();
00226   int iupper = Map.MaxMyGID();
00227 
00228   //printf("Proc[%d], ilower = %d, iupper = %d.\n", MyPID, ilower, iupper);
00229   HYPRE_IJMatrixCreate(MPI_COMM_WORLD, ilower, iupper, ilower, iupper, &Matrix);
00230   HYPRE_IJMatrixSetObjectType(Matrix, HYPRE_PARCSR);
00231   HYPRE_IJMatrixInitialize(Matrix);
00232   
00233   for(i = 0; i < A.NumMyRows(); i++){
00234     int numElements;
00235     A.NumMyRowEntries(i, numElements);
00236     vector<int> my_indices; my_indices.resize(numElements);
00237     vector<double> my_values; my_values.resize(numElements);
00238     int numEntries;
00239     A.ExtractMyRowCopy(i, numElements, numEntries, &my_values[0], &my_indices[0]);
00240     for(int j = 0; j < numEntries; j++) {
00241       my_indices[j] = A.GCID(my_indices[j]);
00242     }
00243     int GlobalRow[1];
00244     GlobalRow[0] = A.GRID(i);
00245     HYPRE_IJMatrixSetValues(Matrix, 1, &numEntries, GlobalRow, &my_indices[0], &my_values[0]);
00246   } 
00247   HYPRE_IJMatrixAssemble(Matrix);
00248 
00249   EpetraExt_HypreIJMatrix JadA(Matrix);
00250  
00251   JadA.SetMaps(JadA.RowMatrixRowMap(), A.RowMatrixColMap());
00252   // Create vectors for Power method
00253 
00254   Epetra_Vector q(Map);
00255   Epetra_Vector z(Map); z.Random();
00256   Epetra_Vector resid(Map);
00257 
00258   Epetra_Flops flopcounter;
00259   A.SetFlopCounter(flopcounter);
00260   q.SetFlopCounter(A);
00261   z.SetFlopCounter(A);
00262   resid.SetFlopCounter(A);
00263   JadA.SetFlopCounter(A);
00264 
00265   if (verbose) cout << "=======================================" << endl
00266         << "Testing Jad using CrsMatrix as input..." << endl
00267         << "=======================================" << endl;
00268 
00269   A.ResetFlops();
00270   powerMethodTests(A, JadA, Map, q, z, resid, verbose);
00271 
00272 #ifdef EPETRA_MPI
00273   MPI_Finalize() ;
00274 #endif
00275 
00276 return ierr ;
00277 }
00278 
00279 int powerMethodTests(Epetra_RowMatrix & A, Epetra_RowMatrix & JadA, Epetra_Map & Map, 
00280          Epetra_Vector & q, Epetra_Vector & z, Epetra_Vector & resid, bool verbose) {
00281 
00282   // variable needed for iteration
00283   double lambda = 0.0;
00284   // int niters = 10000;
00285   int niters = 300;
00286   double tolerance = 1.0e-2;
00287   int ierr = 0;
00288 
00290   
00291   // Iterate
00292 
00293   Epetra_Time timer(Map.Comm());
00294   
00295   double startTime = timer.ElapsedTime();
00296   EPETRA_TEST_ERR(power_method(false, A, q, z, resid, &lambda, niters, tolerance, verbose),ierr);
00297   double elapsed_time = timer.ElapsedTime() - startTime;
00298   double total_flops = q.Flops();
00299   double MFLOPs = total_flops/elapsed_time/1000000.0;
00300   double lambdaref = lambda;
00301   double flopsref = total_flops;
00302 
00303   if (verbose) 
00304     cout << "\n\nTotal MFLOPs for reference first solve = " << MFLOPs << endl
00305       <<     "Total FLOPS                            = " <<total_flops <<endl<<endl;
00306 
00307   lambda = 0.0;
00308   startTime = timer.ElapsedTime();
00309   
00310   EPETRA_TEST_ERR(power_method(false, JadA, q, z, resid, &lambda, niters, tolerance, verbose),ierr);
00311   elapsed_time = timer.ElapsedTime() - startTime;
00312   total_flops = q.Flops();
00313   MFLOPs = total_flops/elapsed_time/1000000.0;
00314 
00315   if (verbose) 
00316     cout << "\n\nTotal MFLOPs for candidate first solve = " << MFLOPs << endl
00317       <<     "Total FLOPS                            = " <<total_flops <<endl<<endl;
00318 
00319   EPETRA_TEST_ERR(checkValues(lambda,lambdaref," No-transpose Power Method result", verbose),ierr);
00320   //EPETRA_TEST_ERR(checkValues(total_flops,flopsref," No-transpose Power Method flop count", verbose),ierr);
00321 
00323   
00324   // Solve transpose problem
00325 
00326   if (verbose) cout << "\n\nUsing transpose of matrix and solving again (should give same result).\n\n"
00327         << endl;
00328   // Iterate
00329   lambda = 0.0;
00330   startTime = timer.ElapsedTime();
00331   EPETRA_TEST_ERR(power_method(true, A, q, z, resid, &lambda, niters, tolerance, verbose),ierr);
00332   elapsed_time = timer.ElapsedTime() - startTime;
00333   total_flops = q.Flops();
00334   MFLOPs = total_flops/elapsed_time/1000000.0;
00335   lambdaref = lambda;
00336   flopsref = total_flops;
00337 
00338   if (verbose) 
00339    cout << "\n\nTotal MFLOPs for reference transpose solve = " << MFLOPs << endl
00340      <<     "Total FLOPS                                = " <<total_flops <<endl<<endl;
00341 
00342   lambda = 0.0;
00343   startTime = timer.ElapsedTime();
00344   EPETRA_TEST_ERR(power_method(true, JadA, q, z, resid, &lambda, niters, tolerance, verbose),ierr);
00345   elapsed_time = timer.ElapsedTime() - startTime;
00346   total_flops = q.Flops();
00347   MFLOPs = total_flops/elapsed_time/1000000.0;
00348 
00349   if (verbose) 
00350     cout << "\n\nTotal MFLOPs for candidate transpose solve = " << MFLOPs << endl
00351       <<     "Total FLOPS                                = " <<total_flops <<endl<<endl;
00352 
00353   EPETRA_TEST_ERR(checkValues(lambda,lambdaref,"Transpose Power Method result", verbose),ierr);
00354   //EPETRA_TEST_ERR(checkValues(total_flops,flopsref,"Transpose Power Method flop count", verbose),ierr);
00355 
00356   EPETRA_TEST_ERR(check(A, JadA, verbose),ierr);
00357 
00358   return(0);
00359 }
00360 int power_method(bool TransA, Epetra_RowMatrix& A, Epetra_Vector& q, Epetra_Vector& z0, 
00361      Epetra_Vector& resid, double* lambda, int niters, double tolerance, bool verbose) 
00362 {  
00363   
00364   // Fill z with random Numbers
00365   Epetra_Vector z(z0);
00366   
00367   // variable needed for iteration
00368   double normz, residual;
00369 
00370   int ierr = 1;
00371   
00372   for(int iter = 0; iter < niters; iter++) {
00373     z.Norm2(&normz); // Compute 2-norm of z
00374     q.Scale(1.0/normz, z);
00375     A.Multiply(TransA, q, z); // Compute z = A*q // SEGFAULT HAPPENS HERE
00376     q.Dot(z, lambda); // Approximate maximum eigenvaluE
00377     if(iter%100==0 || iter+1==niters) {
00378       resid.Update(1.0, z, -(*lambda), q, 0.0); // Compute A*q - lambda*q
00379       resid.Norm2(&residual);
00380       if(verbose) cout << "Iter = " << iter << "  Lambda = " << *lambda 
00381                        << "  Residual of A*q - lambda*q = " << residual << endl;
00382     }
00383     if(residual < tolerance) {
00384       ierr = 0;
00385       break;
00386     }
00387   }
00388   return(ierr);
00389 }
00390 
00391 int check(Epetra_RowMatrix& A, Epetra_RowMatrix & B, bool verbose)  {  
00392 
00393   int ierr = 0;
00394   EPETRA_TEST_ERR(!A.Comm().NumProc()==B.Comm().NumProc(),ierr);
00395   EPETRA_TEST_ERR(!A.Comm().MyPID()==B.Comm().MyPID(),ierr);
00396   EPETRA_TEST_ERR(!A.Filled()==B.Filled(),ierr);
00397   EPETRA_TEST_ERR(!A.HasNormInf()==B.HasNormInf(),ierr);
00398   EPETRA_TEST_ERR(!A.LowerTriangular()==B.LowerTriangular(),ierr);
00399   EPETRA_TEST_ERR(!A.Map().SameAs(B.Map()),ierr);
00400   EPETRA_TEST_ERR(!A.MaxNumEntries()==B.MaxNumEntries(),ierr);
00401   EPETRA_TEST_ERR(!A.NumGlobalCols()==B.NumGlobalCols(),ierr);
00402   EPETRA_TEST_ERR(!A.NumGlobalDiagonals()==B.NumGlobalDiagonals(),ierr);
00403   EPETRA_TEST_ERR(!A.NumGlobalNonzeros()==B.NumGlobalNonzeros(),ierr);
00404   EPETRA_TEST_ERR(!A.NumGlobalRows()==B.NumGlobalRows(),ierr);
00405   EPETRA_TEST_ERR(!A.NumMyCols()==B.NumMyCols(),ierr);
00406   EPETRA_TEST_ERR(!A.NumMyDiagonals()==B.NumMyDiagonals(),ierr);
00407   EPETRA_TEST_ERR(!A.NumMyNonzeros()==B.NumMyNonzeros(),ierr);
00408   for (int i=0; i<A.NumMyRows(); i++) {
00409     int nA, nB;
00410     A.NumMyRowEntries(i,nA); B.NumMyRowEntries(i,nB);
00411     EPETRA_TEST_ERR(!nA==nB,ierr);
00412   }
00413   EPETRA_TEST_ERR(!A.NumMyRows()==B.NumMyRows(),ierr);
00414   EPETRA_TEST_ERR(!A.OperatorDomainMap().SameAs(B.OperatorDomainMap()),ierr);
00415   EPETRA_TEST_ERR(!A.OperatorRangeMap().SameAs(B.OperatorRangeMap()),ierr);
00416   
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