EpetraExt_Permutation_impl.h

Go to the documentation of this file.

//@HEADER
// ***********************************************************************
//
//     EpetraExt: Epetra Extended - Linear Algebra Services Package
//                 Copyright (2011) Sandia Corporation
//
// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
// the U.S. Government retains certain rights in this software.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are
// met:
//
// 1. Redistributions of source code must retain the above copyright
// notice, this list of conditions and the following disclaimer.
//
// 2. Redistributions in binary form must reproduce the above copyright
// notice, this list of conditions and the following disclaimer in the
// documentation and/or other materials provided with the distribution.
//
// 3. Neither the name of the Corporation nor the names of the
// contributors may be used to endorse or promote products derived from
// this software without specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
//
// ***********************************************************************
//@HEADER
#ifndef EpetraExt_PERMUTATION_IMPL_H
#define EpetraExt_PERMUTATION_IMPL_H

#include <EpetraExt_ConfigDefs.h>

#include <EpetraExt_Permutation.h>

#include <Epetra_Export.h>
#include <Epetra_Map.h>
#include <Epetra_Comm.h>
#include <Epetra_MultiVector.h>
#include <Epetra_CrsGraph.h>
#include <Epetra_CrsMatrix.h>

namespace EpetraExt {

template<class T>
struct Perm_traits {
  static const char* typeName()
  { static const char name[] = "unknown"; return( name ); }

  static T* clone(T* example,
      Epetra_DataAccess CV,
      const Epetra_BlockMap& map,
      int int_argument)
  {  return( NULL ); }

  static void replaceMap(T* obj, const Epetra_BlockMap& map)
  { cerr << "not implemented for unknown type"<<endl; }

  static T*
  produceColumnPermutation(Permutation<T>* perm,
         T* srcObj)
  { cerr << "not implemented for unknown type"<<endl; }

};//struct Perm_traits



template<>
struct Perm_traits<Epetra_CrsMatrix> {

  static const char* typeName()
  { static const char name[] = "Epetra_CrsMatrix"; return( name ); }


  static Epetra_CrsMatrix* clone(Epetra_CrsMatrix* example,
         Epetra_DataAccess CV,
         const Epetra_BlockMap& map,
         int rowLength)
  {
    //don't need the example object currently...
    (void)example;

    //we need a Epetra_Map, rather than a Epetra_BlockMap, to create a
    //Epetra_CrsMatrix.

    const Epetra_Map* pointmap =
      dynamic_cast<const Epetra_Map*>(&map);
    if (pointmap == NULL) {
      cerr << "dynamic_cast<const Epetra_Map*> failed."<<endl;
      return(NULL);
    }

    return( new Epetra_CrsMatrix(CV, *pointmap, rowLength) );
  }


  static void replaceMap(Epetra_CrsMatrix* mat, const Epetra_BlockMap& map)
  { mat->ReplaceRowMap(map); }

  static Epetra_CrsMatrix*
  produceColumnPermutation(Permutation<Epetra_CrsMatrix>* perm,
         Epetra_CrsMatrix* srcObj)
  {
    //First we need to export this permutation to match the column-map of the
    //object being column-permuted. (We need to have locally available all
    //elements of the permutation corresponding to the local columns of the
    //object being permuted.)

    const Epetra_Map& origColMap = srcObj->ColMap();

    Permutation<Epetra_CrsMatrix>* colperm =
      new Permutation<Epetra_CrsMatrix>(origColMap);
    colperm->PutValue(0);

    Epetra_Export p_exporter(perm->Map(), origColMap);
    colperm->Export(*perm, p_exporter, Add);

    const Epetra_Map& origRowMap = srcObj->RowMap();
    int numMyRows = origRowMap.NumMyElements();
    const int* myGlobalRows = origRowMap.MyGlobalElements();

    //Create the new object, giving it the same map as the original object.

    Epetra_CrsMatrix* result = new Epetra_CrsMatrix(Copy, origRowMap, 1);

    for(int i=0; i<numMyRows; ++i) {
      int globalRow = myGlobalRows[i];
      int len = srcObj->NumGlobalEntries(globalRow);

      int numIndices;
      double* src_values = new double[len];
      int* src_indices = new int[len];
      int err = srcObj->ExtractGlobalRowCopy(globalRow, len, numIndices,
               src_values, src_indices);
      if (err < 0 || numIndices != len) {
  cerr<<"Perm_traits<CrsMatrix>::produceColumnPermutation err("<<err<<") row "
      <<globalRow<<", len "<<len<<", numIndices "<<numIndices<<endl;
      }

      int* pindices = new int[len];

      const Epetra_BlockMap& pmap = colperm->Map();
      int* p = colperm->Values();

      for(int j=0; j<len; ++j) {
  int old_col = src_indices[j];

  int lid = pmap.LID(old_col);
  if (lid<0) {
    cerr << "Perm_traits<CrsMatrix>::permuteColumnIndices GID("<<old_col
         <<") not found"<<endl;
    break;
  }

  pindices[j] = p[lid];
      }

      err = result->InsertGlobalValues(globalRow, len, src_values, pindices);
      if (err < 0) {
  cerr << "Perm_traits<CrsMatrix>::permuteColumnIndices err("<<err
       <<") row "<<globalRow<<endl;
      }

      delete [] pindices;
      delete [] src_indices;
      delete [] src_values;
    }

    result->FillComplete();

    delete colperm;

    return(result);
  }

};//struct Perm_traits<Epetra_CrsMatrix>



template<>
struct Perm_traits<Epetra_CrsGraph> {

