Teko_BlockingEpetra.cpp

// @HEADER
// *****************************************************************************
//      Teko: A package for block and physics based preconditioning
//
// Copyright 2010 NTESS and the Teko contributors.
// SPDX-License-Identifier: BSD-3-Clause
// *****************************************************************************
// @HEADER

#include "Teko_BlockingEpetra.hpp"

#include "Epetra_IntVector.h"
#include "Epetra_LocalMap.h"

using Teuchos::RCP;
using Teuchos::rcp;

namespace Teko {
namespace Epetra {
namespace Blocking {

const MapPair buildSubMap(const std::vector<int> &gid, const Epetra_Comm &comm) {
  Teuchos::RCP<Epetra_Map> gidMap    = rcp(new Epetra_Map(-1, gid.size(), &gid[0], 0, comm));
  Teuchos::RCP<Epetra_Map> contigMap = rcp(new Epetra_Map(-1, gid.size(), 0, comm));

  return std::make_pair(gidMap, contigMap);
}

const ImExPair buildExportImport(const Epetra_Map &baseMap, const MapPair &maps) {
  return std::make_pair(rcp(new Epetra_Import(*maps.first, baseMap)),
                        rcp(new Epetra_Export(*maps.first, baseMap)));
}

void buildSubVectors(const std::vector<MapPair> &maps,
                     std::vector<RCP<Epetra_MultiVector> > &vectors, int count) {
  std::vector<MapPair>::const_iterator mapItr;

  // loop over all maps
  for (mapItr = maps.begin(); mapItr != maps.end(); mapItr++) {
    // add new elements to vectors
    RCP<Epetra_MultiVector> mv = rcp(new Epetra_MultiVector(*(*mapItr).second, count));
    vectors.push_back(mv);
  }
}

void one2many(std::vector<RCP<Epetra_MultiVector> > &many, const Epetra_MultiVector &single,
              const std::vector<RCP<Epetra_Import> > &subImport) {
  // std::vector<RCP<Epetra_Vector> >::const_iterator vecItr;
  std::vector<RCP<Epetra_MultiVector> >::const_iterator vecItr;
  std::vector<RCP<Epetra_Import> >::const_iterator impItr;

  // using Importers fill the sub vectors from the mama vector
  for (vecItr = many.begin(), impItr = subImport.begin(); vecItr != many.end();
       ++vecItr, ++impItr) {
    // for ease of access to the destination
    RCP<Epetra_MultiVector> destVec = *vecItr;

    // extract the map with global indicies from the current vector
    const Epetra_BlockMap &globalMap = (*impItr)->TargetMap();

    // build the import vector as a view on the destination
    Epetra_MultiVector importVector(View, globalMap, destVec->Values(), destVec->Stride(),
                                    destVec->NumVectors());

    // perform the import
    importVector.Import(single, **impItr, Insert);
  }
}

void many2one(Epetra_MultiVector &one, const std::vector<RCP<const Epetra_MultiVector> > &many,
              const std::vector<RCP<Epetra_Export> > &subExport) {
  // std::vector<RCP<const Epetra_Vector> >::const_iterator vecItr;
  std::vector<RCP<const Epetra_MultiVector> >::const_iterator vecItr;
  std::vector<RCP<Epetra_Export> >::const_iterator expItr;

  // using Exporters fill the empty vector from the sub-vectors
  for (vecItr = many.begin(), expItr = subExport.begin(); vecItr != many.end();
       ++vecItr, ++expItr) {
    // for ease of access to the source
    RCP<const Epetra_MultiVector> srcVec = *vecItr;

    // extract the map with global indicies from the current vector
    const Epetra_BlockMap &globalMap = (*expItr)->SourceMap();

    // build the export vector as a view of the destination
    Epetra_MultiVector exportVector(View, globalMap, srcVec->Values(), srcVec->Stride(),
                                    srcVec->NumVectors());
    one.Export(exportVector, **expItr, Insert);
  }
}

RCP<Epetra_IntVector> getSubBlockColumnGIDs(const Epetra_CrsMatrix &A, const MapPair &mapPair) {
  RCP<const Epetra_Map> blkGIDMap    = mapPair.first;
  RCP<const Epetra_Map> blkContigMap = mapPair.second;

  // fill index vector for rows
  Epetra_IntVector rIndex(A.RowMap(), true);
  for (int i = 0; i < A.NumMyRows(); i++) {
    // LID is need to get to contiguous map
    int lid = blkGIDMap->LID(A.GRID(i));  // this LID makes me nervous
    if (lid > -1)
      rIndex[i] = blkContigMap->GID(lid);
    else
      rIndex[i] = -1;
  }

  // get relavant column indices
  Epetra_Import import(A.ColMap(), A.RowMap());
  RCP<Epetra_IntVector> cIndex = rcp(new Epetra_IntVector(A.ColMap(), true));
  cIndex->Import(rIndex, import, Insert);

  return cIndex;
}

// build a single subblock Epetra_CrsMatrix
RCP<Epetra_CrsMatrix> buildSubBlock(int i, int j, const Epetra_CrsMatrix &A,
                                    const std::vector<MapPair> &subMaps) {
  // get the number of variables families
  int numVarFamily = subMaps.size();

  TEUCHOS_ASSERT(i >= 0 && i < numVarFamily);
  TEUCHOS_ASSERT(j >= 0 && j < numVarFamily);

  const Epetra_Map &gRowMap = *subMaps[i].first;   // new GIDs
  const Epetra_Map &rowMap  = *subMaps[i].second;  // contiguous GIDs
  const Epetra_Map &colMap  = *subMaps[j].second;

  const RCP<Epetra_IntVector> plocal2ContigGIDs = getSubBlockColumnGIDs(A, subMaps[j]);
  Epetra_IntVector &local2ContigGIDs            = *plocal2ContigGIDs;

  RCP<Epetra_CrsMatrix> mat = rcp(new Epetra_CrsMatrix(Copy, rowMap, 0));

  // get entry information
  int numMyRows     = rowMap.NumMyElements();
  int maxNumEntries = A.MaxNumEntries();

  // for extraction
  std::vector<int> indices(maxNumEntries);
  std::vector<double> values(maxNumEntries);

  // for insertion
  std::vector<int> colIndices(maxNumEntries);
  std::vector<double> colValues(maxNumEntries);

  // insert each row into subblock
  // let FillComplete handle column distribution
  for (int localRow = 0; localRow < numMyRows; localRow++) {
    int numEntries = -1;
    int globalRow  = gRowMap.GID(localRow);
    int lid        = A.RowMap().LID(globalRow);
    int contigRow  = rowMap.GID(localRow);
    TEUCHOS_ASSERT(lid > -1);

    int err = A.ExtractMyRowCopy(lid, maxNumEntries, numEntries, &values[0], &indices[0]);
    TEUCHOS_ASSERT(err == 0);

    int numOwnedCols = 0;
    for (int localCol = 0; localCol < numEntries; localCol++) {
      TEUCHOS_ASSERT(indices[localCol] > -1);

      // if global id is not owned by this column
      int gid = local2ContigGIDs[indices[localCol]];
      if (gid == -1) continue;  // in contiguous row

      colIndices[numOwnedCols] = gid;
      colValues[numOwnedCols]  = values[localCol];
      numOwnedCols++;
    }

    // insert it into the new matrix
    mat->InsertGlobalValues(contigRow, numOwnedCols, &colValues[0], &colIndices[0]);
  }

  // fill it and automagically optimize the storage
  mat->FillComplete(colMap, rowMap);

  return mat;
}

// build a single subblock Epetra_CrsMatrix
void rebuildSubBlock(int i, int j, const Epetra_CrsMatrix &A, const std::vector<MapPair> &subMaps,
                     Epetra_CrsMatrix &mat) {
  // get the number of variables families
  int numVarFamily = subMaps.size();

  TEUCHOS_ASSERT(i >= 0 && i < numVarFamily);
  TEUCHOS_ASSERT(j >= 0 && j < numVarFamily);

  const Epetra_Map &gRowMap = *subMaps[i].first;   // new GIDs
  const Epetra_Map &rowMap  = *subMaps[i].second;  // contiguous GIDs

  const RCP<Epetra_IntVector> plocal2ContigGIDs = getSubBlockColumnGIDs(A, subMaps[j]);
  Epetra_IntVector &local2ContigGIDs            = *plocal2ContigGIDs;

  mat.PutScalar(0.0);

  // get entry information
  int numMyRows     = rowMap.NumMyElements();
  int maxNumEntries = A.MaxNumEntries();

  // for extraction
  std::vector<int> indices(maxNumEntries);
  std::vector<double> values(maxNumEntries);

  // for insertion
  std::vector<int> colIndices(maxNumEntries);
  std::vector<double> colValues(maxNumEntries);

  // insert each row into subblock
  // let FillComplete handle column distribution
  for (int localRow = 0; localRow < numMyRows; localRow++) {
    int numEntries = -1;
    int globalRow  = gRowMap.GID(localRow);
    int lid        = A.RowMap().LID(globalRow);
    int contigRow  = rowMap.GID(localRow);
    TEUCHOS_ASSERT(lid > -1);

    int err = A.ExtractMyRowCopy(lid, maxNumEntries, numEntries, &values[0], &indices[0]);
    TEUCHOS_ASSERT(err == 0);

    int numOwnedCols = 0;
    for (int localCol = 0; localCol < numEntries; localCol++) {
      TEUCHOS_ASSERT(indices[localCol] > -1);

      // if global id is not owned by this column
      int gid = local2ContigGIDs[indices[localCol]];
      if (gid == -1) continue;  // in contiguous row

      colIndices[numOwnedCols] = gid;
      colValues[numOwnedCols]  = values[localCol];
      numOwnedCols++;
    }

    // insert it into the new matrix
    mat.SumIntoGlobalValues(contigRow, numOwnedCols, &colValues[0], &colIndices[0]);
  }
}

}  // namespace Blocking
}  // namespace Epetra
}  // namespace Teko