Tpetra_Import_Util.hpp

Go to the documentation of this file.

// @HEADER
// *****************************************************************************
//          Tpetra: Templated Linear Algebra Services Package
//
// Copyright 2008 NTESS and the Tpetra contributors.
// SPDX-License-Identifier: BSD-3-Clause
// *****************************************************************************
// @HEADER
 
#ifndef TPETRA_IMPORT_UTIL_HPP
#define TPETRA_IMPORT_UTIL_HPP

 
#include "Tpetra_ConfigDefs.hpp"
#include "Tpetra_Import.hpp"
#include "Tpetra_HashTable.hpp"
#include "Tpetra_Map.hpp"
#include "Tpetra_Util.hpp"
#include "Tpetra_Distributor.hpp"
#include <Teuchos_Array.hpp>
#include <utility>
 
namespace Tpetra {
  namespace Import_Util {
    template <typename LocalOrdinal, typename GlobalOrdinal, typename Node>
    void
    getPidGidPairs (const Tpetra::Import<LocalOrdinal,GlobalOrdinal,Node>& Importer,
                    Teuchos::Array< std::pair<int,GlobalOrdinal> >& gpids,
                    bool use_minus_one_for_local);

    template <typename LocalOrdinal, typename GlobalOrdinal, typename Node>
    void
    getPids (const Tpetra::Import<LocalOrdinal,GlobalOrdinal,Node>& Importer,
             Teuchos::Array<int>& pids,
             bool use_minus_one_for_local);

    // Like the above, but without the resize
    template <typename LocalOrdinal, typename GlobalOrdinal, typename Node>
    void
    getPids (const Tpetra::Import<LocalOrdinal,GlobalOrdinal,Node>& Importer,
             Teuchos::ArrayView<int>& pids,
             bool use_minus_one_for_local);
 

    template <typename LocalOrdinal, typename GlobalOrdinal, typename Node>
    void
    getRemotePIDs (const Tpetra::Import<LocalOrdinal,GlobalOrdinal,Node>& Importer,
                   Teuchos::Array<int>& RemotePIDs);
 
  } // namespace Import_Util
} // namespace Tpetra
 
