EpetraExt Package Browser (Single Doxygen Collection) Development
test/Permutation/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 //Permutation Test routine
00043 #include <Epetra_ConfigDefs.h>
00044 #include "EpetraExt_Version.h"
00045 
00046 #ifdef EPETRA_MPI
00047 #include "Epetra_MpiComm.h"
00048 #include <mpi.h>
00049 #endif
00050 
00051 #include "Epetra_SerialComm.h"
00052 #include "Epetra_Time.h"
00053 #include "Epetra_Map.h"
00054 #include "Epetra_CrsGraph.h"
00055 #include "Epetra_CrsMatrix.h"
00056 #include "Epetra_Vector.h"
00057 #include "EpetraExt_Permutation.h"
00058 #include "../epetra_test_err.h"
00059 
00060 int check_rowpermute_crsmatrix_local_diagonal(Epetra_Comm& Comm, bool verbose);
00061 int check_rowpermute_crsmatrix_global_diagonal(Epetra_Comm& Comm, bool verbose);
00062 int check_rowpermute_crsgraph_local_diagonal(Epetra_Comm& Comm, bool verbose);
00063 int check_colpermute_crsgraph(Epetra_Comm& Comm, bool verbose);
00064 int check_colpermute_crsmatrix(Epetra_Comm& Comm, bool verbose);
00065 int check_rowpermute_multivector_local(Epetra_Comm& Comm, bool verbose);
00066 
00067 int main(int argc, char *argv[]) {
00068 
00069   int returnierr=0;
00070 
00071   bool verbose = false;
00072 
00073 #ifdef EPETRA_MPI
00074 
00075   // Initialize MPI
00076 
00077   MPI_Init(&argc,&argv);
00078   Epetra_MpiComm Comm(MPI_COMM_WORLD);
00079 
00080 #else
00081   Epetra_SerialComm Comm;
00082 #endif
00083 
00084   // Check if we should print results to standard out
00085   if (argc>1) {
00086     if (argv[1][0]=='-' && argv[1][1]=='v') verbose = true;
00087   }
00088 
00089   //Make sure the value of verbose is consistent across processors.
00090   int verbose_int = verbose ? 1 : 0;
00091   Comm.Broadcast(&verbose_int, 1, 0);
00092   verbose = verbose_int==1 ? true : false;
00093 
00094   if (!verbose) {
00095     Comm.SetTracebackMode(0); // This should shut down error traceback reporting
00096   }
00097 
00098   if (verbose && Comm.MyPID()==0)
00099     cout << EpetraExt::EpetraExt_Version() << endl << endl;
00100 
00101   EPETRA_CHK_ERR( check_rowpermute_crsmatrix_local_diagonal( Comm, verbose ) );
00102 
00103   EPETRA_CHK_ERR( check_rowpermute_crsmatrix_global_diagonal( Comm, verbose) );
00104 
00105   EPETRA_CHK_ERR( check_rowpermute_crsgraph_local_diagonal( Comm, verbose) );
00106 
00107   EPETRA_CHK_ERR( check_colpermute_crsgraph( Comm, verbose) );
00108 
00109   EPETRA_CHK_ERR( check_colpermute_crsmatrix( Comm, verbose) );
00110 
00111   EPETRA_CHK_ERR( check_rowpermute_multivector_local( Comm, verbose) );
00112 
00113 
00114 #ifdef EPETRA_MPI
00115   MPI_Finalize();
00116 #endif
00117 
00118   return returnierr;
00119 }
00120 
00121 //------------------------------------------------------------------------------
00122 int check_rowpermute_crsmatrix_local_diagonal(Epetra_Comm& Comm,
00123              bool verbose)
00124 {
00125   int MyPID = Comm.MyPID();
00126   int NumProc = Comm.NumProc();
00127 
00128   Comm.Barrier();
00129   bool verbose1 = verbose;
00130 
00131   if (verbose) verbose = (MyPID==0);
00132 
00133   if (verbose) {
00134     cerr << "================check_rowpermute_crsmatrix_local_diagonal=========="
00135    <<endl;
00136   }
00137 
00138   int NumMyElements = 5;
00139   int NumGlobalElements = NumMyElements*NumProc;
00140  
00141   Epetra_Map Map(NumGlobalElements, NumMyElements, 0, Comm);
00142 
00143   int* p = new int[NumMyElements];
00144   int firstGlobalRow = MyPID*NumMyElements;
00145 
00146   //Set up a permutation that will reverse the order of all LOCAL rows. (i.e.,
00147   //this test won't cause any inter-processor data movement.)
00148 
00149   if (verbose) {
00150     cout << "Permutation P:"<<endl;
00151   }
00152 
00153   int i;
00154 
00155   for(i=0; i<NumMyElements; ++i) {
00156     p[i] = firstGlobalRow+NumMyElements-1-i;
00157     if (verbose1) {
00158       cout << "p["<<firstGlobalRow+i<<"]: "<<p[i]<<endl;
00159     }
00160   }
00161 
00162   Epetra_CrsMatrix A(Copy, Map, 1);
00163 
00164   int col;
00165   double val;
00166 
00167   //set up a diagonal matrix A. It's diagonal because that's the easiest
00168   //to fill and to examine output before and after permutation...
00169 
00170   for(i=0; i<NumMyElements; ++i) {
00171     int row = firstGlobalRow+i;
00172     val = 1.0*row;
00173     col = row;
00174 
00175     A.InsertGlobalValues(row, 1, &val, &col);
00176   }
00177  
00178   A.FillComplete();
00179 
00180   if (verbose1) {
00181     cout << "********** matrix A: **************"<<endl;
00182     cout << A << endl;
00183   }
00184 
00185   EpetraExt::Permutation<Epetra_CrsMatrix> P(Copy, Map, p);
00186 
00187   Epetra_CrsMatrix& B = P(A);
00188 
00189   if (verbose1) {
00190     cout <<"************ permuted matrix B: ***************"<<endl;
00191     cout << B << endl;
00192   }
00193 
00194   return(0);
00195 }
00196 
00197 //------------------------------------------------------------------------------
00198 int check_rowpermute_crsgraph_local_diagonal(Epetra_Comm& Comm,
00199             bool verbose)
00200 {
00201   int MyPID = Comm.MyPID();
00202   int NumProc = Comm.NumProc();
00203 
00204   Comm.Barrier();
00205   bool verbose1 = verbose;
00206 
00207   if (verbose) verbose = (MyPID==0);
00208 
00209   if (verbose) {
00210     cerr << "================check_rowpermute_crsgraph_local_diagonal=========="
00211    <<endl;
00212   }
00213 
00214   int NumMyElements = 5;
00215   int NumGlobalElements = NumMyElements*NumProc;
00216  
00217   Epetra_Map Map(NumGlobalElements, NumMyElements, 0, Comm);
00218 
00219   int* p = new int[NumMyElements];
00220   int firstGlobalRow = MyPID*NumMyElements;
00221 
00222   //Set up a permutation that will reverse the order of all LOCAL rows. (i.e.,
00223   //this test won't cause any inter-processor data movement.)
00224 
00225   if (verbose) {
00226     cout << "Permutation P:"<<endl;
00227   }
00228 
00229   int i;
00230 
00231   for(i=0; i<NumMyElements; ++i) {
00232     p[i] = firstGlobalRow+NumMyElements-1-i;
00233     if (verbose1) {
00234       cout << "p["<<firstGlobalRow+i<<"]: "<<p[i]<<endl;
00235     }
00236   }
00237 
00238   Epetra_CrsGraph Agrph(Copy, Map, 1);
00239 
00240   int col;
00241 
00242   //set up a diagonal graph. It's diagonal because that's the easiest
00243   //to fill and to examine output before and after permutation...
00244 
00245   for(i=0; i<NumMyElements; ++i) {
00246     int row = firstGlobalRow+i;
00247     col = row;
00248 
00249     Agrph.InsertGlobalIndices(row, 1, &col);
00250   }
00251  
00252   Agrph.FillComplete();
00253 
00254   if (verbose1) {
00255     cout << "*************** graph Agrph: ********************"<<endl;
00256     cout << Agrph << endl;
00257   }
00258 
00259   EpetraExt::Permutation<Epetra_CrsGraph> P(Copy, Map, p);
00260 
00261   Epetra_CrsGraph& Bgrph = P(Agrph);
00262 
00263   if (verbose1) {
00264     cout <<"************* permuted graph Bgrph: ****************"<<endl;
00265     cout << Bgrph << endl;
00266   }
00267 
00268   return(0);
00269 }
00270 
00271 //------------------------------------------------------------------------------
00272 int check_colpermute_crsgraph(Epetra_Comm& Comm,
00273             bool verbose)
00274 {
00275   int MyPID = Comm.MyPID();
00276   int NumProc = Comm.NumProc();
00277 
00278   Comm.Barrier();
00279   bool verbose1 = verbose;
00280 
00281   if (verbose) verbose = (MyPID==0);
00282 
00283   if (verbose) {
00284     cerr << "================check_colpermute_crsgraph=========="
00285    <<endl;
00286   }
00287 
00288   int NumMyElements = 5;
00289   int NumGlobalElements = NumMyElements*NumProc;
00290  
00291   Epetra_Map Map(NumGlobalElements, NumMyElements, 0, Comm);
00292 
00293   int* p = new int[NumMyElements];
00294   int firstGlobalRow = MyPID*NumMyElements;
00295 
00296   if (verbose) {
00297     cout << "Permutation P:"<<endl;
00298   }
00299 
00300   int i;
00301 
00302   for(i=0; i<NumMyElements; ++i) {
00303     int row = firstGlobalRow+i;
00304     p[i] = NumGlobalElements - row - 1;
00305     if (verbose1) {
00306       cout << "p["<<firstGlobalRow+i<<"]: "<<p[i]<<endl;
00307     }
00308   }
00309 
00310   Epetra_CrsGraph Agrph(Copy, Map, 1);
00311 
00312   int col;
00313 
00314   //set up a tri-diagonal graph.
00315 
00316   for(i=0; i<NumMyElements; ++i) {
00317     int row = firstGlobalRow+i;
00318     col = NumGlobalElements - row - 1;
00319 
00320     Agrph.InsertGlobalIndices(row, 1, &col);
00321 
00322     if (col > 0) {
00323       int colm1 = col-1;
00324       Agrph.InsertGlobalIndices(row, 1, &colm1);
00325     }
00326 
00327     if (col < NumGlobalElements-1) {
00328       int colp1 = col+1;
00329       Agrph.InsertGlobalIndices(row, 1, &colp1);
00330     }
00331   }
00332  
00333   Agrph.FillComplete();
00334 
00335   if (verbose1) {
00336     cout << "*************** graph Agrph: ********************"<<endl;
00337     cout << Agrph << endl;
00338   }
00339 
00340   EpetraExt::Permutation<Epetra_CrsGraph> P(Copy, Map, p);
00341 
00342   bool column_permutation = true;
00343   Epetra_CrsGraph& Bgrph = P(Agrph, column_permutation);
00344 
00345   if (verbose1) {
00346     cout <<"************* column-permuted graph Bgrph: ****************"<<endl;
00347     cout << Bgrph << endl;
00348   }
00349 
00350   delete [] p;
00351 
00352   return(0);
00353 }
00354 
00355 //-------------------------------------------------------------------------------
00356 int check_rowpermute_crsmatrix_global_diagonal(Epetra_Comm& Comm,
00357       bool verbose)
00358 {
00359   int MyPID = Comm.MyPID();
00360   int NumProc = Comm.NumProc();
00361 
00362   Comm.Barrier();
00363   bool verbose1 = verbose;
00364 
00365   if (verbose) verbose = (MyPID==0);
00366 
00367   if (verbose) {
00368     cerr << "================check_rowpermute_crsmatrix_global_diagonal=========="
00369    <<endl;
00370   }
00371 
00372   int NumMyElements = 5;
00373   int NumGlobalElements = NumMyElements*NumProc;
00374  
00375   Epetra_Map Map(NumGlobalElements, NumMyElements, 0, Comm);
00376 
00377   int* p = new int[NumMyElements];
00378   int firstGlobalRow = MyPID*NumMyElements;
00379 
00380   //Now set up a permutation that will GLOBALLY reverse the order of all rows.
00381   //(i.e., if there are multiple processors, there will be inter-processor
00382   //data movement as rows are migrated.)
00383 
00384   int i;
00385 
00386   Epetra_CrsMatrix A(Copy, Map, 1);
00387 
00388   int col;
00389   double val;
00390 
00391   //set up a diagonal matrix A. It's diagonal because that's the easiest
00392   //to fill and to examine output before and after permutation...
00393 
00394   for(i=0; i<NumMyElements; ++i) {
00395     int row = firstGlobalRow+i;
00396     val = 1.0*row;
00397     col = row;
00398 
00399     A.InsertGlobalValues(row, 1, &val, &col);
00400   }
00401  
00402   A.FillComplete();
00403 
00404   if (verbose1) {
00405     cout << "******************* matrix A: ****************************"<<endl;
00406     cout << A << endl;
00407   }
00408 
00409   if (verbose) {
00410     cout << "Permutation P:"<<endl;
00411   }
00412 
00413   for(i=0; i<NumMyElements; ++i) {
00414     int globalrow = NumGlobalElements-(firstGlobalRow+i)-1;
00415     p[i] = globalrow;
00416     if (verbose1) {
00417       cout << "p["<<firstGlobalRow+i<<"]: "<<p[i]<<endl;
00418     }
00419   }
00420 
00421   EpetraExt::Permutation<Epetra_CrsMatrix> Pglobal(Copy, Map, p);
00422 
00423   Epetra_CrsMatrix& Bglobal = Pglobal(A);
00424 
00425   if (verbose1) {
00426     cout << "******************* permuted matrix Bglobal: *******************" <<endl;
00427     cout << Bglobal << endl;
00428   }
00429 
00430   return(0);
00431 }
00432 
00433 //-------------------------------------------------------------------------------
00434 int check_colpermute_crsmatrix(Epetra_Comm& Comm,
00435              bool verbose)
00436 {
00437   int MyPID = Comm.MyPID();
00438   int NumProc = Comm.NumProc();
00439 
00440   Comm.Barrier();
00441   bool verbose1 = verbose;
00442 
00443   if (verbose) verbose = (MyPID==0);
00444 
00445   if (verbose) {
00446     cerr << "================check_colpermute_crsmatrix=========="
00447    <<endl;
00448   }
00449 
00450   int NumMyElements = 5;
00451   int NumGlobalElements = NumMyElements*NumProc;
00452  
00453   Epetra_Map Map(NumGlobalElements, NumMyElements, 0, Comm);
00454 
00455   int* p = new int[NumMyElements];
00456   int firstGlobalRow = MyPID*NumMyElements;
00457 
00458   if (verbose) {
00459     cout << "Permutation P:"<<endl;
00460   }
00461 
00462   int i;
00463 
00464   for(i=0; i<NumMyElements; ++i) {
00465     int row = firstGlobalRow+i;
00466     p[i] = NumGlobalElements - row - 1;
00467     if (verbose1) {
00468       cout << "p["<<firstGlobalRow+i<<"]: "<<p[i]<<endl;
00469     }
00470   }
00471 
00472   Epetra_CrsMatrix A(Copy, Map, 1);
00473 
00474   int col;
00475   double val;
00476 
00477   //set up a tri-diagonal graph.
00478 
00479   for(i=0; i<NumMyElements; ++i) {
00480     int row = firstGlobalRow+i;
00481     col = NumGlobalElements - row - 1;
00482     val = 1.0*col;
00483 
00484     A.InsertGlobalValues(row, 1, &val, &col);
00485 
00486     if (col > 0) {
00487       int colm1 = col-1;
00488       val = 1.0*colm1;
00489       A.InsertGlobalValues(row, 1, &val, &colm1);
00490     }
00491 
00492     if (col < NumGlobalElements-1) {
00493       int colp1 = col+1;
00494       val = 1.0*colp1;
00495       A.InsertGlobalValues(row, 1, &val, &colp1);
00496     }
00497   }
00498  
00499   A.FillComplete();
00500 
00501   if (verbose1) {
00502     cout << "*************** matrix A: ********************"<<endl;
00503     cout << A << endl;
00504   }
00505 
00506   EpetraExt::Permutation<Epetra_CrsMatrix> P(Copy, Map, p);
00507 
00508   bool column_permutation = true;
00509   Epetra_CrsMatrix& B = P(A, column_permutation);
00510 
00511   if (verbose1) {
00512     cout <<"************* column-permuted matrix B: ****************"<<endl;
00513     cout << B << endl;
00514   }
00515 
00516   delete [] p;
00517 
00518   return(0);
00519 }
00520 
00521 //------------------------------------------------------------------------------
00522 int check_rowpermute_multivector_local(Epetra_Comm& Comm,
00523                bool verbose)
00524 {
00525   int MyPID = Comm.MyPID();
00526   int NumProc = Comm.NumProc();
00527 
00528   Comm.Barrier();
00529   bool verbose1 = verbose;
00530 
00531   if (verbose) verbose = (MyPID==0);
00532 
00533   if (verbose) {
00534     cerr << "================check_rowpermute_multivector_local=========="
00535    <<endl;
00536   }
00537 
00538   int NumMyElements = 5;
00539   int NumGlobalElements = NumMyElements*NumProc;
00540  
00541   Epetra_Map Map(NumGlobalElements, NumMyElements, 0, Comm);
00542 
00543   int* p = new int[NumMyElements];
00544   int firstGlobalRow = MyPID*NumMyElements;
00545 
00546   //Set up a permutation that will reverse the order of all LOCAL rows. (i.e.,
00547   //this test won't cause any inter-processor data movement.)
00548 
00549   if (verbose) {
00550     cout << "Permutation P:"<<endl;
00551   }
00552 
00553   int i;
00554 
00555   for(i=0; i<NumMyElements; ++i) {
00556     p[i] = firstGlobalRow+NumMyElements-1-i;
00557     if (verbose1) {
00558       cout << "p["<<firstGlobalRow+i<<"]: "<<p[i]<<endl;
00559     }
00560   }
00561 
00562   Epetra_MultiVector v(Map, 3);
00563 
00564   double* v0 = v[0];
00565   double* v1 = v[1];
00566   double* v2 = v[2];
00567 
00568   for(i=0; i<NumMyElements; ++i) {
00569     v0[i] = 1.0*(firstGlobalRow+i) + 0.1;
00570     v1[i] = 1.0*(firstGlobalRow+i) + 0.2;
00571     v2[i] = 1.0*(firstGlobalRow+i) + 0.3;
00572   }
00573  
00574   if (verbose1) {
00575     cout << "*************** MultiVector v: ********************"<<endl;
00576     cout << v << endl;
00577   }
00578 
00579   EpetraExt::Permutation<Epetra_MultiVector> P(Copy, Map, p);
00580 
00581   Epetra_MultiVector& Pv = P(v);
00582 
00583   if (verbose1) {
00584     cout <<"************* permuted MultiVector Pv: ****************"<<endl;
00585     cout << Pv << endl;
00586   }
00587 
00588   return(0);
00589 }
00590 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines