part_redist.cpp

This example demonstrates repartitioning and redistributing the contents of an Epetra_LinearProblem object, using the Isorropia::Partitioner and Isorropia::Redistributor classes. This program does not use user-specified weights/costs.

00001 //@HEADER
00002 // ************************************************************************
00003 //
00004 //               Isorropia: Partitioning and Load Balancing Package
00005 //                 Copyright (2006) 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 //
00025 // ************************************************************************
00026 //@HEADER
00027 
00028 //--------------------------------------------------------------------
00029 //This file is a self-contained example of creating an Epetra_LinearProblem
00030 //object, and using Isorropia to create a rebalanced copy of it.
00031 //--------------------------------------------------------------------
00032 
00033 //Include Isorropia_Exception.hpp only because the helper functions at
00034 //the bottom of this file (which create the epetra objects) can
00035 //potentially throw exceptions.
00036 #include <Isorropia_Exception.hpp>
00037 
00038 //The Isorropia symbols being demonstrated are declared
00039 //in these headers:
00040 #include <Isorropia_Epetra.hpp>
00041 #include <Isorropia_EpetraRedistributor.hpp>
00042 #include <Isorropia_EpetraPartitioner.hpp>
00043 #include <Isorropia_EpetraCostDescriber.hpp>
00044 
00045 #ifdef HAVE_MPI
00046 #include <mpi.h>
00047 #endif
00048 
00049 #ifdef HAVE_EPETRA
00050 #ifdef HAVE_MPI
00051 #include <Epetra_MpiComm.h>
00052 #else
00053 #include <Epetra_SerialComm.h>
00054 #endif
00055 #include <Epetra_Map.h>
00056 #include <Epetra_Vector.h>
00057 #include <Epetra_CrsMatrix.h>
00058 #include <Epetra_LinearProblem.h>
00059 #endif
00060 
00061 #include "ispatest_lbeval_utils.hpp"
00062 
00063 //Declaration for helper-function that creates epetra objects. This
00064 //function is implemented at the bottom of this file.
00065 #ifdef HAVE_EPETRA
00066 Epetra_LinearProblem* create_epetra_problem(int numProcs,
00067                                             int localProc,
00068                                             int local_n);
00069 #endif
00070 
00071 int main(int argc, char** argv) {
00072 #if defined(HAVE_MPI) && defined(HAVE_EPETRA)
00073 
00074   int p, numProcs = 1;
00075   int localProc = 0;
00076 
00077   //first, set up our MPI environment...
00078   MPI_Init(&argc, &argv);
00079   MPI_Comm_rank(MPI_COMM_WORLD, &localProc);
00080   MPI_Comm_size(MPI_COMM_WORLD, &numProcs);
00081 
00082   int local_n = 600;
00083 
00084   //Create a Epetra_LinearProblem object.
00085 
00086   Epetra_LinearProblem* linprob = 0;
00087   try {
00088     linprob = create_epetra_problem(numProcs, localProc, local_n);
00089   }
00090   catch(std::exception& exc) {
00091     std::cout << "linsys example: create_epetra_problem threw exception '"
00092           << exc.what() << "' on proc " << localProc << std::endl;
00093     MPI_Finalize();
00094     return(-1);
00095   }
00096 
00097   //We'll need a Teuchos::ParameterList object to pass to the
00098   //Isorropia::Epetra::Partitioner class.
00099   Teuchos::ParameterList paramlist;
00100 
00101   // If Zoltan is available, the Zoltan package will be used for
00102   // the partitioning operation. By default, Isorropia selects Zoltan's
00103   // Hypergraph partitioner. If a method other than Hypergraph is
00104   // desired, it can be specified by first creating a parameter sublist
00105   // named "Zoltan", and then setting appropriate Zoltan parameters in
00106   // that sublist. A sublist is created like this:
00107   //      Teuchos::ParameterList& sublist = paramlist.sublist("Zoltan");
00108   //
00109 
00110   // If Zoltan is not available, a simple linear partitioner will be
00111   // used to partition such that the number of nonzeros is equal (or
00112   // close to equal) on each processor.
00113 
00114 
00115   Epetra_RowMatrix* rowmatrix = linprob->GetMatrix();
00116   Teuchos::RCP<const Epetra_RowMatrix> rowmat =
00117     Teuchos::rcp(rowmatrix, false);
00118 
00119 
00120   //Now create the partitioner 
00121 
00122   Teuchos::RCP<Isorropia::Epetra::Partitioner> partitioner =
00123     Teuchos::rcp(new Isorropia::Epetra::Partitioner(rowmat, paramlist));
00124 
00125   //Next create a Redistributor object and use it to create balanced
00126   //copies of the objects in linprob.
00127 
00128   Isorropia::Epetra::Redistributor rd(partitioner);
00129 
00130   Teuchos::RCP<Epetra_CrsMatrix> bal_matrix;
00131   Teuchos::RCP<Epetra_MultiVector> bal_x;
00132   Teuchos::RCP<Epetra_MultiVector> bal_b;
00133 
00134   //Use a try-catch block because Isorropia will throw an exception
00135   //if it encounters an error.
00136 
00137   if (localProc == 0) {
00138     std::cout << " calling Isorropia::Epetra::Redistributor::redistribute..."
00139         << std::endl;
00140   }
00141 
00142   try {
00143     bal_matrix = rd.redistribute(*linprob->GetMatrix());
00144     bal_x = rd.redistribute(*linprob->GetLHS());
00145     bal_b = rd.redistribute(*linprob->GetRHS());
00146   }
00147   catch(std::exception& exc) {
00148     std::cout << "linsys example: Isorropia::Epetra::Redistributor threw "
00149          << "exception '" << exc.what() << "' on proc "
00150          << localProc << std::endl;
00151     MPI_Finalize();
00152     return(-1);
00153   }
00154 
00155   Epetra_LinearProblem balanced_problem(bal_matrix.get(),
00156                                         bal_x.get(), bal_b.get());
00157 
00158 
00159   // Results
00160 
00161   double goalWeight = 1.0 / (double)numProcs;
00162   double bal0, bal1, cutn0, cutn1, cutl0, cutl1;
00163   Isorropia::Epetra::CostDescriber default_costs;
00164 
00165   // Balance and cut quality before partitioning
00166 
00167   ispatest::compute_hypergraph_metrics(*(linprob->GetMatrix()), default_costs, goalWeight,
00168                      bal0, cutn0, cutl0);
00169 
00170   // Balance and cut quality after partitioning
00171 
00172   ispatest::compute_hypergraph_metrics(*bal_matrix, default_costs, goalWeight,
00173                      bal1, cutn1, cutl1);
00174 
00175   if (localProc == 0){
00176     std::cout << "Before partitioning: ";
00177     std::cout << "Balance " << bal0 << " cutN " << cutn0 << " cutL " << cutl0;
00178     std::cout << std::endl;
00179 
00180     std::cout << "After partitioning:  ";
00181     std::cout << "Balance " << bal1 << " cutN " << cutn1 << " cutL " << cutl1;
00182     std::cout << std::endl;
00183   }
00184 
00185   //Finally, delete the pointer objects that we asked to be created.
00186   delete linprob->GetMatrix();
00187   delete linprob->GetLHS();
00188   delete linprob->GetRHS();
00189   delete linprob;
00190 
00191   if (localProc == 0) {
00192     std::cout << std::endl;
00193   }
00194 
00195   MPI_Finalize();
00196 
00197 #else
00198   std::cout << "part_redist: must have both MPI and EPETRA. Make sure Trilinos "
00199     << "is configured with --enable-mpi and --enable-epetra." << std::endl;
00200 #endif
00201 
00202   return(0);
00203 }
00204 
00205 //Below is the implementation of the helper-function that creates the
00206 //poorly-balanced epetra objects for use in the above example program.
00207 
00208 #if defined(HAVE_MPI) && defined(HAVE_EPETRA)
00209 
00210 Epetra_LinearProblem* create_epetra_problem(int numProcs,
00211                                             int localProc,
00212                                             int local_n)
00213 {
00214   if (localProc == 0) {
00215     std::cout << " creating Epetra_CrsMatrix with un-even distribution..."
00216             << std::endl;
00217   }
00218 
00219   //create an Epetra_CrsMatrix with rows spread un-evenly over
00220   //processors.
00221   Epetra_MpiComm comm(MPI_COMM_WORLD);
00222   int global_num_rows = numProcs*local_n;
00223 
00224   int mid_proc = numProcs/2;
00225   bool num_procs_even = numProcs%2==0 ? true : false;
00226 
00227   int adjustment = local_n/2;
00228 
00229   //adjust local_n so that it's not equal on all procs.
00230   if (localProc < mid_proc) {
00231     local_n -= adjustment;
00232   }
00233   else {
00234     local_n += adjustment;
00235   }
00236 
00237   //if numProcs is not an even number, undo the local_n adjustment
00238   //on one proc so that the total will still be correct.
00239   if (localProc == numProcs-1) {
00240     if (num_procs_even == false) {
00241       local_n -= adjustment;
00242     }
00243   }
00244 
00245   //now we're ready to create a row-map.
00246   Epetra_Map rowmap(global_num_rows, local_n, 0, comm);
00247 
00248   //create a matrix
00249   int nnz_per_row = 9;
00250   Epetra_CrsMatrix* matrix =
00251     new Epetra_CrsMatrix(Copy, rowmap, nnz_per_row);
00252 
00253   // Add  rows one-at-a-time
00254   double negOne = -1.0;
00255   double posTwo = 4.0;
00256   for (int i=0; i<local_n; i++) {
00257     int GlobalRow = matrix->GRID(i);
00258     int RowLess1 = GlobalRow - 1;
00259     int RowPlus1 = GlobalRow + 1;
00260     int RowLess2 = GlobalRow - 2;
00261     int RowPlus2 = GlobalRow + 2;
00262     int RowLess3 = GlobalRow - 3;
00263     int RowPlus3 = GlobalRow + 3;
00264     int RowLess4 = GlobalRow - 4;
00265     int RowPlus4 = GlobalRow + 4;
00266 
00267     if (RowLess4>=0) {
00268       matrix->InsertGlobalValues(GlobalRow, 1, &negOne, &RowLess4);
00269     }
00270     if (RowLess3>=0) {
00271       matrix->InsertGlobalValues(GlobalRow, 1, &negOne, &RowLess3);
00272     }
00273     if (RowLess2>=0) {
00274       matrix->InsertGlobalValues(GlobalRow, 1, &negOne, &RowLess2);
00275     }
00276     if (RowLess1>=0) {
00277       matrix->InsertGlobalValues(GlobalRow, 1, &negOne, &RowLess1);
00278     }
00279     if (RowPlus1<global_num_rows) {
00280       matrix->InsertGlobalValues(GlobalRow, 1, &negOne, &RowPlus1);
00281     }
00282     if (RowPlus2<global_num_rows) {
00283       matrix->InsertGlobalValues(GlobalRow, 1, &negOne, &RowPlus2);
00284     }
00285     if (RowPlus3<global_num_rows) {
00286       matrix->InsertGlobalValues(GlobalRow, 1, &negOne, &RowPlus3);
00287     }
00288     if (RowPlus4<global_num_rows) {
00289       matrix->InsertGlobalValues(GlobalRow, 1, &negOne, &RowPlus4);
00290     }
00291 
00292     matrix->InsertGlobalValues(GlobalRow, 1, &posTwo, &GlobalRow);
00293   }
00294 
00295   int err = matrix->FillComplete();
00296   if (err != 0) {
00297     throw Isorropia::Exception("create_epetra_matrix: error in matrix.FillComplete()");
00298   }
00299 
00300   Epetra_Vector* x = new Epetra_Vector(rowmap);
00301   Epetra_Vector* b = new Epetra_Vector(rowmap);
00302   return(new Epetra_LinearProblem(matrix, x, b));
00303 }
00304 
00305 #endif //HAVE_MPI && HAVE_EPETRA
00306