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 #if defined(EpetraExt_SHOW_DEPRECATED_WARNINGS)
45 #ifdef __GNUC__
46 #warning "The EpetraExt package is deprecated"
47 #endif
48 #endif
49 
50 #include <EpetraExt_ConfigDefs.h>
51 
52 #include <EpetraExt_Permutation.h>
53 
54 #include <Epetra_Export.h>
55 #include <Epetra_Map.h>
56 #include <Epetra_Comm.h>
57 #include <Epetra_MultiVector.h>
58 #include <Epetra_CrsGraph.h>
59 #include <Epetra_CrsMatrix.h>
60 #include <Epetra_GIDTypeVector.h>
61 
62 namespace EpetraExt {
63 
82 template<class T>
83 struct Perm_traits {
85  static const char* typeName()
86  { static const char name[] = "unknown"; return( name ); }
87 
95  static T* clone(T* example,
97  const Epetra_BlockMap& map,
98  int int_argument)
99  { return( NULL ); }
100 
104  static void replaceMap(T* obj, const Epetra_BlockMap& map)
105  { std::cerr << "not implemented for unknown type"<<std::endl; }
106 
108  template<typename int_type>
109  static T*
111  T* srcObj)
112  { std::cerr << "not implemented for unknown type"<<std::endl; }
113 
114 };//struct Perm_traits
115 
116 
117 
121 template<>
123 
125  static const char* typeName()
126  { static const char name[] = "Epetra_CrsMatrix"; return( name ); }
127 
128 
132  const Epetra_BlockMap& map,
133  int rowLength)
134  {
135  //don't need the example object currently...
136  (void)example;
137 
138  //we need a Epetra_Map, rather than a Epetra_BlockMap, to create a
139  //Epetra_CrsMatrix.
140 
141  const Epetra_Map* pointmap =
142  dynamic_cast<const Epetra_Map*>(&map);
143  if (pointmap == NULL) {
144  std::cerr << "dynamic_cast<const Epetra_Map*> failed."<<std::endl;
145  return(NULL);
146  }
147 
148  return( new Epetra_CrsMatrix(CV, *pointmap, rowLength) );
149  }
150 
151 
153  static void replaceMap(Epetra_CrsMatrix* mat, const Epetra_BlockMap& map)
154  { mat->ReplaceRowMap(map); }
155 
157  template<typename int_type>
158  static Epetra_CrsMatrix*
160  Epetra_CrsMatrix* srcObj)
161  {
162  //First we need to export this permutation to match the column-map of the
163  //object being column-permuted. (We need to have locally available all
164  //elements of the permutation corresponding to the local columns of the
165  //object being permuted.)
166 
167  const Epetra_Map& origColMap = srcObj->ColMap();
168 
171  colperm->PutValue(0);
172 
173  Epetra_Export p_exporter(perm->Map(), origColMap);
174  colperm->Export(*perm, p_exporter, Add);
175 
176  const Epetra_Map& origRowMap = srcObj->RowMap();
177  int numMyRows = origRowMap.NumMyElements();
178  int_type* myGlobalRows = 0;
179  origRowMap.MyGlobalElementsPtr(myGlobalRows);
180 
181  //Create the new object, giving it the same map as the original object.
182 
183  Epetra_CrsMatrix* result = new Epetra_CrsMatrix(Copy, origRowMap, 1);
184 
185  for(int i=0; i<numMyRows; ++i) {
186  int_type globalRow = myGlobalRows[i];
187  int len = srcObj->NumGlobalEntries(globalRow);
188 
189  int numIndices;
190  double* src_values = new double[len];
191  int_type* src_indices = new int_type[len];
192  int err = srcObj->ExtractGlobalRowCopy(globalRow, len, numIndices,
193  src_values, src_indices);
194  if (err < 0 || numIndices != len) {
195  std::cerr<<"Perm_traits<CrsMatrix>::produceColumnPermutation err("<<err<<") row "
196  <<globalRow<<", len "<<len<<", numIndices "<<numIndices<<std::endl;
197  }
198 
199  int_type* pindices = new int_type[len];
200 
201  const Epetra_BlockMap& pmap = colperm->Map();
202  int_type* p = colperm->Values();
203 
204  for(int j=0; j<len; ++j) {
205  int_type old_col = src_indices[j];
206 
207  int lid = pmap.LID(old_col);
208  if (lid<0) {
209  std::cerr << "Perm_traits<CrsMatrix>::permuteColumnIndices GID("<<old_col
210  <<") not found"<<std::endl;
211  break;
212  }
213 
214  pindices[j] = p[lid];
215  }
216 
217  err = result->InsertGlobalValues(globalRow, len, src_values, pindices);
218  if (err < 0) {
219  std::cerr << "Perm_traits<CrsMatrix>::permuteColumnIndices err("<<err
220  <<") row "<<globalRow<<std::endl;
221  }
222 
223  delete [] pindices;
224  delete [] src_indices;
225  delete [] src_values;
226  }
227 
228  result->FillComplete();
229 
230  delete colperm;
231 
232  return(result);
233  }
234 
235 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
236 
237  static Epetra_CrsMatrix*
239  Epetra_CrsMatrix* srcObj)
240  {
241  return TproduceColumnPermutation<int>(perm, srcObj);
242  }
243 #endif
244 
245 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
246 
247  static Epetra_CrsMatrix*
249  Epetra_CrsMatrix* srcObj)
250  {
251  return TproduceColumnPermutation<long long>(perm, srcObj);
252  }
253 #endif
254 };//struct Perm_traits<Epetra_CrsMatrix>
255 
256 
257 
261 template<>
263 
265  static const char* typeName()
266  { static const char name[] = "Epetra_CrsGraph"; return( name ); }
267 
268 
272  const Epetra_BlockMap& map,
273  int rowLength)
274  {
275  //don't need the example object currently...
276  (void)example;
277 
278  return( new Epetra_CrsGraph(CV, map, rowLength) );
279  }
280 
281 
283  static void replaceMap(Epetra_CrsGraph* graph, const Epetra_BlockMap& map)
284  { graph->ReplaceRowMap(map); }
285 
287  template<typename int_type>
288  static Epetra_CrsGraph*
290  Epetra_CrsGraph* srcObj)
291  {
292  //First we need to export this permutation to match the column-map of the
293  //object being column-permuted. (We need to have locally available all
294  //elements of the permutation corresponding to the local columns of the
295  //object being permuted.)
296 
297  const Epetra_BlockMap& origColMap = srcObj->ColMap();
298 
301  colperm->PutValue(0);
302 
303  Epetra_Export p_exporter(perm->Map(), origColMap);
304  colperm->Export(*perm, p_exporter, Add);
305 
306  const Epetra_BlockMap& origRowMap = srcObj->RowMap();
307  int numMyRows = origRowMap.NumMyElements();
308  int_type* myGlobalRows = 0;
309  origRowMap.MyGlobalElementsPtr(myGlobalRows);
310 
311  //Create the new object, giving it the same map as the original object.
312 
313  Epetra_CrsGraph* result = new Epetra_CrsGraph(Copy, origRowMap, 1);
314 
315  for(int i=0; i<numMyRows; ++i) {
316  int_type globalRow = myGlobalRows[i];
317  int len = srcObj->NumGlobalIndices(globalRow);
318 
319  int numIndices;
320  int_type* src_indices = new int_type[len];
321  int err = srcObj->ExtractGlobalRowCopy(globalRow, len, numIndices, src_indices);
322  if (err < 0 || numIndices != len) {
323  std::cerr<<"Perm_traits<CrsGraph>::produceColumnPermutation err("<<err<<") row "
324  <<globalRow<<", len "<<len<<", numIndices "<<numIndices<<std::endl;
325  }
326 
327  int_type* pindices = new int_type[len];
328 
329  const Epetra_BlockMap& pmap = colperm->Map();
330  int_type* p = colperm->Values();
331 
332  for(int j=0; j<len; ++j) {
333  int_type old_col = src_indices[j];
334 
335  int lid = pmap.LID(old_col);
336  if (lid<0) {
337  std::cerr << "Perm_traits<CrsGraph>::permuteColumnIndices GID("<<old_col
338  <<") not found"<<std::endl;
339  break;
340  }
341 
342  pindices[j] = p[lid];
343  }
344 
345  err = result->InsertGlobalIndices(globalRow, len, pindices);
346  if (err < 0) {
347  std::cerr << "Perm_traits<CrsGraph>::produceColumnPermutation err("<<err
348  <<") row "<<globalRow<<std::endl;
349  }
350 
351  delete [] pindices;
352  delete [] src_indices;
353  }
354 
355  result->FillComplete();
356 
357  delete colperm;
358 
359  return(result);
360  }
361 
362 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
363 
364  static Epetra_CrsGraph*
366  Epetra_CrsGraph* srcObj)
367  {
368  return TproduceColumnPermutation<int>(perm, srcObj);
369  }
370 #endif
371 
372 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
373 
374  static Epetra_CrsGraph*
376  Epetra_CrsGraph* srcObj)
377  {
378  return TproduceColumnPermutation<long long>(perm, srcObj);
379  }
380 #endif
381 };//struct Perm_traits<Epetra_CrsGraph>
382 
383 
387 template<>
389 
391  static const char* typeName()
392  { static const char name[] = "Epetra_MultiVector"; return( name ); }
393 
394 
397  Epetra_DataAccess /* CV */,
398  const Epetra_BlockMap& map,
399  int /* numVectors */)
400  {
401  return( new Epetra_MultiVector(map, example->NumVectors()) );
402  }
403 
404 
406  static void replaceMap(Epetra_MultiVector* mvec, const Epetra_BlockMap& map)
407  { mvec->ReplaceMap(map); }
408 
409 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
410 
411  static Epetra_MultiVector*
413  Epetra_MultiVector* /* srcObj */)
414  {
415  std::cerr << "col-permutation not implemented for Epetra_MultiVector"<<std::endl;
416  return(NULL);
417  }
418 #endif
419 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
420 
421  static Epetra_MultiVector*
423  Epetra_MultiVector* /* srcObj */)
424  {
425  std::cerr << "col-permutation not implemented for Epetra_MultiVector"<<std::endl;
426  return(NULL);
427  }
428 #endif
429 };//struct Perm_traits<Epetra_CrsGraph>
430 
431 
432 //-------------------------------------------------------------------------
433 //Now the method definitions for the EpetraExt::Permutation class.
434 //-------------------------------------------------------------------------
435 
436 template<typename T, typename int_type>
438  const Epetra_BlockMap& map,
439  int_type* permutation)
440  : Epetra_GIDTypeVector<int_type>::impl(CV, map, permutation),
441  newObj_(NULL),
442  origObj_(NULL)
443 {
444  if (!isTypeSupported()) {
445  std::cerr << "unsupported type for permutation, aborting" << std::endl;
446  abort();
447  }
448 }
449 
450 template<typename T, typename int_type>
452  : Epetra_GIDTypeVector<int_type>::impl(map),
453  newObj_(NULL),
454  origObj_(NULL)
455 {
456  if (!isTypeSupported()) {
457  std::cerr << "unsupported type for permutation, aborting" << std::endl;
458  abort();
459  }
460 }
461 
462 template<typename T, typename int_type>
464  : Epetra_GIDTypeVector<int_type>::impl((const typename Epetra_GIDTypeVector<int_type>::impl&)src),
465  newObj_(NULL),
466  origObj_(NULL)
467 {
468  if (!isTypeSupported()) {
469  std::cerr << "unsupported type for permutation, aborting" << std::endl;
470  abort();
471  }
472 }
473 
474 template<typename T, typename int_type>
476 {
477  if (newObj_ != NULL) delete newObj_;
478 }
479 
480 template<typename T, typename int_type>
482 {
483  const char* type_name = Perm_traits<T>::typeName();
484  if (!strcmp(type_name, "unknown")) {
485  return(false);
486  }
487 
488  return( true );
489 }
490 
491 template<typename T, typename int_type>
494 {
495  //In this function we're going to produce a new object which is a
496  //row-permutation of the input object (orig).
497  //
498  //Our permutation inherits IntVector, and the permutation is defined by the
499  //contents of the integer vector 'p', such that if p[i] = j then row i of
500  //the input object becomes row j of the permuted object.
501  //
502  //The permutation is accomplished by creating a map defined by the
503  //permutation, then using an Epetra_Export operation to move data from the
504  //input object into the permuted object.
505  //
506  //The permutation may be global. In other words, the rows of the object may
507  //be arbitrarily rearranged, including across processors.
508  //
509 
510  origObj_ = &orig;
511 
512  //The 'Map()' accessor returns Epetra_DistObject::Map() for CrsGraph and
513  //CrsMatrix, which turns out to be the RowMap() for those objects. For
514  //MultiVector it returns the correct object because MultiVectors only have
515  //one map.
516 
517  const Epetra_BlockMap& origMap = orig.Map();
518 
519  //Create an Epetra_Map representing the permutation.
520 
521  Epetra_Map* pmap = new Epetra_Map((int_type) Epetra_DistObject::Map().NumGlobalPoints64(),
522  Epetra_DistObject::Map().NumMyPoints(),
524  (int_type) Epetra_DistObject::Map().IndexBase64(),
525  Epetra_DistObject::Map().Comm());
526 
527  TPermutation* p = this;
528 
529  //Next check that the maps are compatible. If they aren't, we'll redistribute
530  //the permutation to match the distribution of the input object.
531 
532  if (!pmap->PointSameAs(origMap)) {
533  Epetra_Export p_exporter(Epetra_DistObject::Map(), origMap);
534  TPermutation* newp = new TPermutation(origMap);
535  newp->Export(*p, p_exporter, Add);
536  p = newp;
537 
538  delete pmap;
539  pmap = new Epetra_Map((int_type) p->Map().NumGlobalPoints64(),
540  p->Map().NumMyPoints(),
541  p->Values(),
542  (int_type) p->Map().IndexBase64(),
543  p->Map().Comm());
544  }
545 
546  //Create the new object, initially giving it the map defined by the
547  //permutation.
548 
549  newObj_ = Perm_traits<T>::clone(origObj_, Copy, *pmap, 1);
550 
551  //Create an exporter which will export data from the original object to the
552  //permuted object.
553 
554  Epetra_Export exporter(origMap, *pmap);
555 
556  //Now export the original object to the permuted object.
557 
558  newObj_->Export(orig, exporter, Add);
559 
560  //Now, since the export operation moved not only row-contents but also
561  //row-numbering, we need to replace the permuted row-numbering with the
562  //original row-numbering. We do this by replacing the permuted map with
563  //the original row-map.
564 
566 
567  delete pmap;
568 
569  if (p != this) {
570  delete p; //delete "newp" created if the PointSameAs test failed above
571  }
572 
573  return( *newObj_ );
574 }
575 
576 template<typename T, typename int_type>
577 typename TPermutation<T, int_type>::OutputRef
578 TPermutation<T, int_type>::operator()( typename TPermutation<T, int_type>::InputRef orig,
579  bool column_permutation )
580 {
581  origObj_ = &orig;
582  newObj_ = NULL;
583 
584  if (!column_permutation) {
585  return( operator()(orig) );
586  }
587 
588  if (strcmp("Epetra_CrsMatrix", Perm_traits<T>::typeName()) &&
589  strcmp("Epetra_CrsGraph", Perm_traits<T>::typeName())) {
590  std::cerr << "Permutation: column-permutation only implemented for"
591  << "CrsMatrix and CrsGraph." << std::endl;
592  assert(0);
593  }
594 
596 
597  return( *newObj_ );
598 }
599 
600 } // namespace EpetraExt
601 
602 #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