EpetraExt Development
Loading...
Searching...
No Matches
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
45
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
56namespace EpetraExt {
57
76template<class T>
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*
104 produceColumnPermutation(TPermutation<T, int_type>* perm,
105 T* srcObj)
106 { std::cerr << "not implemented for unknown type"<<std::endl; }
107
108};//struct Perm_traits
109
110
111
115template<>
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;
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);
141
142 return( new Epetra_CrsMatrix(CV, *pointmap, rowLength) );
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*
153 TproduceColumnPermutation(TPermutation<Epetra_CrsMatrix, int_type>* perm,
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
163 TPermutation<Epetra_CrsMatrix, int_type>* colperm =
164 new TPermutation<Epetra_CrsMatrix, int_type>(origColMap);
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
231 static Epetra_CrsMatrix*
232 produceColumnPermutation(TPermutation<Epetra_CrsMatrix, int>* perm,
233 Epetra_CrsMatrix* srcObj)
234 {
235 return TproduceColumnPermutation<int>(perm, srcObj);
236 }
237#endif
238
239#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
241 static Epetra_CrsMatrix*
242 produceColumnPermutation(TPermutation<Epetra_CrsMatrix, long long>* perm,
243 Epetra_CrsMatrix* srcObj)
244 {
245 return TproduceColumnPermutation<long long>(perm, srcObj);
246 }
247#endif
248};//struct Perm_traits<Epetra_CrsMatrix>
249
250
251
255template<>
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*
283 TproduceColumnPermutation(TPermutation<Epetra_CrsGraph, int_type>* perm,
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
293 TPermutation<Epetra_CrsGraph, int_type>* colperm =
294 new TPermutation<Epetra_CrsGraph, int_type>(origColMap);
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
358 static Epetra_CrsGraph*
359 produceColumnPermutation(TPermutation<Epetra_CrsGraph, int>* perm,
360 Epetra_CrsGraph* srcObj)
361 {
362 return TproduceColumnPermutation<int>(perm, srcObj);
363 }
364#endif
365
366#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
368 static Epetra_CrsGraph*
369 produceColumnPermutation(TPermutation<Epetra_CrsGraph, long long>* perm,
370 Epetra_CrsGraph* srcObj)
371 {
372 return TproduceColumnPermutation<long long>(perm, srcObj);
373 }
374#endif
375};//struct Perm_traits<Epetra_CrsGraph>
376
377
381template<>
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
405 static Epetra_MultiVector*
406 produceColumnPermutation(Permutation<Epetra_MultiVector>* /* perm */,
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
415 static Epetra_MultiVector*
416 produceColumnPermutation(Permutation64<Epetra_MultiVector>* /* perm */,
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
430template<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
444template<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
456template<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
468template<typename T, typename int_type>
470{
471 if (newObj_ != NULL) delete newObj_;
472}
473
474template<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
485template<typename T, typename int_type>
486typename TPermutation<T, int_type>::OutputRef
487TPermutation<T, int_type>::operator()( typename TPermutation<T, int_type>::InputRef orig )
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
559 Perm_traits<T>::replaceMap(newObj_, origMap);
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
570template<typename T, typename int_type>
571typename TPermutation<T, int_type>::OutputRef
572TPermutation<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
589 newObj_ = Perm_traits<T>::produceColumnPermutation(this, &orig);
590
591 return( *newObj_ );
592}
593
594} // namespace EpetraExt
595
596#endif //EpetraExt_PERMUTATION_IMPL_H
Add
Epetra_DataAccess
Copy
Permutation stores and describes a permutation matrix P.
TPermutation(Epetra_DataAccess CV, const Epetra_BlockMap &map, int_type *permutation)
Constructor.
int LID(int GID) const
int NumMyElements() const
const Epetra_BlockMap & Map() const
EpetraExt::BlockCrsMatrix: A class for constructing a distributed block matrix.
static void replaceMap(Epetra_CrsGraph *graph, const Epetra_BlockMap &map)
replaceMap implementation
static Epetra_CrsGraph * produceColumnPermutation(TPermutation< Epetra_CrsGraph, long long > *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 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 * produceColumnPermutation(TPermutation< Epetra_CrsGraph, int > *perm, Epetra_CrsGraph *srcObj)
return new object which is a column-permutation of srcObj
static const char * typeName()
typeName implementation
static Epetra_CrsMatrix * produceColumnPermutation(TPermutation< Epetra_CrsMatrix, long long > *perm, Epetra_CrsMatrix *srcObj)
return new object, which is a column-permutation of srcObj
static Epetra_CrsMatrix * TproduceColumnPermutation(TPermutation< Epetra_CrsMatrix, int_type > *perm, Epetra_CrsMatrix *srcObj)
return new object, which is a column-permutation of srcObj
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
static Epetra_CrsMatrix * produceColumnPermutation(TPermutation< Epetra_CrsMatrix, int > *perm, Epetra_CrsMatrix *srcObj)
return new object, which is a column-permutation of srcObj
static const char * typeName()
typeName implementation
static void replaceMap(Epetra_MultiVector *mvec, const Epetra_BlockMap &map)
replaceMap implementation
static const char * typeName()
typeName implementation
static Epetra_MultiVector * clone(Epetra_MultiVector *example, Epetra_DataAccess, const Epetra_BlockMap &map, int)
clone implementation
static Epetra_MultiVector * produceColumnPermutation(Permutation64< Epetra_MultiVector > *, Epetra_MultiVector *)
permute column-indices within a specified row, if applicable
static Epetra_MultiVector * produceColumnPermutation(Permutation< Epetra_MultiVector > *, Epetra_MultiVector *)
permute column-indices within a specified row, if applicable
Define some traits to make it easier to deal with template-parameters which are objects to be permute...
static void replaceMap(T *obj, const Epetra_BlockMap &map)
replace the object's row-map (or if it's not a matrix, replace its only map)
static T * produceColumnPermutation(TPermutation< T, int_type > *perm, T *srcObj)
return new object, which is a column-permutation of srcObj
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 ...
static const char * typeName()
return a std::string name for the object type