namespace Tpetra {
namespace Import_Util {
 
template <typename LocalOrdinal, typename GlobalOrdinal, typename Node>
void
getPidGidPairs (const Tpetra::Import<LocalOrdinal,GlobalOrdinal,Node>& Importer,
                Teuchos::Array< std::pair<int,GlobalOrdinal> >& gpids,
                bool use_minus_one_for_local)
{
  // Put the (PID,GID) pair in member of Importer.TargetMap() in
  // gpids.  If use_minus_one_for_local==true, put in -1 instead of
  // MyPID.
  const Tpetra::Distributor& D = Importer.getDistributor();
 
  LocalOrdinal ii;
  size_t  i,j,k;
  int mypid = Importer.getTargetMap()->getComm()->getRank();
  size_t N  = Importer.getTargetMap()->getLocalNumElements();
 
  // Get the importer's data
  Teuchos::ArrayView<const LocalOrdinal> RemoteLIDs  = Importer.getRemoteLIDs();
 
  // Get the distributor's data
  size_t NumReceives                           = D.getNumReceives();
  Teuchos::ArrayView<const int> ProcsFrom      = D.getProcsFrom();
  Teuchos::ArrayView<const size_t> LengthsFrom = D.getLengthsFrom();
 
  // Resize the outgoing data structure
  gpids.resize(N);
 
  // Start by claiming that I own all the data
  LocalOrdinal lzero = Teuchos::ScalarTraits<LocalOrdinal>::zero();
  if(use_minus_one_for_local)
    for(ii=lzero; Teuchos::as<size_t>(ii)<N; ii++) gpids[ii]=std::make_pair(-1,Importer.getTargetMap()->getGlobalElement(ii));
  else
    for(ii=lzero; Teuchos::as<size_t>(ii)<N; ii++) gpids[ii]=std::make_pair(mypid,Importer.getTargetMap()->getGlobalElement(ii));
 
  // Now, for each remote ID, record who actually owns it.  This loop follows the operation order in the
  // MpiDistributor so it ought to duplicate that effect.
  for(i=0,j=0; i<NumReceives; i++){
    int pid=ProcsFrom[i];
    for(k=0; k<LengthsFrom[i]; k++){
      if(pid!=mypid) gpids[RemoteLIDs[j]].first=pid;
      j++;
    }
  }
}
 
template <typename LocalOrdinal, typename GlobalOrdinal, typename Node>
void
getPids (const Tpetra::Import<LocalOrdinal,GlobalOrdinal,Node>& Importer,
         Teuchos::Array<int>& pids,
         bool use_minus_one_for_local)
{
  // Resize the outgoing data structure
  pids.resize(Importer.getTargetMap()->getLocalNumElements());
  Teuchos::ArrayView<int> v_pids = pids();
  getPids(Importer,v_pids,use_minus_one_for_local);
}
 
 
template <typename LocalOrdinal, typename GlobalOrdinal, typename Node>
void
getPids (const Tpetra::Import<LocalOrdinal,GlobalOrdinal,Node>& Importer,
         Teuchos::ArrayView<int>& pids,
         bool use_minus_one_for_local)
{
  const Tpetra::Distributor & D=Importer.getDistributor();
 
  LocalOrdinal ii;
  size_t  i,j,k;
  int mypid = Importer.getTargetMap()->getComm()->getRank();
  size_t N  = Importer.getTargetMap()->getLocalNumElements();
  if(N!=(size_t)pids.size()) throw std::runtime_error("Tpetra::Import_Util::getPids(): Incorrect size for output array");
 
  // Get the importer's data
  Teuchos::ArrayView<const LocalOrdinal> RemoteLIDs  = Importer.getRemoteLIDs();
 
  // Get the distributor's data
  size_t NumReceives                           = D.getNumReceives();
  Teuchos::ArrayView<const int> ProcsFrom      = D.getProcsFrom();
  Teuchos::ArrayView<const size_t> LengthsFrom = D.getLengthsFrom();
 
  // Start by claiming that I own all the data
  LocalOrdinal lzero = Teuchos::ScalarTraits<LocalOrdinal>::zero();
  if(use_minus_one_for_local)
    for(ii=lzero; Teuchos::as<size_t>(ii)<N; ii++) pids[ii]=-1;
  else
    for(ii=lzero; Teuchos::as<size_t>(ii)<N; ii++) pids[ii]=mypid;
 
  // Now, for each remote ID, record who actually owns it.  This loop follows the operation order in the
  // MpiDistributor so it ought to duplicate that effect.
  for(i=0,j=0; i<NumReceives; i++){
    int pid=ProcsFrom[i];
    for(k=0; k<LengthsFrom[i]; k++){
      if(pid!=mypid) pids[RemoteLIDs[j]]=pid;
      j++;
    }
  }
}
 
template <typename LocalOrdinal, typename GlobalOrdinal, typename Node>
void
getRemotePIDs (const Tpetra::Import<LocalOrdinal,GlobalOrdinal,Node>& Importer,
               Teuchos::Array<int>& RemotePIDs)
{
  const Tpetra::Distributor& D = Importer.getDistributor();
 
  // Get the importer's data
  Teuchos::ArrayView<const LocalOrdinal> RemoteLIDs  = Importer.getRemoteLIDs();
 
  // Get the distributor's data
  size_t NumReceives                           = D.getNumReceives();
  Teuchos::ArrayView<const int> ProcsFrom      = D.getProcsFrom();
  Teuchos::ArrayView<const size_t> LengthsFrom = D.getLengthsFrom();
 
  // Resize the outgoing data structure
  RemotePIDs.resize(Importer.getNumRemoteIDs());
 
  // Now, for each remote ID, record who actually owns it.  This loop
  // follows the operation order in the MpiDistributor so it ought to
  // duplicate that effect.
  size_t i,j,k;
  for (i = 0, j = 0; i < NumReceives; ++i) {
    const int pid = ProcsFrom[i];
    for (k = 0; k < LengthsFrom[i]; ++k) {
      RemotePIDs[j] = pid;
      j++;
    }
  }
}
 
 
/* Check some of the validity of an Import object
   WARNING: This is a debugging routine only. */
template <typename LocalOrdinal, typename GlobalOrdinal, typename Node>
bool
checkImportValidity (const Tpetra::Import<LocalOrdinal,GlobalOrdinal,Node>& Importer)
{
  using Teuchos::RCP;
  RCP<const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > source = Importer.getSourceMap();
  RCP<const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > target = Importer.getTargetMap();
  RCP<const Teuchos::Comm<int> > comm = source->getComm();
 
  // For now, do not check validity of a locally replicated source map (just return true)
  if (!source->isDistributed()) return true;
 
  int global_is_valid=0;
  bool is_valid=true;
 
  // We check validity by going through each ID in the source map one by one, broadcasting the sender's PID and then all
  // receivers check it.
  LocalOrdinal LO_INVALID = Teuchos::OrdinalTraits<LocalOrdinal>::invalid();
  const int MyPID    = comm->getRank();
  const int NumProcs = comm->getSize();
 
  GlobalOrdinal minSourceGID = source->getMinAllGlobalIndex();
  GlobalOrdinal maxSourceGID = source->getMaxAllGlobalIndex();
  GlobalOrdinal minTargetGID = target->getMinAllGlobalIndex();
  GlobalOrdinal maxTargetGID = target->getMaxAllGlobalIndex();
 
  std::ostringstream os;
 
  /***********************************************/
  /*              Check recv side                */
  /***********************************************/
 
  Teuchos::ArrayView<const LocalOrdinal> permuteTarget = Importer.getPermuteToLIDs();
  Teuchos::ArrayView<const LocalOrdinal> remoteLIDs    = Importer.getRemoteLIDs();
  Teuchos::ArrayView<const LocalOrdinal> exportLIDs    = Importer.getExportLIDs();
  Teuchos::ArrayView<const LocalOrdinal> exportPIDs    = Importer.getExportPIDs();
  Teuchos::Array<int> remotePIDs; getRemotePIDs(Importer,remotePIDs);
 
  // Generate remoteGIDs
  Teuchos::Array<GlobalOrdinal> remoteGIDs(remoteLIDs.size());
  for(size_t i=0; i<(size_t)remoteLIDs.size(); i++) {
    remoteGIDs[i] = target->getGlobalElement(remoteLIDs[i]);
    if(remoteGIDs[i]<0) {
      os<<MyPID<<"ERROR3: source->getGlobalElement(remoteLIDs[l]) is invalid GID="<<remoteGIDs[i]<<" LID= "<<remoteLIDs[i]<<std::endl;
      is_valid=false;
    }
}
  // Generate exportGIDs
  Teuchos::Array<GlobalOrdinal> exportGIDs(exportLIDs.size(),-1);
  for(size_t i=0; i<(size_t)exportLIDs.size(); i++) {
    exportGIDs[i] = source->getGlobalElement(exportLIDs[i]);
    exportGIDs[i]=source->getGlobalElement(exportLIDs[i]);
    if(exportGIDs[i]<0) {
      os<<MyPID<<"ERROR3: source->getGlobalElement(exportLIDs[l]) is invalid GID="<<exportGIDs[i]<<" LID= "<<exportLIDs[i]<<std::endl;
      is_valid=false;
    }
  }
  
 // Zeroth order test: Remote *** GID *** and Export **GID**'s should be disjoint.  
  for( auto &&rgid : remoteGIDs) {
    if(std::find(exportGIDs.begin(),exportGIDs.end(),rgid) != exportGIDs.end()) {
        is_valid = false;
        os<<MyPID<<"ERROR0: Overlap between remoteGIDs and exportGIDs "<<rgid<<std::endl;
      }
  }
 
  int TempPID , OwningPID;
  for(GlobalOrdinal i=minSourceGID; i<maxSourceGID; i++) {
    // Get the (source) owner.
    // Abuse reductions to make up for the fact we don't know the owner is.
    // NOTE: If nobody owns this guy, it we'll get -1.
    LocalOrdinal slid = source->getLocalElement(i);    
    if(slid == LO_INVALID) TempPID = -1;
    else TempPID = MyPID;
    Teuchos::reduceAll<int, int> (*comm, Teuchos::REDUCE_MAX,TempPID, Teuchos::outArg(OwningPID));
 
    // Check to see if I have this guy in the target.  If so, make sure I am receiving him from the owner
    LocalOrdinal tlid = target->getLocalElement(i);    
 
    if(tlid != LO_INVALID) {
      // This guy is in my target map, now to check if I'm receiving him from the owner (which I now know)
      bool is_ok = false;
      
      // This guy is not in the SourceMap at all.  Weird, but acceptable.
      if(OwningPID == -1) continue;
 
      if (OwningPID == MyPID) {
  // I own this guy
  if((size_t) tlid < Importer.getNumSameIDs()) {
    // Check sames
    is_ok = true;
  }
  else {
    // Check permutes
    for (size_t j=0; j<(size_t)permuteTarget.size(); j++) {
      if(tlid == permuteTarget[j]) {
        is_ok=true; 
        break;
      }
    }
  }
      }
      else {
  // Check remotes
  bool already_hit = false;
  for(size_t j=0; j<(size_t)remoteGIDs.size(); j++) {
    if(i == remoteGIDs[j]) {
      // No double hits please
      if(already_hit) {
        is_ok=false; 
        break;
      }
      // GID's match:  Do procs?
      if(OwningPID == remotePIDs[j]) {
        is_ok = true;
        already_hit = true;
      }
    }
  }
      }
      if(!is_ok) {
        os<<MyPID<<"  ERROR1: GID "<<i<<" should be remoted from PID "<<OwningPID<<" but isn't."<<std::endl;
  is_valid=false;
      }
    }
 
  }//end for loop
 
  /***********************************************/
  /*              Check send side                */
  /***********************************************/
  Teuchos::Array<int> local_proc_mask(NumProcs,0), global_proc_mask(NumProcs,0);
 
 
  for(GlobalOrdinal i=minTargetGID; i<maxTargetGID; i++) {
 
    // If I have the target guy, set the proc mask
    LocalOrdinal tlid = target->getLocalElement(i);    
    LocalOrdinal slid = source->getLocalElement(i);    
 
    if(tlid==LO_INVALID) local_proc_mask[MyPID] = 0;
    else local_proc_mask[MyPID] = 1;
 
    Teuchos::reduceAll<int,int>(*comm,Teuchos::REDUCE_MAX,NumProcs, &local_proc_mask[0],&global_proc_mask[0]);
 
 
    if(slid !=LO_INVALID) {
      // If I own this unknown on the src I should check to make sure I'm exporting it to each guy in the global_proc_mask who wants it
      for(int j=0; j<NumProcs; j++) {
  if(j==MyPID) continue; // skip self unknowns
  if(global_proc_mask[j]==1) {
    bool is_ok = false;
    // This guy needs the unknown
    bool already_hit = false;
    for(size_t k=0; k<(size_t)exportPIDs.size(); k++) {
      if (exportPIDs[k] == j && source->getGlobalElement(exportLIDs[k]) == i) {
      // No double hits please
      if(already_hit) {
        is_ok=false; 
        break;
      }
      else {
        is_ok=true;
        already_hit=true;
      }
      }
    }
    if(!is_ok) {
      os<<MyPID<<" ERROR2: GID "<<i<<" should be sent to PID "<<j<<" but isn't"<<std::endl;
      is_valid=false;
    }
  }
      }
    }
  }
  
  // cbl check that for each of my remote GIDs I receive a corresponding export id. 
 
  Teuchos::Array<int> proc_num_exports_recv(NumProcs,0);
 
  Teuchos::Array<int> remoteGIDcount(remoteGIDs.size(),0);
 
  int allexpsiz=0;
  Teuchos::reduceAll<int,int>(*comm,Teuchos::REDUCE_MAX,exportGIDs.size(),  Teuchos::outArg(allexpsiz));
  
  for(int i=0;i<allexpsiz;++i) {
    Teuchos::Array<GlobalOrdinal> myexpgid(NumProcs,-2), yourexpgid(NumProcs,-2);
    Teuchos::Array<int> myexppid(NumProcs,-2), yourexppid(NumProcs,-2);
    if(i<exportGIDs.size()) {
      myexpgid[MyPID] = exportGIDs[i];
      myexppid[MyPID] = exportPIDs[i];
    }
    Teuchos::reduceAll<int,GlobalOrdinal>(*comm,Teuchos::REDUCE_MAX,NumProcs, &myexpgid[0],&yourexpgid[0]);
    Teuchos::reduceAll<int,int>(*comm,Teuchos::REDUCE_MAX,NumProcs, &myexppid[0],&yourexppid[0]);
    for(int p=0;p<NumProcs;++p) { // check one to one and onto
      GlobalOrdinal cgid = yourexpgid[p];
      // ignore -2's. 
      if(cgid == -2) continue;
      if(cgid < 0) {
        os<<MyPID<<" ERROR4: received exportGID is invalid "<<cgid<<std::endl;
        is_valid=false;
      }
      bool foundit=false;
      for(size_t k=0;k<(size_t)remoteGIDs.size();++k) {
        if(cgid == remoteGIDs[k] && yourexppid[p] == MyPID ) {
          if(p != remotePIDs[k]) {
            os<<MyPID<<" ERROR5: receive export from wrong pid: got "<<p<<" expected: "<<remotePIDs[k]<<std::endl;
            is_valid = false;
          }
          remoteGIDcount[k]++;
          if(foundit) {
            os<<MyPID<<" ERROR6: found multiple GIDs from correct pid: GID  "<<remoteGIDs[k]<<std::endl;
            is_valid = false;
          }
          foundit = true;
        }
      }
      if(!foundit &&  yourexppid[p] == MyPID ) {
        os<<MyPID<<" ERROR7: receive gid  "<<cgid<<" that is not in my remote gid list, from pid  "<<p<<std::endl;
        is_valid = false;
      }
 
    }
  }
  // now check that remoteGIDcount is only 1's.
  for(size_t i = 0; i< (size_t) remoteGIDcount.size(); ++i) {
    int rc = remoteGIDcount[i];
    if(rc == 1) continue;
    os<<MyPID<<" ERROR8: my remote at "<<i<<" gid "<<remoteGIDs[i]<<" has count "<<rc<<std::endl;
    is_valid = false;
  }
 
 
  // Do a reduction on the final bool status
  Teuchos::reduceAll<int,int> (*comm, Teuchos::REDUCE_MIN,(int)is_valid, Teuchos::outArg(global_is_valid));
 
  if(!global_is_valid) {
    std::cerr<<os.str()<<std::flush;
    Importer.print(std::cout);
  }
 
  return global_is_valid>0;
}
 
 
/* Check to see that the import object's target map respects AztecOO/ML ordering */
template <typename LocalOrdinal, typename GlobalOrdinal, typename Node>
bool
checkAztecOOMLOrdering (const Tpetra::Import<LocalOrdinal,GlobalOrdinal,Node>& Importer) 
{
 
  Teuchos::RCP<const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > source = Importer.getSourceMap();
  Teuchos::RCP<const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > target = Importer.getTargetMap();
  
  // Get the (pid, gid) pairs (with locals flagged as -1)
  Teuchos::Array<std::pair<int, GlobalOrdinal> > gpids;
  getPidGidPairs (Importer, gpids, true);
 
  bool is_ok=true;
  bool is_local = true;
  int last_PID = -1;
  GlobalOrdinal INVALID = Teuchos::OrdinalTraits<GlobalOrdinal>::invalid();
  GlobalOrdinal last_GID = Teuchos::OrdinalTraits<GlobalOrdinal>::invalid();
  
  for(size_t i=0; i<(size_t) gpids.size(); i++) {
    int pid = gpids[i].first;
    GlobalOrdinal gid = gpids[i].second;
     
    if(is_local == false && pid == -1) {
      // Out-of-order local
      is_ok = false;
      break;
    }
    else if(pid == -1) {
      // Local, matching PID
      if(source->getGlobalElement(i) != target->getGlobalElement(i)) {
        // GIDs don't match, though the PIDs do
        is_ok = false;
        break;
      }
    }
    else {
      // Off-rank section
      is_local = false;
      if(last_PID == -1) {
        last_PID = pid;
        last_GID = gid;
      }
      else if(pid < last_PID) {
        // PIDs out of order
        is_ok=false;
        break;
      }
      else if(pid == last_PID) {
        if(gid < last_GID) {
          // Same rank, but gids out of order
                is_ok=false;
                break;
        }
      }
      else {
        // New rank
        last_PID = pid;
        last_GID  = gid;
      }
    }
  }
 
  return is_ok;
}
 
 
/* Check to see that the import object's target map respects AztecOO/ML ordering */
template <typename LocalOrdinal, typename GlobalOrdinal, typename Node>
bool 
checkBlockConsistency(const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node>& map, size_t block_size) {
  bool consistent_block = true;
 
  for (size_t lid = 0; lid < map.getLocalNumElements (); lid += block_size) {
    auto lbid = floor(double(lid)/block_size);
    auto gbid = floor(double(map.getGlobalElement(lid))/block_size);
 
    for (size_t lid2 = 1; lid2 < block_size; ++lid2) {
      auto lbid_n = floor(double(lid+lid2)/block_size);
      auto gbid_n = floor(double(map.getGlobalElement(lid+lid2))/block_size);
      if (consistent_block)
        consistent_block = (lbid == lbid_n);
      if (consistent_block)
        consistent_block = (gbid == gbid_n);
    }
  }
 
  return consistent_block;
}
 
 
} // namespace Import_Util
} // namespace Tpetra
 
#endif // TPETRA_IMPORT_UTIL_HPP