EpetraExt Package Browser (Single Doxygen Collection) Development
EpetraExt_Permutation_impl.h
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 #ifndef EpetraExt_PERMUTATION_IMPL_H
00042 #define EpetraExt_PERMUTATION_IMPL_H
00043 
00044 #include <EpetraExt_ConfigDefs.h>
00045 
00046 #include <EpetraExt_Permutation.h>
00047 
00048 #include <Epetra_Export.h>
00049 #include <Epetra_Map.h>
00050 #include <Epetra_Comm.h>
00051 #include <Epetra_MultiVector.h>
00052 #include <Epetra_CrsGraph.h>
00053 #include <Epetra_CrsMatrix.h>
00054 
00055 namespace EpetraExt {
00056 
00075 template<class T>
00076 struct Perm_traits {
00078   static const char* typeName()
00079   { static const char name[] = "unknown"; return( name ); }
00080 
00088   static T* clone(T* example,
00089       Epetra_DataAccess CV,
00090       const Epetra_BlockMap& map,
00091       int int_argument)
00092   {  return( NULL ); }
00093 
00097   static void replaceMap(T* obj, const Epetra_BlockMap& map)
00098   { std::cerr << "not implemented for unknown type"<<std::endl; }
00099 
00101   template<typename int_type>
00102   static T*
00103   produceColumnPermutation(TPermutation<T, int_type>* perm,
00104          T* srcObj)
00105   { std::cerr << "not implemented for unknown type"<<std::endl; }
00106 
00107 };//struct Perm_traits
00108 
00109 
00110 
00114 template<>
00115 struct Perm_traits<Epetra_CrsMatrix> {
00116 
00118   static const char* typeName()
00119   { static const char name[] = "Epetra_CrsMatrix"; return( name ); }
00120 
00121 
00123   static Epetra_CrsMatrix* clone(Epetra_CrsMatrix* example,
00124          Epetra_DataAccess CV,
00125          const Epetra_BlockMap& map,
00126          int rowLength)
00127   {
00128     //don't need the example object currently...
00129     (void)example;
00130 
00131     //we need a Epetra_Map, rather than a Epetra_BlockMap, to create a
00132     //Epetra_CrsMatrix.
00133 
00134     const Epetra_Map* pointmap =
00135       dynamic_cast<const Epetra_Map*>(&map);
00136     if (pointmap == NULL) {
00137       std::cerr << "dynamic_cast<const Epetra_Map*> failed."<<std::endl;
00138       return(NULL);
00139     }
00140 
00141     return( new Epetra_CrsMatrix(CV, *pointmap, rowLength) );
00142   }
00143 
00144 
00146   static void replaceMap(Epetra_CrsMatrix* mat, const Epetra_BlockMap& map)
00147   { mat->ReplaceRowMap(map); }
00148 
00150   template<typename int_type>
00151   static Epetra_CrsMatrix*
00152   TproduceColumnPermutation(TPermutation<Epetra_CrsMatrix, int_type>* perm,
00153          Epetra_CrsMatrix* srcObj)
00154   {
00155     //First we need to export this permutation to match the column-map of the
00156     //object being column-permuted. (We need to have locally available all
00157     //elements of the permutation corresponding to the local columns of the
00158     //object being permuted.)
00159 
00160     const Epetra_Map& origColMap = srcObj->ColMap();
00161 
00162     TPermutation<Epetra_CrsMatrix, int_type>* colperm =
00163       new TPermutation<Epetra_CrsMatrix, int_type>(origColMap);
00164     colperm->PutValue(0);
00165 
00166     Epetra_Export p_exporter(perm->Map(), origColMap);
00167     colperm->Export(*perm, p_exporter, Add);
00168 
00169     const Epetra_Map& origRowMap = srcObj->RowMap();
00170     int numMyRows = origRowMap.NumMyElements();
00171     int_type* myGlobalRows = 0;
00172     origRowMap.MyGlobalElementsPtr(myGlobalRows);
00173 
00174     //Create the new object, giving it the same map as the original object.
00175 
00176     Epetra_CrsMatrix* result = new Epetra_CrsMatrix(Copy, origRowMap, 1);
00177 
00178     for(int i=0; i<numMyRows; ++i) {
00179       int_type globalRow = myGlobalRows[i];
00180       int len = srcObj->NumGlobalEntries(globalRow);
00181 
00182       int numIndices;
00183       double* src_values = new double[len];
00184       int_type* src_indices = new int_type[len];
00185       int err = srcObj->ExtractGlobalRowCopy(globalRow, len, numIndices,
00186                src_values, src_indices);
00187       if (err < 0 || numIndices != len) {
00188   std::cerr<<"Perm_traits<CrsMatrix>::produceColumnPermutation err("<<err<<") row "
00189       <<globalRow<<", len "<<len<<", numIndices "<<numIndices<<std::endl;
00190       }
00191 
00192       int_type* pindices = new int_type[len];
00193 
00194       const Epetra_BlockMap& pmap = colperm->Map();
00195       int_type* p = colperm->Values();
00196 
00197       for(int j=0; j<len; ++j) {
00198   int_type old_col = src_indices[j];
00199 
00200   int lid = pmap.LID(old_col);
00201   if (lid<0) {
00202     std::cerr << "Perm_traits<CrsMatrix>::permuteColumnIndices GID("<<old_col
00203          <<") not found"<<std::endl;
00204     break;
00205   }
00206 
00207   pindices[j] = p[lid];
00208       }
00209 
00210       err = result->InsertGlobalValues(globalRow, len, src_values, pindices);
00211       if (err < 0) {
00212   std::cerr << "Perm_traits<CrsMatrix>::permuteColumnIndices err("<<err
00213        <<") row "<<globalRow<<std::endl;
00214       }
00215 
00216       delete [] pindices;
00217       delete [] src_indices;
00218       delete [] src_values;
00219     }
00220 
00221     result->FillComplete();
00222 
00223     delete colperm;
00224 
00225     return(result);
00226   }
00227 
00228 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
00229 
00230   static Epetra_CrsMatrix*
00231   produceColumnPermutation(TPermutation<Epetra_CrsMatrix, int>* perm,
00232          Epetra_CrsMatrix* srcObj)
00233   {
00234     return TproduceColumnPermutation<int>(perm, srcObj);
00235   }
00236 #endif
00237 
00238 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
00239 
00240   static Epetra_CrsMatrix*
00241   produceColumnPermutation(TPermutation<Epetra_CrsMatrix, long long>* perm,
00242          Epetra_CrsMatrix* srcObj)
00243   {
00244     return TproduceColumnPermutation<long long>(perm, srcObj);
00245   }
00246 #endif
00247 };//struct Perm_traits<Epetra_CrsMatrix>
00248 
00249 
00250 
00254 template<>
00255 struct Perm_traits<Epetra_CrsGraph> {
00256 
00258   static const char* typeName()
00259   { static const char name[] = "Epetra_CrsGraph"; return( name ); }
00260 
00261 
00263   static Epetra_CrsGraph* clone(Epetra_CrsGraph* example,
00264         Epetra_DataAccess CV,
00265         const Epetra_BlockMap& map,
00266         int rowLength)
00267   {
00268     //don't need the example object currently...
00269     (void)example;
00270 
00271     return( new Epetra_CrsGraph(CV, map, rowLength) );
00272   }
00273 
00274 
00276   static void replaceMap(Epetra_CrsGraph* graph, const Epetra_BlockMap& map)
00277   { graph->ReplaceRowMap(map); }
00278 
00280   template<typename int_type>
00281   static Epetra_CrsGraph*
00282   TproduceColumnPermutation(TPermutation<Epetra_CrsGraph, int_type>* perm,
00283          Epetra_CrsGraph* srcObj)
00284   {
00285     //First we need to export this permutation to match the column-map of the
00286     //object being column-permuted. (We need to have locally available all
00287     //elements of the permutation corresponding to the local columns of the
00288     //object being permuted.)
00289 
00290     const Epetra_BlockMap& origColMap = srcObj->ColMap();
00291 
00292     TPermutation<Epetra_CrsGraph, int_type>* colperm =
00293       new TPermutation<Epetra_CrsGraph, int_type>(origColMap);
00294     colperm->PutValue(0);
00295 
00296     Epetra_Export p_exporter(perm->Map(), origColMap);
00297     colperm->Export(*perm, p_exporter, Add);
00298 
00299     const Epetra_BlockMap& origRowMap = srcObj->RowMap();
00300     int numMyRows = origRowMap.NumMyElements();
00301     int_type* myGlobalRows = 0;
00302     origRowMap.MyGlobalElementsPtr(myGlobalRows);
00303 
00304     //Create the new object, giving it the same map as the original object.
00305 
00306     Epetra_CrsGraph* result = new Epetra_CrsGraph(Copy, origRowMap, 1);
00307 
00308     for(int i=0; i<numMyRows; ++i) {
00309       int_type globalRow = myGlobalRows[i];
00310       int len = srcObj->NumGlobalIndices(globalRow);
00311 
00312       int numIndices;
00313       int_type* src_indices = new int_type[len];
00314       int err = srcObj->ExtractGlobalRowCopy(globalRow, len, numIndices, src_indices);
00315       if (err < 0 || numIndices != len) {
00316   std::cerr<<"Perm_traits<CrsGraph>::produceColumnPermutation err("<<err<<") row "
00317     <<globalRow<<", len "<<len<<", numIndices "<<numIndices<<std::endl;
00318       }
00319 
00320       int_type* pindices = new int_type[len];
00321 
00322       const Epetra_BlockMap& pmap = colperm->Map();
00323       int_type* p = colperm->Values();
00324 
00325       for(int j=0; j<len; ++j) {
00326   int_type old_col = src_indices[j];
00327 
00328   int lid = pmap.LID(old_col);
00329   if (lid<0) {
00330     std::cerr << "Perm_traits<CrsGraph>::permuteColumnIndices GID("<<old_col
00331          <<") not found"<<std::endl;
00332     break;
00333   }
00334 
00335   pindices[j] = p[lid];
00336       }
00337 
00338       err = result->InsertGlobalIndices(globalRow, len, pindices);
00339       if (err < 0) {
00340   std::cerr << "Perm_traits<CrsGraph>::produceColumnPermutation err("<<err
00341        <<") row "<<globalRow<<std::endl;
00342       }
00343 
00344       delete [] pindices;
00345       delete [] src_indices;
00346     }
00347 
00348     result->FillComplete();
00349 
00350     delete colperm;
00351 
00352     return(result);
00353   }
00354 
00355 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
00356 
00357   static Epetra_CrsGraph*
00358   produceColumnPermutation(TPermutation<Epetra_CrsGraph, int>* perm,
00359          Epetra_CrsGraph* srcObj)
00360   {
00361     return TproduceColumnPermutation<int>(perm, srcObj);
00362   }
00363 #endif
00364 
00365 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
00366 
00367   static Epetra_CrsGraph*
00368   produceColumnPermutation(TPermutation<Epetra_CrsGraph, long long>* perm,
00369          Epetra_CrsGraph* srcObj)
00370   {
00371     return TproduceColumnPermutation<long long>(perm, srcObj);
00372   }
00373 #endif
00374 };//struct Perm_traits<Epetra_CrsGraph>
00375 
00376 
00380 template<>
00381 struct Perm_traits<Epetra_MultiVector> {
00382 
00384   static const char* typeName()
00385   { static const char name[] = "Epetra_MultiVector"; return( name ); }
00386 
00387 
00389   static Epetra_MultiVector* clone(Epetra_MultiVector* example,
00390            Epetra_DataAccess CV,
00391            const Epetra_BlockMap& map,
00392            int numVectors)
00393   {
00394     return( new Epetra_MultiVector(map, example->NumVectors()) );
00395   }
00396 
00397 
00399   static void replaceMap(Epetra_MultiVector* mvec, const Epetra_BlockMap& map)
00400   { mvec->ReplaceMap(map); }
00401 
00402 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
00403 
00404   static Epetra_MultiVector*
00405   produceColumnPermutation(Permutation<Epetra_MultiVector>* perm,
00406          Epetra_MultiVector* srcObj)
00407   {
00408     std::cerr << "col-permutation not implemented for Epetra_MultiVector"<<std::endl;
00409     return(NULL);
00410   }
00411 #endif
00412 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
00413 
00414   static Epetra_MultiVector*
00415   produceColumnPermutation(Permutation64<Epetra_MultiVector>* perm,
00416          Epetra_MultiVector* srcObj)
00417   {
00418     std::cerr << "col-permutation not implemented for Epetra_MultiVector"<<std::endl;
00419     return(NULL);
00420   }
00421 #endif
00422 };//struct Perm_traits<Epetra_CrsGraph>
00423 
00424 
00425 //-------------------------------------------------------------------------
00426 //Now the method definitions for the EpetraExt::Permutation class.
00427 //-------------------------------------------------------------------------
00428 
00429 template<typename T, typename int_type>
00430 TPermutation<T, int_type>::TPermutation(Epetra_DataAccess CV,
00431                          const Epetra_BlockMap& map,
00432                          int_type* permutation)
00433   : Epetra_GIDTypeVector<int_type>::impl(CV, map, permutation),
00434     newObj_(NULL),
00435     origObj_(NULL)
00436 {
00437   if (!isTypeSupported()) {
00438     std::cerr << "unsupported type for permutation, aborting" << std::endl;
00439     abort();
00440   }
00441 }
00442 
00443 template<typename T, typename int_type>
00444 TPermutation<T, int_type>::TPermutation(const Epetra_BlockMap& map)
00445   : Epetra_GIDTypeVector<int_type>::impl(map),
00446     newObj_(NULL),
00447     origObj_(NULL)
00448 {
00449   if (!isTypeSupported()) {
00450     std::cerr << "unsupported type for permutation, aborting" << std::endl;
00451     abort();
00452   }
00453 }
00454 
00455 template<typename T, typename int_type>
00456 TPermutation<T, int_type>::TPermutation(const TPermutation& src)
00457   : Epetra_GIDTypeVector<int_type>::impl((const typename Epetra_GIDTypeVector<int_type>::impl&)src),
00458     newObj_(NULL),
00459     origObj_(NULL)
00460 {
00461   if (!isTypeSupported()) {
00462     std::cerr << "unsupported type for permutation, aborting" << std::endl;
00463     abort();
00464   }
00465 }
00466 
00467 template<typename T, typename int_type>
00468 TPermutation<T, int_type>::~TPermutation()
00469 {
00470   if (newObj_ != NULL) delete newObj_;
00471 }
00472 
00473 template<typename T, typename int_type>
00474 bool TPermutation<T, int_type>::isTypeSupported()
00475 {
00476   const char* type_name = Perm_traits<T>::typeName();
00477   if (!strcmp(type_name, "unknown")) {
00478     return(false);
00479   }
00480 
00481   return( true );
00482 }
00483 
00484 template<typename T, typename int_type>
00485 typename TPermutation<T, int_type>::OutputRef
00486 TPermutation<T, int_type>::operator()( typename TPermutation<T, int_type>::InputRef orig )
00487 {
00488   //In this function we're going to produce a new object which is a
00489   //row-permutation of the input object (orig).
00490   //
00491   //Our permutation inherits IntVector, and the permutation is defined by the
00492   //contents of the integer vector 'p', such that if p[i] = j then row i of
00493   //the input object becomes row j of the permuted object.
00494   //
00495   //The permutation is accomplished by creating a map defined by the
00496   //permutation, then using an Epetra_Export operation to move data from the
00497   //input object into the permuted object.
00498   //
00499   //The permutation may be global. In other words, the rows of the object may
00500   //be arbitrarily rearranged, including across processors.
00501   //
00502 
00503   origObj_ = &orig;
00504 
00505   //The 'Map()' accessor returns Epetra_DistObject::Map() for CrsGraph and
00506   //CrsMatrix, which turns out to be the RowMap() for those objects. For
00507   //MultiVector it returns the correct object because MultiVectors only have
00508   //one map.
00509 
00510   const Epetra_BlockMap& origMap = orig.Map();
00511 
00512   //Create an Epetra_Map representing the permutation.
00513 
00514   Epetra_Map* pmap = new Epetra_Map((int_type) Epetra_DistObject::Map().NumGlobalPoints64(),
00515             Epetra_DistObject::Map().NumMyPoints(),
00516             Epetra_GIDTypeVector<int_type>::impl::Values(),
00517             (int_type) Epetra_DistObject::Map().IndexBase64(),
00518             Epetra_DistObject::Map().Comm());
00519 
00520   TPermutation* p = this;
00521 
00522   //Next check that the maps are compatible. If they aren't, we'll redistribute
00523   //the permutation to match the distribution of the input object.
00524 
00525   if (!pmap->PointSameAs(origMap)) {
00526     Epetra_Export p_exporter(Epetra_DistObject::Map(), origMap);
00527     TPermutation* newp = new TPermutation(origMap);
00528     newp->Export(*p, p_exporter, Add);
00529     p = newp;
00530 
00531     delete pmap;
00532     pmap = new Epetra_Map((int_type) p->Map().NumGlobalPoints64(),
00533         p->Map().NumMyPoints(),
00534         p->Values(),
00535         (int_type) p->Map().IndexBase64(),
00536         p->Map().Comm());
00537   }
00538 
00539   //Create the new object, initially giving it the map defined by the
00540   //permutation.
00541 
00542   newObj_ = Perm_traits<T>::clone(origObj_, Copy, *pmap, 1);
00543 
00544   //Create an exporter which will export data from the original object to the
00545   //permuted object.
00546 
00547   Epetra_Export exporter(origMap, *pmap);
00548 
00549   //Now export the original object to the permuted object.
00550 
00551   newObj_->Export(orig, exporter, Add);
00552 
00553   //Now, since the export operation moved not only row-contents but also
00554   //row-numbering, we need to replace the permuted row-numbering with the
00555   //original row-numbering. We do this by replacing the permuted map with
00556   //the original row-map.
00557 
00558   Perm_traits<T>::replaceMap(newObj_, origMap);
00559 
00560   delete pmap;
00561 
00562   if (p != this) {
00563     delete p; //delete "newp" created if the PointSameAs test failed above
00564   }
00565 
00566   return( *newObj_ );
00567 }
00568 
00569 template<typename T, typename int_type>
00570 typename TPermutation<T, int_type>::OutputRef
00571 TPermutation<T, int_type>::operator()( typename TPermutation<T, int_type>::InputRef orig,
00572           bool column_permutation )
00573 {
00574   origObj_ = &orig;
00575   newObj_ = NULL;
00576 
00577   if (!column_permutation) {
00578     return( operator()(orig) );
00579   }
00580 
00581   if (strcmp("Epetra_CrsMatrix", Perm_traits<T>::typeName()) &&
00582       strcmp("Epetra_CrsGraph", Perm_traits<T>::typeName())) {
00583     std::cerr << "Permutation: column-permutation only implemented for"
00584    << "CrsMatrix and CrsGraph." << std::endl;
00585     assert(0);
00586   }
00587 
00588   newObj_ = Perm_traits<T>::produceColumnPermutation(this, &orig);
00589 
00590   return( *newObj_ );
00591 }
00592 
00593 } // namespace EpetraExt
00594 
00595 #endif //EpetraExt_PERMUTATION_IMPL_H
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines