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