Epetra Package Browser (Single Doxygen Collection) Development
test/CrsRectMatrix/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_LocalMap.h"
00044 #include "Epetra_Map.h"
00045 #include "Epetra_Time.h"
00046 #include "Epetra_Vector.h"
00047 #include "Epetra_CrsMatrix.h"
00048 #include "Epetra_Flops.h"
00049 #include "Epetra_Export.h"
00050 #include "Epetra_Import.h"
00051 #include "Epetra_SerialDenseVector.h"
00052 #include "Epetra_IntSerialDenseVector.h"
00053 #ifdef EPETRA_MPI
00054 #include "Epetra_MpiComm.h"
00055 #include "mpi.h"
00056 #else
00057 #include "Epetra_SerialComm.h"
00058 #endif
00059 #include "../epetra_test_err.h"
00060 #include "Epetra_Version.h"
00061  
00062 int main(int argc, char *argv[])
00063 {
00064   int ierr = 0, i, j, forierr = 0;
00065   int ntrials = 1;
00066 
00067 #ifdef EPETRA_MPI
00068 
00069   // Initialize MPI
00070 
00071   MPI_Init(&argc,&argv);
00072   Epetra_MpiComm Comm( MPI_COMM_WORLD );
00073 
00074 #else
00075 
00076   Epetra_SerialComm Comm;
00077 
00078 #endif
00079 
00080   bool verbose = false;
00081 
00082   // Check if we should print results to standard out
00083   if (argc>1) if (argv[1][0]=='-' && argv[1][1]=='v') verbose = true;
00084 
00085 
00086 
00087   // char tmp; if (Comm.MyPID()==0) { cout << "Press any key to continue..."<< endl;  cin >> tmp;} Comm.Barrier();
00088 
00089   Comm.SetTracebackMode(0); // This should shut down any error traceback reporting
00090   int MyPID = Comm.MyPID();
00091   int NumProc = Comm.NumProc();
00092 
00093   if (verbose && Comm.MyPID()==0)
00094     cout << Epetra_Version() << endl << endl;
00095 
00096   if (verbose) cout << "Processor "<<MyPID<<" of "<< NumProc
00097               << " is alive."<<endl;
00098 
00099   bool verbose1 = verbose;
00100 
00101   // Redefine verbose to only print on PE 0
00102   if (verbose && Comm.MyPID()!=0) verbose = false;
00103 
00104   int NumMyEquations = 10000;
00105 
00106   int NumGlobalEquations = NumMyEquations*NumProc;
00107   int NumGlobalVariables = 2 * NumGlobalEquations+1;
00108 
00109   // Construct a Map that puts approximately the same Number of equations on each processor
00110 
00111   Epetra_Map RowMap(NumGlobalEquations, 0, Comm);
00112   Epetra_Map XMap(NumGlobalVariables, 0, Comm);
00113   Epetra_Map& YMap = RowMap;
00114   
00115   // Get update list and number of local equations from newly created Map
00116   int * MyGlobalElements = new int[RowMap.NumMyElements()];
00117   RowMap.MyGlobalElements(MyGlobalElements);
00118 
00119   // Get update list and number of local equations from newly created XMap
00120   int * XGlobalElements = new int[XMap.NumMyElements()];
00121   XMap.MyGlobalElements(XGlobalElements);
00122 
00123   // Get update list and number of local variables from newly created YMap
00124   int * YGlobalElements = new int[YMap.NumMyElements()];
00125   YMap.MyGlobalElements(YGlobalElements);
00126 
00127   // We need vectors to compute:
00128   // X = A^T*Y
00129   // AATY = A*A^T*Y = A*X
00130   //  and 
00131   // BY = B*Y
00132 
00133   Epetra_Vector Y(YMap);
00134   Epetra_Vector X(XMap);
00135   Epetra_Vector AATY(YMap);
00136   Epetra_Vector BY(YMap);
00137 
00138 
00139   // Fill Y Vector
00140   Y.Random();
00141   //Y.PutScalar(1.0);
00142 
00143   // To create A^T explicitly we need an assembly map that is two elements longer than
00144   // the XMap, because each processor will be making contributions to two rows beyond what
00145   // it will own.
00146   int ATAssemblyNumMyElements = 2*MyGlobalElements[NumMyEquations-1] + 2 - 2*MyGlobalElements[0] + 1;
00147   int * ATAssemblyGlobalElements = new int[ATAssemblyNumMyElements];
00148 
00149   for (i=0; i<ATAssemblyNumMyElements; i++) ATAssemblyGlobalElements[i] = 2*MyGlobalElements[0] + i;
00150   Epetra_Map ATAssemblyMap(-1, ATAssemblyNumMyElements, ATAssemblyGlobalElements, 0, Comm);
00151 
00152   // Create a Epetra_Matrix with the values of A
00153   // A is a simple 1D weighted average operator that mimics a restriction operator
00154   // that might be found in a multigrid code.
00155   // Also create A^T explicitly
00156 
00157   Epetra_CrsMatrix A(Copy, RowMap, 3);
00158   Epetra_CrsMatrix ATAssembly(Copy, ATAssemblyMap, 2);
00159   Epetra_CrsMatrix AT(Copy, XMap, 2);
00160   
00161   //cout << "ATAssemblyMap = "<< endl<< ATAssemblyMap << endl
00162   //     << "XMap = " << endl << XMap << endl
00163   //     << "RowMap = " << endl << RowMap << endl;
00164   // Add  rows one-at-a-time
00165   // Need some vectors to help
00166   // Off diagonal Values will always be -1
00167 
00168 
00169   double *Values = new double[3];
00170   int *Indices = new int[3];
00171   int NumEntries;
00172   /*
00173   Values[0] = 0.25;
00174   Values[1] = 0.5;
00175   Values[2] = 0.25;
00176   */
00177   Values[0] = 0.5;
00178   Values[1] = 0.25;
00179   Values[2] = 0.25;
00180   forierr = 0;
00181   for (i=0; i<NumMyEquations; i++)
00182     {
00183   /*
00184       Indices[0] = 2*MyGlobalElements[i];
00185       Indices[1] = 2*MyGlobalElements[i]+1;
00186       Indices[2] = 2*MyGlobalElements[i]+2;
00187    */
00188       Indices[0] = 2*MyGlobalElements[i]+1;
00189       Indices[1] = 2*MyGlobalElements[i]+2;
00190       Indices[2] = 2*MyGlobalElements[i];
00191       NumEntries = 3;
00192       forierr += !(A.InsertGlobalValues(MyGlobalElements[i], NumEntries, Values, Indices)==0);
00193       for (j=0; j<3; j++)
00194   forierr += !(ATAssembly.InsertGlobalValues(Indices[j],1, &(Values[j]), &(MyGlobalElements[i]))>=0);
00195     }
00196   EPETRA_TEST_ERR(forierr,ierr);
00197 
00198 
00199   EPETRA_TEST_ERR(!(ATAssembly.FillComplete()==0),ierr);
00200   // Gather AT values from ATAssembly matrix
00201   Epetra_Export Exporter(ATAssemblyMap, XMap);
00202   EPETRA_TEST_ERR(!(AT.Export(ATAssembly, Exporter, Add)==0),ierr);
00203 
00204   // Finish up
00205   EPETRA_TEST_ERR(!(A.FillComplete(XMap, YMap)==0),ierr);
00206   EPETRA_TEST_ERR(!(AT.FillComplete(YMap, XMap)==0),ierr);
00207 
00208 
00209   if (verbose1 && NumGlobalEquations<20) { 
00210     if (verbose) cout << "\n\n Matrix A\n" << endl;
00211     cout << A << endl;
00212     if (verbose) cout << " \n\n Matrix A Transpose\n" << endl;
00213     cout <<  AT << endl;
00214   }
00215 
00216 
00217   // Create a Epetra_Matrix containing B = A*A^T.
00218   // This matrix will be a square tridiagonal matrix.  We will use it to compare the results
00219   // of A*(A^T*X) using two methods: (1) with two calls to Multiply using A^T and then A and
00220   // (2) using B directly.
00221 
00222   Epetra_CrsMatrix B(Copy, RowMap, 3);
00223 
00224   Values[0] = 1.0/16.0;
00225   Values[1] = 3.0/8.0;
00226   Values[2] = 1.0/16.0;
00227   int Valstart;
00228   forierr = 0;
00229   for (i=0; i<NumMyEquations; i++)
00230     {
00231       if (MyGlobalElements[i] == 0) {
00232       Indices[0] = MyGlobalElements[i];
00233       Indices[1] = MyGlobalElements[i]+1;
00234       NumEntries = 2;
00235       Valstart = 1;
00236       }
00237       else {
00238   Indices[0] = MyGlobalElements[i]-1;
00239   Indices[1] = MyGlobalElements[i];
00240   Indices[2] = MyGlobalElements[i]+1;
00241   NumEntries = 3;
00242   Valstart = 0;
00243       }
00244       if (MyGlobalElements[i] == NumGlobalEquations-1) NumEntries--;
00245       forierr += !(B.InsertGlobalValues(MyGlobalElements[i], NumEntries, Values+Valstart, Indices)==0);
00246     }
00247   EPETRA_TEST_ERR(forierr,ierr);
00248 
00249   // Finish up
00250   EPETRA_TEST_ERR(!(B.FillComplete()==0),ierr);
00251   if (verbose && NumGlobalEquations<20) cout << "\n\nMatrix B \n" << endl;
00252   if (verbose1 && NumGlobalEquations<20) cout << B << endl;
00253 
00254 
00255   Epetra_Flops counter;
00256   A.SetFlopCounter(counter);
00257   B.SetFlopCounter(A);
00258   Epetra_Time timer(Comm);
00259   for (i=0; i<ntrials; i++) {
00260     EPETRA_TEST_ERR(!(B.Multiply(false, Y, BY)==0),ierr); // Compute BY = B*Y
00261   }
00262   double elapsed_time = timer.ElapsedTime();
00263   double total_flops = B.Flops();
00264   counter.ResetFlops();
00265   
00266   double MFLOPs = total_flops/elapsed_time/1000000.0;
00267 
00268   if (verbose) cout << "\n\nTotal MFLOPs for B*Y = " << MFLOPs << endl<< endl;
00269   if (verbose && NumGlobalEquations<20) cout << "\n\nVector Z = B*Y \n";
00270   if (verbose1 && NumGlobalEquations<20) cout << BY << endl;
00271  
00272 
00274 
00275   timer.ResetStartTime();
00276   for (i=0; i<ntrials; i++) {
00277     EPETRA_TEST_ERR(!(A.Multiply(true, Y, X)==0),ierr); // Compute X = A^T*Y
00278   }
00279   elapsed_time = timer.ElapsedTime();
00280   total_flops = A.Flops();
00281   counter.ResetFlops();
00282   MFLOPs = total_flops/elapsed_time/1000000.0;
00283 
00284   if (verbose) cout << "\n\nTotal MFLOPs for A^T*Y using A and trans=true = " << MFLOPs << endl<< endl;
00285   if (verbose && NumGlobalEquations<20) cout << "\n\nVector Z = AT*Y \n";
00286   if (verbose1 && NumGlobalEquations<20) cout << X << endl;
00287 
00289 
00290   timer.ResetStartTime();
00291   EPETRA_TEST_ERR(!(A.Multiply(false, X, AATY)==0),ierr); // Compute AATY = A*X
00292   elapsed_time = timer.ElapsedTime();
00293   total_flops = A.Flops();
00294   MFLOPs = total_flops/elapsed_time/1000000.0;
00295   counter.ResetFlops();
00296   Epetra_Vector resid(YMap);
00297   resid.Update(1.0, BY, -1.0, AATY, 0.0);
00298   double residual;
00299   resid.Norm2(&residual);
00300 
00301   if (verbose) cout << "\n\nTotal MFLOPs for A*X using A and trans=false = " << MFLOPs << endl<< endl;
00302   if (verbose) cout << "Residual = " << residual << endl<< endl;
00303   if (verbose && NumGlobalEquations<20) cout << "\n\nVector Z = A*ATY \n";
00304   if (verbose1 && NumGlobalEquations<20) cout << AATY << endl;
00305 
00307 
00308   AT.SetFlopCounter(counter);
00309   timer.ResetStartTime();
00310   for (i=0; i<ntrials; i++) {
00311     EPETRA_TEST_ERR(!(AT.Multiply(false, Y, X)==0),ierr); // Compute X = A^T*Y
00312   }
00313   elapsed_time = timer.ElapsedTime();
00314   total_flops = AT.Flops();
00315   counter.ResetFlops();
00316   MFLOPs = total_flops/elapsed_time/1000000.0;
00317 
00318   if (verbose) cout << "\n\nTotal MFLOPs for A^T*Y using AT and trans=false = " << MFLOPs << endl<< endl;
00319   if (verbose && NumGlobalEquations<20) cout << "\n\nVector Z = AT*Y \n";
00320   if (verbose1 && NumGlobalEquations<20) cout << X << endl;
00321 
00323 
00324   timer.ResetStartTime();
00325   for (i=0; i<ntrials; i++) {
00326     EPETRA_TEST_ERR(!(AT.Multiply(true, X, AATY)==0),ierr); // Compute AATY = A*X
00327   }
00328   elapsed_time = timer.ElapsedTime();
00329   total_flops = AT.Flops();
00330   MFLOPs = total_flops/elapsed_time/1000000.0;
00331   counter.ResetFlops();
00332   resid.Update(1.0, BY, -1.0, AATY, 0.0);
00333   resid.Norm2(&residual);
00334 
00335   if (verbose) cout << "\n\nTotal MFLOPs for A*X using AT and trans=true = " << MFLOPs << endl<< endl;
00336   if (verbose) cout << "Residual = " << residual << endl<< endl;
00337   if (verbose && NumGlobalEquations<20) cout << "\n\nVector Z = A*ATY \n";
00338   if (verbose1 && NumGlobalEquations<20) cout <<AATY << endl;
00339 
00340 
00342   // Now test use of Epetra_LocalMap vectors: First test case of local replicated range vector
00343 
00344   {
00345     Epetra_CrsMatrix AL(Copy, RowMap, 3);
00346     for (i=0; i<NumMyEquations; i++)
00347       {
00348   forierr += !(A.ExtractGlobalRowCopy(MyGlobalElements[i], 3, NumEntries, Values, Indices)==0);
00349   forierr += !(AL.InsertGlobalValues(MyGlobalElements[i], NumEntries, Values, Indices)==0);
00350       }
00351     EPETRA_TEST_ERR(forierr,ierr);
00352 
00353     Epetra_LocalMap YLMap(NumGlobalEquations, 0, Comm);
00354     EPETRA_TEST_ERR(!(AL.FillComplete(XMap, YLMap)==0),ierr);
00355     AL.SetFlopCounter(A);
00356     Epetra_Vector YL(YLMap);
00357     Epetra_Vector ALX(YLMap);
00358   
00359     timer.ResetStartTime();
00360     for (i=0; i<ntrials; i++) {
00361       EPETRA_TEST_ERR(!(A.Multiply(false, X, Y)==0),ierr); // Compute Y= A*X
00362     }
00363     elapsed_time = timer.ElapsedTime();
00364     total_flops = A.Flops();
00365     counter.ResetFlops();
00366     MFLOPs = total_flops/elapsed_time/1000000.0;
00367 
00368 
00369 
00370     if (verbose) cout << "\n\nTotal MFLOPs for Y=A*X using global distributed Y = " << MFLOPs << endl<< endl;
00371     if (verbose && NumGlobalEquations<20) cout << "\n\nVector Y = A*X using distributed Y \n";
00372     if (verbose1 && NumGlobalEquations<20) cout << Y << endl;
00373     if (verbose && NumGlobalEquations<20) cout << "\n\nA using dist Y range map\n";
00374     if (verbose1 && NumGlobalEquations<20) cout << A << endl;
00375 
00376     timer.ResetStartTime();
00377     for (i=0; i<ntrials; i++) {
00378       EPETRA_TEST_ERR(!(AL.Multiply(false, X, ALX)==0),ierr); // Compute YL= A*X
00379     }
00380     elapsed_time = timer.ElapsedTime();
00381     total_flops = AL.Flops();
00382     counter.ResetFlops();
00383     MFLOPs = total_flops/elapsed_time/1000000.0;
00384 
00385     if (verbose) cout << "\n\nTotal MFLOPs for Y=A*X using Local replicated Y = " << MFLOPs << endl<< endl;
00386     if (verbose && NumGlobalEquations<20) cout << "\n\nVector YL = AL*X using local replicated Y \n";
00387     if (verbose1 && NumGlobalEquations<20) cout << ALX << endl;
00388     if (verbose && NumGlobalEquations<20) cout << "\n\nA using local Y range map\n";
00389     if (verbose1 && NumGlobalEquations<20) cout << AL << endl;
00390 
00391     // Now gather Y values from the distributed Y and compare them to the local replicated Y values
00392     Epetra_Import g2limporter(YLMap, YMap); // Import from distributed Y map to local Y map
00393     EPETRA_TEST_ERR(!(YL.Import(Y, g2limporter, Insert)==0),ierr);
00394     if (verbose && NumGlobalEquations<20) cout << "\n\nVector YL = imported from distributed Y \n";
00395     if (verbose1 && NumGlobalEquations<20) cout << YL << endl;
00396     EPETRA_TEST_ERR(!(YL.Update(-1.0, ALX, 1.0)==0),ierr);
00397     EPETRA_TEST_ERR(!(YL.Norm2(&residual)==0),ierr);
00398     if (verbose) cout << "Residual = " << residual << endl<< endl;
00399       
00400 
00401     // 
00402     // Multiply by transpose
00403     //
00404 
00405     timer.ResetStartTime();
00406     for (i=0; i<ntrials; i++) {
00407       EPETRA_TEST_ERR(!(A.Multiply(true, Y, X)==0),ierr); // Compute X = A^TY
00408     }
00409     elapsed_time = timer.ElapsedTime();
00410     total_flops = A.Flops();
00411     counter.ResetFlops();
00412     MFLOPs = total_flops/elapsed_time/1000000.0;
00413 
00414 
00415 
00416     if (verbose) cout << "\n\nTotal MFLOPs for X=A^TY using global distributed Y = " << MFLOPs << endl<< endl;
00417     if (verbose && NumGlobalEquations<20) cout << "\n\nVector X using distributed Y \n";
00418     if (verbose1 && NumGlobalEquations<20) cout << Y << endl;
00419     if (verbose && NumGlobalEquations<20) cout << "\n\nA using dist Y range map\n";
00420     if (verbose1 && NumGlobalEquations<20) cout << A << endl;
00421 
00422     Epetra_Vector X1(XMap);
00423 
00424     timer.ResetStartTime();
00425     for (i=0; i<ntrials; i++) {
00426       EPETRA_TEST_ERR(!(AL.Multiply(true, ALX, X1)==0),ierr); // Compute X1 = AL^T*Y
00427     }
00428     elapsed_time = timer.ElapsedTime();
00429     total_flops = AL.Flops();
00430     counter.ResetFlops();
00431     MFLOPs = total_flops/elapsed_time/1000000.0;
00432 
00433     if (verbose) cout << "\n\nTotal MFLOPs for X1=A^T*Y using Local replicated Y = " << MFLOPs << endl<< endl;
00434     if (verbose && NumGlobalEquations<20) cout << "\n\nVector X1 using local replicated Y \n";
00435     if (verbose1 && NumGlobalEquations<20) cout << X1 << endl;
00436 
00437     EPETRA_TEST_ERR(!(X1.Update(-1.0, X, 1.0)==0),ierr);
00438     EPETRA_TEST_ERR(!(X1.Norm2(&residual)==0),ierr);
00439     if (verbose) cout << "Residual = " << residual << endl<< endl;
00440   }
00441       
00443   // Finally test use of Epetra_LocalMap vectors using local replicated domain vector
00444 
00445   {
00446     Epetra_CrsMatrix AL(Copy, RowMap, 3);
00447     for (i=0; i<NumMyEquations; i++)
00448       {
00449   forierr += !(A.ExtractGlobalRowCopy(MyGlobalElements[i], 3, NumEntries, Values, Indices)==0);
00450   forierr += !(AL.InsertGlobalValues(MyGlobalElements[i], NumEntries, Values, Indices)==0);
00451       }
00452     EPETRA_TEST_ERR(forierr,ierr);
00453 
00454     Epetra_LocalMap XLMap(NumGlobalVariables, 0, Comm);
00455     EPETRA_TEST_ERR(!(AL.FillComplete(XLMap, YMap)==0),ierr);
00456     AL.SetFlopCounter(A);
00457     Epetra_Vector XL(XLMap);
00458     Epetra_Vector ALX(XLMap);
00459   
00460     timer.ResetStartTime();
00461     for (i=0; i<ntrials; i++) {
00462       EPETRA_TEST_ERR(!(A.Multiply(false, X, Y)==0),ierr); // Compute Y= A*X
00463     }
00464     elapsed_time = timer.ElapsedTime();
00465     total_flops = A.Flops();
00466     counter.ResetFlops();
00467     MFLOPs = total_flops/elapsed_time/1000000.0;
00468 
00469 
00470 
00471     if (verbose) cout << "\n\nTotal MFLOPs for Y=A*X using global distributed X = " << MFLOPs << endl<< endl;
00472     if (verbose && NumGlobalEquations<20) cout << "\n\nVector Y = A*X using distributed X \n";
00473     if (verbose1 && NumGlobalEquations<20) cout << Y << endl;
00474     //if (verbose && NumGlobalEquations<20) cout << "\n\nA using dist X range map\n";
00475     //if (verbose1 && NumGlobalEquations<20) cout << A << endl;
00476 
00477     // Now gather X values from the distributed X 
00478     Epetra_Import g2limporter(XLMap, XMap); // Import from distributed X map to local X map
00479     EPETRA_TEST_ERR(!(XL.Import(X, g2limporter, Insert)==0),ierr);
00480     if (verbose && NumGlobalEquations<20) cout << "\n\nVector XL = imported from distributed X \n";
00481     if (verbose1 && NumGlobalEquations<20) cout << XL << endl;
00482     Epetra_Vector Y1(Y);
00483     timer.ResetStartTime();
00484     for (i=0; i<ntrials; i++) {
00485       EPETRA_TEST_ERR(!(AL.Multiply(false, XL, Y1)==0),ierr); // Compute Y1= AL*XL
00486     }
00487     elapsed_time = timer.ElapsedTime();
00488     total_flops = AL.Flops();
00489     counter.ResetFlops();
00490     MFLOPs = total_flops/elapsed_time/1000000.0;
00491 
00492     if (verbose) cout << "\n\nTotal MFLOPs for Y1=A*XL using Local replicated X = " << MFLOPs << endl<< endl;
00493     if (verbose && NumGlobalEquations<20) cout << "\n\nVector Y1 = AL*XL using local replicated X \n";
00494     if (verbose1 && NumGlobalEquations<20) cout << Y1 << endl;
00495     //if (verbose && NumGlobalEquations<20) cout << "\n\nA using local X domain map\n";
00496     //if (verbose1 && NumGlobalEquations<20) cout << AL << endl;
00497 
00498     EPETRA_TEST_ERR(!(Y1.Update(-1.0, Y, 1.0)==0),ierr);
00499     EPETRA_TEST_ERR(!(Y1.Norm2(&residual)==0),ierr);
00500     if (verbose) cout << "Residual = " << residual << endl<< endl;
00501       
00502 
00503     // 
00504     // Multiply by transpose
00505     //
00506 
00507     timer.ResetStartTime();
00508     for (i=0; i<ntrials; i++) {
00509       EPETRA_TEST_ERR(!(A.Multiply(true, Y, X)==0),ierr); // Compute X = A^TY
00510     }
00511     elapsed_time = timer.ElapsedTime();
00512     total_flops = A.Flops();
00513     counter.ResetFlops();
00514     MFLOPs = total_flops/elapsed_time/1000000.0;
00515 
00516 
00517 
00518     if (verbose) cout << "\n\nTotal MFLOPs for X=A^TY using global distributed X = " << MFLOPs << endl<< endl;
00519     if (verbose && NumGlobalEquations<20) cout << "\n\nVector X using distributed X \n";
00520     if (verbose1 && NumGlobalEquations<20) cout << X << endl;
00521     //if (verbose && NumGlobalEquations<20) cout << "\n\nA using dist X domain map\n";
00522     //if (verbose1 && NumGlobalEquations<20) cout << A << endl;
00523 
00524     timer.ResetStartTime();
00525     for (i=0; i<ntrials; i++) {
00526       EPETRA_TEST_ERR(!(AL.Multiply(true, Y, XL)==0),ierr); // Compute XL = AL^T*Y
00527     }
00528     elapsed_time = timer.ElapsedTime();
00529     total_flops = AL.Flops();
00530     counter.ResetFlops();
00531     MFLOPs = total_flops/elapsed_time/1000000.0;
00532 
00533     if (verbose) cout << "\n\nTotal MFLOPs for XL=A^T*Y1 using Local replicated X = " << MFLOPs << endl<< endl;
00534     if (verbose && NumGlobalEquations<20) cout << "\n\nVector XL using local replicated X \n";
00535     if (verbose1 && NumGlobalEquations<20) cout << XL << endl;
00536 
00537     Epetra_Vector XL1(XLMap);
00538     EPETRA_TEST_ERR(!(XL1.Import(X, g2limporter, Insert)==0),ierr);
00539     EPETRA_TEST_ERR(!(XL1.Update(-1.0, XL, 1.0)==0),ierr);
00540     EPETRA_TEST_ERR(!(XL1.Norm2(&residual)==0),ierr);
00541     if (verbose) cout << "Residual = " << residual << endl<< endl;
00542   }     
00543   // Release all objects
00544   delete [] Values;
00545   delete [] Indices;
00546   delete [] MyGlobalElements;
00547   delete [] XGlobalElements;
00548   delete [] YGlobalElements;
00549   delete [] ATAssemblyGlobalElements;
00550 
00551 
00552 #ifdef EPETRA_MPI
00553   MPI_Finalize() ;
00554 #endif
00555 
00556 /* end main
00557 */
00558 return ierr ;
00559 }
00560 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines