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