Teko_BlockingEpetra.cpp

/*
// @HEADER
// 
// ***********************************************************************
// 
//      Teko: A package for block and physics based preconditioning
//                  Copyright 2010 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 Eric C. Cyr (eccyr@sandia.gov)
// 
// ***********************************************************************
// 
// @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 & gColMap = *subMaps[j].first;
   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 & gColMap = *subMaps[j].first;
   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;

   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]);
   }
}

} // end Blocking
} // end Epetra
} // end Teko