  static const char* typeName()
  { static const char name[] = "Epetra_CrsGraph"; return( name ); }


  static Epetra_CrsGraph* clone(Epetra_CrsGraph* example,
        Epetra_DataAccess CV,
        const Epetra_BlockMap& map,
        int rowLength)
  {
    //don't need the example object currently...
    (void)example;

    return( new Epetra_CrsGraph(CV, map, rowLength) );
  }


  static void replaceMap(Epetra_CrsGraph* graph, const Epetra_BlockMap& map)
  { graph->ReplaceRowMap(map); }

  static Epetra_CrsGraph*
  produceColumnPermutation(Permutation<Epetra_CrsGraph>* perm,
         Epetra_CrsGraph* srcObj)
  {
    //First we need to export this permutation to match the column-map of the
    //object being column-permuted. (We need to have locally available all
    //elements of the permutation corresponding to the local columns of the
    //object being permuted.)

    const Epetra_BlockMap& origColMap = srcObj->ColMap();

    Permutation<Epetra_CrsGraph>* colperm =
      new Permutation<Epetra_CrsGraph>(origColMap);
    colperm->PutValue(0);

    Epetra_Export p_exporter(perm->Map(), origColMap);
    colperm->Export(*perm, p_exporter, Add);

    const Epetra_BlockMap& origRowMap = srcObj->RowMap();
    int numMyRows = origRowMap.NumMyElements();
    const int* myGlobalRows = origRowMap.MyGlobalElements();

    //Create the new object, giving it the same map as the original object.

    Epetra_CrsGraph* result = new Epetra_CrsGraph(Copy, origRowMap, 1);

    for(int i=0; i<numMyRows; ++i) {
      int globalRow = myGlobalRows[i];
      int len = srcObj->NumGlobalIndices(globalRow);

      int numIndices;
      int* src_indices = new int[len];
      int err = srcObj->ExtractGlobalRowCopy(globalRow, len, numIndices, src_indices);
      if (err < 0 || numIndices != len) {
  cerr<<"Perm_traits<CrsGraph>::produceColumnPermutation err("<<err<<") row "
    <<globalRow<<", len "<<len<<", numIndices "<<numIndices<<endl;
      }

      int* pindices = new int[len];

      const Epetra_BlockMap& pmap = colperm->Map();
      int* p = colperm->Values();

      for(int j=0; j<len; ++j) {
  int old_col = src_indices[j];

  int lid = pmap.LID(old_col);
  if (lid<0) {
    cerr << "Perm_traits<CrsGraph>::permuteColumnIndices GID("<<old_col
         <<") not found"<<endl;
    break;
  }

  pindices[j] = p[lid];
      }

      err = result->InsertGlobalIndices(globalRow, len, pindices);
      if (err < 0) {
  cerr << "Perm_traits<CrsGraph>::produceColumnPermutation err("<<err
       <<") row "<<globalRow<<endl;
      }

      delete [] pindices;
      delete [] src_indices;
    }

    result->FillComplete();

    delete colperm;

    return(result);
  }

};//struct Perm_traits<Epetra_CrsGraph>


template<>
struct Perm_traits<Epetra_MultiVector> {

  static const char* typeName()
  { static const char name[] = "Epetra_MultiVector"; return( name ); }


  static Epetra_MultiVector* clone(Epetra_MultiVector* example,
           Epetra_DataAccess CV,
           const Epetra_BlockMap& map,
           int numVectors)
  {
    return( new Epetra_MultiVector(map, example->NumVectors()) );
  }


  static void replaceMap(Epetra_MultiVector* mvec, const Epetra_BlockMap& map)
  { mvec->ReplaceMap(map); }

