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