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