EpetraExt Package Browser (Single Doxygen Collection)  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
EpetraExt_Permutation_impl.h
Go to the documentation of this file.
1 //@HEADER
2 // ***********************************************************************
3 //
4 // EpetraExt: Epetra Extended - Linear Algebra Services Package
5 // Copyright (2011) Sandia Corporation
6 //
7 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8 // the U.S. Government retains certain rights in this software.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38 //
39 // ***********************************************************************
40 //@HEADER
41 #ifndef EpetraExt_PERMUTATION_IMPL_H
42 #define EpetraExt_PERMUTATION_IMPL_H
43 
44 #include <EpetraExt_ConfigDefs.h>
45 
46 #include <EpetraExt_Permutation.h>
47 
48 #include <Epetra_Export.h>
49 #include <Epetra_Map.h>
50 #include <Epetra_Comm.h>
51 #include <Epetra_MultiVector.h>
52 #include <Epetra_CrsGraph.h>
53 #include <Epetra_CrsMatrix.h>
54 #include <Epetra_GIDTypeVector.h>
55 
56 namespace EpetraExt {
57 
76 template<class T>
77 struct Perm_traits {
79  static const char* typeName()
80  { static const char name[] = "unknown"; return( name ); }
81 
89  static T* clone(T* example,
91  const Epetra_BlockMap& map,
92  int int_argument)
93  { return( NULL ); }
94 
98  static void replaceMap(T* obj, const Epetra_BlockMap& map)
99  { std::cerr << "not implemented for unknown type"<<std::endl; }
100 
102  template<typename int_type>
103  static T*
105  T* srcObj)
106  { std::cerr << "not implemented for unknown type"<<std::endl; }
107 
108 };//struct Perm_traits
109 
110 
111 
115 template<>
117 
119  static const char* typeName()
120  { static const char name[] = "Epetra_CrsMatrix"; return( name ); }
121 
122 
126  const Epetra_BlockMap& map,
127  int rowLength)
128  {
129  //don't need the example object currently...
130  (void)example;
131 
132  //we need a Epetra_Map, rather than a Epetra_BlockMap, to create a
133  //Epetra_CrsMatrix.
134 
135  const Epetra_Map* pointmap =
136  dynamic_cast<const Epetra_Map*>(&map);
137  if (pointmap == NULL) {
138  std::cerr << "dynamic_cast<const Epetra_Map*> failed."<<std::endl;
139  return(NULL);
140  }
141 
142  return( new Epetra_CrsMatrix(CV, *pointmap, rowLength) );
143  }
144 
145 
147  static void replaceMap(Epetra_CrsMatrix* mat, const Epetra_BlockMap& map)
148  { mat->ReplaceRowMap(map); }
149 
151  template<typename int_type>
152  static Epetra_CrsMatrix*
154  Epetra_CrsMatrix* srcObj)
155  {
156  //First we need to export this permutation to match the column-map of the
157  //object being column-permuted. (We need to have locally available all
158  //elements of the permutation corresponding to the local columns of the
159  //object being permuted.)
160 
161  const Epetra_Map& origColMap = srcObj->ColMap();
162 
165  colperm->PutValue(0);
166 
167  Epetra_Export p_exporter(perm->Map(), origColMap);
168  colperm->Export(*perm, p_exporter, Add);
169 
170  const Epetra_Map& origRowMap = srcObj->RowMap();
171  int numMyRows = origRowMap.NumMyElements();
172  int_type* myGlobalRows = 0;
173  origRowMap.MyGlobalElementsPtr(myGlobalRows);
174 
175  //Create the new object, giving it the same map as the original object.
176 
177  Epetra_CrsMatrix* result = new Epetra_CrsMatrix(Copy, origRowMap, 1);
178 
179  for(int i=0; i<numMyRows; ++i) {
180  int_type globalRow = myGlobalRows[i];
181  int len = srcObj->NumGlobalEntries(globalRow);
182 
183  int numIndices;
184  double* src_values = new double[len];
185  int_type* src_indices = new int_type[len];
186  int err = srcObj->ExtractGlobalRowCopy(globalRow, len, numIndices,
187  src_values, src_indices);
188  if (err < 0 || numIndices != len) {
189  std::cerr<<"Perm_traits<CrsMatrix>::produceColumnPermutation err("<<err<<") row "
190  <<globalRow<<", len "<<len<<", numIndices "<<numIndices<<std::endl;
191  }
192 
193  int_type* pindices = new int_type[len];
194 
195  const Epetra_BlockMap& pmap = colperm->Map();
196  int_type* p = colperm->Values();
197 
198  for(int j=0; j<len; ++j) {
199  int_type old_col = src_indices[j];
200 
201  int lid = pmap.LID(old_col);
202  if (lid<0) {
203  std::cerr << "Perm_traits<CrsMatrix>::permuteColumnIndices GID("<<old_col
204  <<") not found"<<std::endl;
205  break;
206  }
207 
208  pindices[j] = p[lid];
209  }
210 
211  err = result->InsertGlobalValues(globalRow, len, src_values, pindices);
212  if (err < 0) {
213  std::cerr << "Perm_traits<CrsMatrix>::permuteColumnIndices err("<<err
214  <<") row "<<globalRow<<std::endl;
215  }
216 
217  delete [] pindices;
218  delete [] src_indices;
219  delete [] src_values;
220  }
221 
222  result->FillComplete();
223 
224  delete colperm;
225 
226  return(result);
227  }
228 
229 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
230 
231  static Epetra_CrsMatrix*
233  Epetra_CrsMatrix* srcObj)
234  {
235  return TproduceColumnPermutation<int>(perm, srcObj);
236  }
237 #endif
238 
239 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
240 
241  static Epetra_CrsMatrix*
243  Epetra_CrsMatrix* srcObj)
244  {
245  return TproduceColumnPermutation<long long>(perm, srcObj);
246  }
247 #endif
248 };//struct Perm_traits<Epetra_CrsMatrix>
249 
250 
251 
255 template<>
257 
259  static const char* typeName()
260  { static const char name[] = "Epetra_CrsGraph"; return( name ); }
261 
262 
266  const Epetra_BlockMap& map,
267  int rowLength)
268  {
269  //don't need the example object currently...
270  (void)example;
271 
272  return( new Epetra_CrsGraph(CV, map, rowLength) );
273  }
274 
275 
277  static void replaceMap(Epetra_CrsGraph* graph, const Epetra_BlockMap& map)
278  { graph->ReplaceRowMap(map); }
279 
281  template<typename int_type>
282  static Epetra_CrsGraph*
284  Epetra_CrsGraph* srcObj)
285  {
286  //First we need to export this permutation to match the column-map of the
287  //object being column-permuted. (We need to have locally available all
288  //elements of the permutation corresponding to the local columns of the
289  //object being permuted.)
290 
291  const Epetra_BlockMap& origColMap = srcObj->ColMap();
292 
295  colperm->PutValue(0);
296 
297  Epetra_Export p_exporter(perm->Map(), origColMap);
298  colperm->Export(*perm, p_exporter, Add);
299 
300  const Epetra_BlockMap& origRowMap = srcObj->RowMap();
301  int numMyRows = origRowMap.NumMyElements();
302  int_type* myGlobalRows = 0;
303  origRowMap.MyGlobalElementsPtr(myGlobalRows);
304 
305  //Create the new object, giving it the same map as the original object.
306 
307  Epetra_CrsGraph* result = new Epetra_CrsGraph(Copy, origRowMap, 1);
308 
309  for(int i=0; i<numMyRows; ++i) {
310  int_type globalRow = myGlobalRows[i];
311  int len = srcObj->NumGlobalIndices(globalRow);
312 
313  int numIndices;
314  int_type* src_indices = new int_type[len];
315  int err = srcObj->ExtractGlobalRowCopy(globalRow, len, numIndices, src_indices);
316  if (err < 0 || numIndices != len) {
317  std::cerr<<"Perm_traits<CrsGraph>::produceColumnPermutation err("<<err<<") row "
318  <<globalRow<<", len "<<len<<", numIndices "<<numIndices<<std::endl;
319  }
320 
321  int_type* pindices = new int_type[len];
322 
323  const Epetra_BlockMap& pmap = colperm->Map();
324  int_type* p = colperm->Values();
325 
326  for(int j=0; j<len; ++j) {
327  int_type old_col = src_indices[j];
328 
329  int lid = pmap.LID(old_col);
330  if (lid<0) {
331  std::cerr << "Perm_traits<CrsGraph>::permuteColumnIndices GID("<<old_col
332  <<") not found"<<std::endl;
333  break;
334  }
335 
336  pindices[j] = p[lid];
337  }
338 
339  err = result->InsertGlobalIndices(globalRow, len, pindices);
340  if (err < 0) {
341  std::cerr << "Perm_traits<CrsGraph>::produceColumnPermutation err("<<err
342  <<") row "<<globalRow<<std::endl;
343  }
344 
345  delete [] pindices;
346  delete [] src_indices;
347  }
348 
349  result->FillComplete();
350 
351  delete colperm;
352 
353  return(result);
354  }
355 
356 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
357 
358  static Epetra_CrsGraph*
360  Epetra_CrsGraph* srcObj)
361  {
362  return TproduceColumnPermutation<int>(perm, srcObj);
363  }
364 #endif
365 
366 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
367 
368  static Epetra_CrsGraph*
370  Epetra_CrsGraph* srcObj)
371  {
372  return TproduceColumnPermutation<long long>(perm, srcObj);
373  }
374 #endif
375 };//struct Perm_traits<Epetra_CrsGraph>
376 
377 
381 template<>
383 
385  static const char* typeName()
386  { static const char name[] = "Epetra_MultiVector"; return( name ); }
387 
388 
391  Epetra_DataAccess /* CV */,
392  const Epetra_BlockMap& map,
393  int /* numVectors */)
394  {
395  return( new Epetra_MultiVector(map, example->NumVectors()) );
396  }
397 
398 
400  static void replaceMap(Epetra_MultiVector* mvec, const Epetra_BlockMap& map)
401  { mvec->ReplaceMap(map); }
402 
403 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
404 
405  static Epetra_MultiVector*
407  Epetra_MultiVector* /* srcObj */)
408  {
409  std::cerr << "col-permutation not implemented for Epetra_MultiVector"<<std::endl;
410  return(NULL);
411  }
412 #endif
413 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
414 
415  static Epetra_MultiVector*
417  Epetra_MultiVector* /* srcObj */)
418  {
419  std::cerr << "col-permutation not implemented for Epetra_MultiVector"<<std::endl;
420  return(NULL);
421  }
422 #endif
423 };//struct Perm_traits<Epetra_CrsGraph>
424 
425 
426 //-------------------------------------------------------------------------
427 //Now the method definitions for the EpetraExt::Permutation class.
428 //-------------------------------------------------------------------------
429 
430 template<typename T, typename int_type>
432  const Epetra_BlockMap& map,
433  int_type* permutation)
434  : Epetra_GIDTypeVector<int_type>::impl(CV, map, permutation),
435  newObj_(NULL),
436  origObj_(NULL)
437 {
438  if (!isTypeSupported()) {
439  std::cerr << "unsupported type for permutation, aborting" << std::endl;
440  abort();
441  }
442 }
443 
444 template<typename T, typename int_type>
446  : Epetra_GIDTypeVector<int_type>::impl(map),
447  newObj_(NULL),
448  origObj_(NULL)
449 {
450  if (!isTypeSupported()) {
451  std::cerr << "unsupported type for permutation, aborting" << std::endl;
452  abort();
453  }
454 }
455 
456 template<typename T, typename int_type>
458  : Epetra_GIDTypeVector<int_type>::impl((const typename Epetra_GIDTypeVector<int_type>::impl&)src),
459  newObj_(NULL),
460  origObj_(NULL)
461 {
462  if (!isTypeSupported()) {
463  std::cerr << "unsupported type for permutation, aborting" << std::endl;
464  abort();
465  }
466 }
467 
468 template<typename T, typename int_type>
470 {
471  if (newObj_ != NULL) delete newObj_;
472 }
473 
474 template<typename T, typename int_type>
476 {
477  const char* type_name = Perm_traits<T>::typeName();
478  if (!strcmp(type_name, "unknown")) {
479  return(false);
480  }
481 
482  return( true );
483 }
484 
485 template<typename T, typename int_type>
488 {
489  //In this function we're going to produce a new object which is a
490  //row-permutation of the input object (orig).
491  //
492  //Our permutation inherits IntVector, and the permutation is defined by the
493  //contents of the integer vector 'p', such that if p[i] = j then row i of
494  //the input object becomes row j of the permuted object.
495  //
496  //The permutation is accomplished by creating a map defined by the
497  //permutation, then using an Epetra_Export operation to move data from the
498  //input object into the permuted object.
499  //
500  //The permutation may be global. In other words, the rows of the object may
501  //be arbitrarily rearranged, including across processors.
502  //
503 
504  origObj_ = &orig;
505 
506  //The 'Map()' accessor returns Epetra_DistObject::Map() for CrsGraph and
507  //CrsMatrix, which turns out to be the RowMap() for those objects. For
508  //MultiVector it returns the correct object because MultiVectors only have
509  //one map.
510 
511  const Epetra_BlockMap& origMap = orig.Map();
512 
513  //Create an Epetra_Map representing the permutation.
514 
515  Epetra_Map* pmap = new Epetra_Map((int_type) Epetra_DistObject::Map().NumGlobalPoints64(),
516  Epetra_DistObject::Map().NumMyPoints(),
518  (int_type) Epetra_DistObject::Map().IndexBase64(),
519  Epetra_DistObject::Map().Comm());
520 
521  TPermutation* p = this;
522 
523  //Next check that the maps are compatible. If they aren't, we'll redistribute
524  //the permutation to match the distribution of the input object.
525 
526  if (!pmap->PointSameAs(origMap)) {
527  Epetra_Export p_exporter(Epetra_DistObject::Map(), origMap);
528  TPermutation* newp = new TPermutation(origMap);
529  newp->Export(*p, p_exporter, Add);
530  p = newp;
531 
532  delete pmap;
533  pmap = new Epetra_Map((int_type) p->Map().NumGlobalPoints64(),
534  p->Map().NumMyPoints(),
535  p->Values(),
536  (int_type) p->Map().IndexBase64(),
537  p->Map().Comm());
538  }
539 
540  //Create the new object, initially giving it the map defined by the
541  //permutation.
542 
543  newObj_ = Perm_traits<T>::clone(origObj_, Copy, *pmap, 1);
544 
545  //Create an exporter which will export data from the original object to the
546  //permuted object.
547 
548  Epetra_Export exporter(origMap, *pmap);
549 
550  //Now export the original object to the permuted object.
551 
552  newObj_->Export(orig, exporter, Add);
553 
554  //Now, since the export operation moved not only row-contents but also
555  //row-numbering, we need to replace the permuted row-numbering with the
556  //original row-numbering. We do this by replacing the permuted map with
557  //the original row-map.
558 
560 
561  delete pmap;
562 
563  if (p != this) {
564  delete p; //delete "newp" created if the PointSameAs test failed above
565  }
566 
567  return( *newObj_ );
568 }
569 
570 template<typename T, typename int_type>
571 typename TPermutation<T, int_type>::OutputRef
572 TPermutation<T, int_type>::operator()( typename TPermutation<T, int_type>::InputRef orig,
573  bool column_permutation )
574 {
575  origObj_ = &orig;
576  newObj_ = NULL;
577 
578  if (!column_permutation) {
579  return( operator()(orig) );
580  }
581 
582  if (strcmp("Epetra_CrsMatrix", Perm_traits<T>::typeName()) &&
583  strcmp("Epetra_CrsGraph", Perm_traits<T>::typeName())) {
584  std::cerr << "Permutation: column-permutation only implemented for"
585  << "CrsMatrix and CrsGraph." << std::endl;
586  assert(0);
587  }
588 
590 
591  return( *newObj_ );
592 }
593 
594 } // namespace EpetraExt
595 
596 #endif //EpetraExt_PERMUTATION_IMPL_H
static Epetra_CrsGraph * produceColumnPermutation(TPermutation< Epetra_CrsGraph, long long > *perm, Epetra_CrsGraph *srcObj)
return new object which is a column-permutation of srcObj
bool PointSameAs(const Epetra_BlockMap &Map) const
int NumGlobalEntries(long long Row) const
static Epetra_MultiVector * produceColumnPermutation(Permutation64< Epetra_MultiVector > *, Epetra_MultiVector *)
permute column-indices within a specified row, if applicable
static Epetra_CrsGraph * TproduceColumnPermutation(TPermutation< Epetra_CrsGraph, int_type > *perm, Epetra_CrsGraph *srcObj)
return new object which is a column-permutation of srcObj
static Epetra_CrsGraph * clone(Epetra_CrsGraph *example, Epetra_DataAccess CV, const Epetra_BlockMap &map, int rowLength)
clone implementation
static void replaceMap(Epetra_MultiVector *mvec, const Epetra_BlockMap &map)
replaceMap implementation
Transform< T, T >::NewTypeRef Transform_Composite< T >typename Transform< T, T >::OriginalTypeRef orig this origObj_
int NumGlobalIndices(long long Row) const
int ExtractGlobalRowCopy(int_type Row, int Length, int &NumEntries, double *values, int_type *Indices) const
static T * clone(T *example, Epetra_DataAccess CV, const Epetra_BlockMap &map, int int_argument)
clone function accepts an example of the object being cloned, and enough constructor arguments to be ...
virtual int InsertGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
static const char * typeName()
typeName implementation
int ReplaceRowMap(const Epetra_BlockMap &newmap)
const Epetra_BlockMap & ColMap() const
static const char * typeName()
typeName implementation
const Epetra_Map & ColMap() const
static Epetra_CrsMatrix * produceColumnPermutation(TPermutation< Epetra_CrsMatrix, long long > *perm, Epetra_CrsMatrix *srcObj)
return new object, which is a column-permutation of srcObj
int InsertGlobalIndices(int_type GlobalRow, int NumIndices, int_type *Indices)
int FillComplete(bool OptimizeDataStorage=true)
static Epetra_CrsMatrix * produceColumnPermutation(TPermutation< Epetra_CrsMatrix, int > *perm, Epetra_CrsMatrix *srcObj)
return new object, which is a column-permutation of srcObj
const Epetra_Map & RowMap() const
int NumMyElements() const
virtual ~TPermutation()
Destructor.
Permutation stores and describes a permutation matrix P.
static void replaceMap(Epetra_CrsMatrix *mat, const Epetra_BlockMap &map)
replaceMap implementation
static Epetra_CrsMatrix * clone(Epetra_CrsMatrix *example, Epetra_DataAccess CV, const Epetra_BlockMap &map, int rowLength)
clone implementation
int ReplaceRowMap(const Epetra_BlockMap &newmap)
static Epetra_MultiVector * produceColumnPermutation(Permutation< Epetra_MultiVector > *, Epetra_MultiVector *)
permute column-indices within a specified row, if applicable
virtual const Epetra_BlockMap & Map() const =0
static Epetra_CrsGraph * produceColumnPermutation(TPermutation< Epetra_CrsGraph, int > *perm, Epetra_CrsGraph *srcObj)
return new object which is a column-permutation of srcObj
static Epetra_MultiVector * clone(Epetra_MultiVector *example, Epetra_DataAccess, const Epetra_BlockMap &map, int)
clone implementation
OutputRef operator()(InputRef orig)
This method creates a new object which is a permuted copy of the input argument.
static const char * typeName()
typeName implementation
static T * produceColumnPermutation(TPermutation< T, int_type > *perm, T *srcObj)
return new object, which is a column-permutation of srcObj
const Epetra_BlockMap & RowMap() const
int MyGlobalElementsPtr(int *&MyGlobalElementList) const
int LID(int GID) const
static void replaceMap(T *obj, const Epetra_BlockMap &map)
replace the object&#39;s row-map (or if it&#39;s not a matrix, replace its only map)
TPermutation(Epetra_DataAccess CV, const Epetra_BlockMap &map, int_type *permutation)
Constructor.
EpetraExt::SameTypeTransform< T >::TransformTypeRef InputRef
static const char * typeName()
return a std::string name for the object type
Define some traits to make it easier to deal with template-parameters which are objects to be permute...
int ExtractGlobalRowCopy(int_type Row, int LenOfIndices, int &NumIndices, int_type *Indices) const
Epetra_DataAccess
EpetraExt::SameTypeTransform< T >::TransformTypeRef OutputRef
static void replaceMap(Epetra_CrsGraph *graph, const Epetra_BlockMap &map)
replaceMap implementation
static Epetra_CrsMatrix * TproduceColumnPermutation(TPermutation< Epetra_CrsMatrix, int_type > *perm, Epetra_CrsMatrix *srcObj)
return new object, which is a column-permutation of srcObj