  static Epetra_MultiVector*
  produceColumnPermutation(Permutation<Epetra_MultiVector>* perm,
         Epetra_MultiVector* srcObj)
  {
    cerr << "col-permutation not implemented for Epetra_MultiVector"<<endl;
    return(NULL);
  }

};//struct Perm_traits<Epetra_CrsGraph>


//-------------------------------------------------------------------------
//Now the method definitions for the EpetraExt::Permutation class.
//-------------------------------------------------------------------------

template<typename T>
Permutation<T>::Permutation(Epetra_DataAccess CV,
                         const Epetra_BlockMap& map,
                         int* permutation)
  : Epetra_IntVector(CV, map, permutation),
    newObj_(NULL),
    origObj_(NULL)
{
  if (!isTypeSupported()) {
    cerr << "unsupported type for permutation, aborting" << endl;
    abort();
  }
}

template<typename T>
Permutation<T>::Permutation(const Epetra_BlockMap& map)
  : Epetra_IntVector(map),
    newObj_(NULL),
    origObj_(NULL)
{
  if (!isTypeSupported()) {
    cerr << "unsupported type for permutation, aborting" << endl;
    abort();
  }
}

template<typename T>
Permutation<T>::Permutation(const Permutation& src)
  : Epetra_IntVector((const Epetra_IntVector&)src),
    newObj_(NULL),
    origObj_(NULL)
{
  if (!isTypeSupported()) {
    cerr << "unsupported type for permutation, aborting" << endl;
    abort();
  }
}

template<typename T>
Permutation<T>::~Permutation()
{
  if (newObj_ != NULL) delete newObj_;
}

template<typename T>
bool Permutation<T>::isTypeSupported()
{
  const char* type_name = Perm_traits<T>::typeName();
  if (!strcmp(type_name, "unknown")) {
    return(false);
  }

  return( true );
}

template<typename T>
typename Permutation<T>::OutputRef
Permutation<T>::operator()( typename Permutation<T>::InputRef orig )
{
  //In this function we're going to produce a new object which is a
  //row-permutation of the input object (orig).
  //
  //Our permutation inherits IntVector, and the permutation is defined by the
  //contents of the integer vector 'p', such that if p[i] = j then row i of
  //the input object becomes row j of the permuted object.
  //
  //The permutation is accomplished by creating a map defined by the
  //permutation, then using an Epetra_Export operation to move data from the
  //input object into the permuted object.
  //
  //The permutation may be global. In other words, the rows of the object may
  //be arbitrarily rearranged, including across processors.
  //

  origObj_ = &orig;

  //The 'Map()' accessor returns Epetra_DistObject::Map() for CrsGraph and
  //CrsMatrix, which turns out to be the RowMap() for those objects. For
  //MultiVector it returns the correct object because MultiVectors only have
  //one map.

  const Epetra_BlockMap& origMap = orig.Map();

  //Create an Epetra_Map representing the permutation.

  Epetra_Map* pmap = new Epetra_Map(Map().NumGlobalPoints(),
            Map().NumMyPoints(),
            Values(),
            Map().IndexBase(),
            Map().Comm());

  Permutation* p = this;

  //Next check that the maps are compatible. If they aren't, we'll redistribute
  //the permutation to match the distribution of the input object.

  if (!pmap->PointSameAs(origMap)) {
    Epetra_Export p_exporter(Map(), origMap);
    Permutation* newp = new Permutation(origMap);
    newp->Export(*p, p_exporter, Add);
    p = newp;

    delete pmap;
    pmap = new Epetra_Map(p->Map().NumGlobalPoints(),
        p->Map().NumMyPoints(),
        p->Values(),
        p->Map().IndexBase(),
        p->Map().Comm());
  }

  //Create the new object, initially giving it the map defined by the
  //permutation.

  newObj_ = Perm_traits<T>::clone(origObj_, Copy, *pmap, 1);

  //Create an exporter which will export data from the original object to the
  //permuted object.

  Epetra_Export exporter(origMap, *pmap);

  //Now export the original object to the permuted object.

  newObj_->Export(orig, exporter, Add);

  //Now, since the export operation moved not only row-contents but also
  //row-numbering, we need to replace the permuted row-numbering with the
  //original row-numbering. We do this by replacing the permuted map with
  //the original row-map.

  Perm_traits<T>::replaceMap(newObj_, origMap);

  delete pmap;

  if (p != this) {
    delete p; //delete "newp" created if the PointSameAs test failed above
  }

  return( *newObj_ );
}

template<typename T>
typename Permutation<T>::OutputRef
Permutation<T>::operator()( typename Permutation<T>::InputRef orig,
          bool column_permutation )
{
  origObj_ = &orig;
  newObj_ = NULL;

  if (!column_permutation) {
    return( operator()(orig) );
  }

  if (strcmp("Epetra_CrsMatrix", Perm_traits<T>::typeName()) &&
      strcmp("Epetra_CrsGraph", Perm_traits<T>::typeName())) {
    cerr << "Permutation: column-permutation only implemented for"
   << "CrsMatrix and CrsGraph." << endl;
    assert(0);
  }

  newObj_ = Perm_traits<T>::produceColumnPermutation(this, &orig);

  return( *newObj_ );
}

} // namespace EpetraExt

#endif //EpetraExt_PERMUTATION_IMPL_H