MueLu_AlgebraicPermutationStrategy_def.hpp

Go to the documentation of this file.

// @HEADER
// *****************************************************************************
//        MueLu: A package for multigrid based preconditioning
//
// Copyright 2012 NTESS and the MueLu contributors.
// SPDX-License-Identifier: BSD-3-Clause
// *****************************************************************************
// @HEADER
 
/*
 * MueLu_AlgebraicPermutationStrategy_def.hpp
 *
 *  Created on: Feb 20, 2013
 *      Author: tobias
 */
 
#ifndef MUELU_ALGEBRAICPERMUTATIONSTRATEGY_DEF_HPP_
#define MUELU_ALGEBRAICPERMUTATIONSTRATEGY_DEF_HPP_
 
#include <queue>
 
#include <Teuchos_ScalarTraits.hpp>
 
#include <Xpetra_MultiVector.hpp>
#include <Xpetra_Matrix.hpp>
#include <Xpetra_CrsGraph.hpp>
#include <Xpetra_Vector.hpp>
#include <Xpetra_MultiVectorFactory.hpp>
#include <Xpetra_VectorFactory.hpp>
#include <Xpetra_CrsMatrixWrap.hpp>
#include <Xpetra_Export.hpp>
#include <Xpetra_ExportFactory.hpp>
#include <Xpetra_Import.hpp>
#include <Xpetra_ImportFactory.hpp>
#include <Xpetra_MatrixMatrix.hpp>
 
#include "MueLu_Utilities.hpp"
#include "MueLu_AlgebraicPermutationStrategy_decl.hpp"
 
namespace MueLu {
 
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void AlgebraicPermutationStrategy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
    BuildPermutation(const Teuchos::RCP<Matrix>& A, const Teuchos::RCP<const Map>& permRowMap,
                     Level& currentLevel, const FactoryBase* genFactory) const {
  const Scalar SC_ZERO = Teuchos::ScalarTraits<Scalar>::zero();
  const Scalar SC_ONE  = Teuchos::ScalarTraits<Scalar>::one();
 
  typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType MT;
  const MT MT_ZERO = Teuchos::ScalarTraits<MT>::zero();
  const MT MT_ONE  = Teuchos::ScalarTraits<MT>::one();
 
  const Teuchos::RCP<const Teuchos::Comm<int> > comm = A->getRowMap()->getComm();
  int numProcs                                       = comm->getSize();
  int myRank                                         = comm->getRank();
 
  size_t nDofsPerNode = 1;
  if (A->IsView("stridedMaps")) {
    Teuchos::RCP<const Map> permRowMapStrided = A->getRowMap("stridedMaps");
    nDofsPerNode                              = Teuchos::rcp_dynamic_cast<const StridedMap>(permRowMapStrided)->getFixedBlockSize();
  }
 
  std::vector<std::pair<GlobalOrdinal, GlobalOrdinal> > permutedDiagCandidates;
  std::vector<std::pair<GlobalOrdinal, GlobalOrdinal> > keepDiagonalEntries;
  std::vector<MT> Weights;
 
  // loop over all local rows in matrix A and keep diagonal entries if corresponding
  // matrix rows are not contained in permRowMap
  for (size_t row = 0; row < A->getRowMap()->getLocalNumElements(); row++) {
    GlobalOrdinal grow = A->getRowMap()->getGlobalElement(row);
 
    if (permRowMap->isNodeGlobalElement(grow) == true) continue;
 
    size_t nnz = A->getNumEntriesInLocalRow(row);
 
    // extract local row information from matrix
    Teuchos::ArrayView<const LocalOrdinal> indices;
    Teuchos::ArrayView<const Scalar> vals;
    A->getLocalRowView(row, indices, vals);
 
    TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::as<size_t>(indices.size()) != nnz, Exceptions::RuntimeError,
                               "MueLu::PermutationFactory::Build: number of nonzeros not equal to number of indices? Error.");
 
    // find column entry with max absolute value
    GlobalOrdinal gMaxValIdx = 0;
    MT norm1                 = MT_ZERO;
    MT maxVal                = MT_ZERO;
    for (size_t j = 0; j < Teuchos::as<size_t>(indices.size()); j++) {
      norm1 += Teuchos::ScalarTraits<Scalar>::magnitude(vals[j]);
      if (Teuchos::ScalarTraits<Scalar>::magnitude(vals[j]) > maxVal) {
        maxVal     = Teuchos::ScalarTraits<Scalar>::magnitude(vals[j]);
        gMaxValIdx = A->getColMap()->getGlobalElement(indices[j]);
      }
    }
 
    if (grow == gMaxValIdx)  // only keep row/col pair if it's diagonal dominant!!!
      keepDiagonalEntries.push_back(std::make_pair(grow, grow));
  }
 
  // handle rows that are marked to be relevant for permutations
  for (size_t row = 0; row < permRowMap->getLocalNumElements(); row++) {
    GlobalOrdinal grow = permRowMap->getGlobalElement(row);
    LocalOrdinal lArow = A->getRowMap()->getLocalElement(grow);
    size_t nnz         = A->getNumEntriesInLocalRow(lArow);
 
    // extract local row information from matrix
    Teuchos::ArrayView<const LocalOrdinal> indices;
    Teuchos::ArrayView<const Scalar> vals;
    A->getLocalRowView(lArow, indices, vals);
 
    TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::as<size_t>(indices.size()) != nnz, Exceptions::RuntimeError,
                               "MueLu::PermutationFactory::Build: number of nonzeros not equal to number of indices? Error.");
 
    // find column entry with max absolute value
    GlobalOrdinal gMaxValIdx = 0;
    MT norm1                 = MT_ZERO;
    MT maxVal                = MT_ZERO;
    for (size_t j = 0; j < Teuchos::as<size_t>(indices.size()); j++) {
      norm1 += Teuchos::ScalarTraits<Scalar>::magnitude(vals[j]);
      if (Teuchos::ScalarTraits<Scalar>::magnitude(vals[j]) > maxVal) {
        maxVal     = Teuchos::ScalarTraits<Scalar>::magnitude(vals[j]);
        gMaxValIdx = A->getColMap()->getGlobalElement(indices[j]);
      }
    }
 
    if (maxVal > Teuchos::ScalarTraits<MT>::zero()) {  // keep only max Entries \neq 0.0
      permutedDiagCandidates.push_back(std::make_pair(grow, gMaxValIdx));
      Weights.push_back(maxVal / (norm1 * Teuchos::as<MT>(nnz)));
    } else {
      std::cout << "ATTENTION: row " << grow << " has only zero entries -> singular matrix!" << std::endl;
    }
  }
 
  // sort Weights in descending order
  std::vector<int> permutation;
  sortingPermutation(Weights, permutation);
 
  // create new vector with exactly one possible entry for each column
 
  // each processor which requests the global column id gcid adds 1 to gColVec
  // gColVec will be summed up over all processors and communicated to gDomVec
  // which is based on the non-overlapping domain map of A.
 
  Teuchos::RCP<Vector> gColVec = VectorFactory::Build(A->getColMap());
  Teuchos::RCP<Vector> gDomVec = VectorFactory::Build(A->getDomainMap());
  gColVec->putScalar(SC_ZERO);
  gDomVec->putScalar(SC_ZERO);
 
  // put in all keep diagonal entries
  for (typename std::vector<std::pair<GlobalOrdinal, GlobalOrdinal> >::const_iterator p = keepDiagonalEntries.begin(); p != keepDiagonalEntries.end(); ++p) {
    gColVec->sumIntoGlobalValue((*p).second, 1.0);
  }
 
  Teuchos::RCP<Export> exporter = ExportFactory::Build(gColVec->getMap(), gDomVec->getMap());
  gDomVec->doExport(*gColVec, *exporter, Xpetra::ADD);  // communicate blocked gcolids to all procs
  gColVec->doImport(*gDomVec, *exporter, Xpetra::INSERT);
 
  std::vector<std::pair<GlobalOrdinal, GlobalOrdinal> > permutedDiagCandidatesFiltered;  // TODO reserve memory
  std::map<GlobalOrdinal, MT> gColId2Weight;
 
  {
    Teuchos::ArrayRCP<Scalar> ddata = gColVec->getDataNonConst(0);
    for (size_t i = 0; i < permutedDiagCandidates.size(); ++i) {
      // loop over all candidates
      std::pair<GlobalOrdinal, GlobalOrdinal> pp = permutedDiagCandidates[permutation[i]];
      GlobalOrdinal grow                         = pp.first;
      GlobalOrdinal gcol                         = pp.second;
 
      LocalOrdinal lcol = A->getColMap()->getLocalElement(gcol);
      if (Teuchos::ScalarTraits<Scalar>::real(ddata[lcol]) > MT_ZERO) {
        continue;  // skip lcol: column already handled by another row
      }
 
      // mark column as already taken
      ddata[lcol] += SC_ONE;
 
      permutedDiagCandidatesFiltered.push_back(std::make_pair(grow, gcol));
      gColId2Weight[gcol] = Weights[permutation[i]];
    }
  }
 
  // communicate how often each column index is requested by the different procs
  gDomVec->doExport(*gColVec, *exporter, Xpetra::ADD);
  gColVec->doImport(*gDomVec, *exporter, Xpetra::INSERT);  // probably not needed // TODO check me
 
  //*****************************************************************************************
  // first communicate ALL global ids of column indices which are requested by more
  // than one proc to all other procs
  // detect which global col indices are requested by more than one proc
  // and store them in the multipleColRequests vector
  std::vector<GlobalOrdinal> multipleColRequests;  // store all global column indices from current processor that are also
  // requested by another processor. This is possible, since they are stored
  // in gDomVec which is based on the nonoverlapping domain map. That is, each
  // global col id is handled by exactly one proc.
  std::queue<GlobalOrdinal> unusedColIdx;  // unused column indices on current processor
 
  for (size_t sz = 0; sz < gDomVec->getLocalLength(); ++sz) {
    Teuchos::ArrayRCP<const Scalar> arrDomVec = gDomVec->getData(0);
    //
    // FIXME (mfh 30 Oct 2015) I _think_ it's OK to check just the
    // real part, because this is a count.  (Shouldn't MueLu use
    // mag_type instead of Scalar here to save space?)
    //
    if (Teuchos::ScalarTraits<Scalar>::real(arrDomVec[sz]) > MT_ONE) {
      multipleColRequests.push_back(gDomVec->getMap()->getGlobalElement(sz));
    } else if (Teuchos::ScalarTraits<Scalar>::real(arrDomVec[sz]) == MT_ZERO) {
      unusedColIdx.push(gDomVec->getMap()->getGlobalElement(sz));
    }
  }
 
  // communicate the global number of column indices which are requested by more than one proc
  LocalOrdinal localMultColRequests  = Teuchos::as<LocalOrdinal>(multipleColRequests.size());
  LocalOrdinal globalMultColRequests = 0;
 
  // sum up all entries in multipleColRequests over all processors
  //
  // FIXME (lbv 11 Feb 2019) why cast localMultColRequests as LocalOrdinal
  // when it was properly declared as such just a line above?
  //
  MueLu_sumAll(gDomVec->getMap()->getComm(), (LocalOrdinal)localMultColRequests, globalMultColRequests);
 
  if (globalMultColRequests > 0) {
    // special handling: two processors request the same global column id.
    // decide which processor gets it
 
    // distribute number of multipleColRequests to all processors
    // each processor stores how many column ids for exchange are handled by the cur proc
    std::vector<GlobalOrdinal> numMyMultColRequests(numProcs, 0);
    std::vector<GlobalOrdinal> numGlobalMultColRequests(numProcs, 0);
    numMyMultColRequests[myRank] = localMultColRequests;
    Teuchos::reduceAll(*comm, Teuchos::REDUCE_MAX, numProcs,
                       numMyMultColRequests.data(),
                       numGlobalMultColRequests.data());
 
    // communicate multipleColRequests entries to all processors
    int nMyOffset = 0;
    for (int i = 0; i < myRank - 1; i++)
      nMyOffset += numGlobalMultColRequests[i];  // calculate offset to store the weights on the corresponding place in procOverlappingWeights
 
    const GlobalOrdinal zero = 0;
    std::vector<GlobalOrdinal> procMultRequestedColIds(globalMultColRequests, zero);
    std::vector<GlobalOrdinal> global_procMultRequestedColIds(globalMultColRequests, zero);
 
    // loop over all local column GIDs that are also requested by other procs
    for (size_t i = 0; i < multipleColRequests.size(); i++) {
      procMultRequestedColIds[nMyOffset + i] = multipleColRequests[i];  // all weights are > 0 ?
    }
 
    // template ordinal, package (double)
    Teuchos::reduceAll(*comm, Teuchos::REDUCE_MAX, Teuchos::as<int>(globalMultColRequests),
                       procMultRequestedColIds.data(),
                       global_procMultRequestedColIds.data());
 
    // loop over global_procOverlappingWeights and eliminate wrong entries...
    for (size_t k = 0; k < global_procMultRequestedColIds.size(); k++) {
      GlobalOrdinal globColId = global_procMultRequestedColIds[k];
 
      std::vector<MT> MyWeightForColId(numProcs, MT_ZERO);
      std::vector<MT> GlobalWeightForColId(numProcs, MT_ZERO);
 
      if (gColVec->getMap()->isNodeGlobalElement(globColId)) {
        MyWeightForColId[myRank] = gColId2Weight[globColId];
      } else {
        MyWeightForColId[myRank] = MT_ZERO;
      }
 
      Teuchos::reduceAll(*comm, Teuchos::REDUCE_MAX, numProcs,
                         MyWeightForColId.data(),
                         GlobalWeightForColId.data());
 
      if (gColVec->getMap()->isNodeGlobalElement(globColId)) {
        // note: 2 procs could have the same weight for a column index.
        // pick the first one.
        MT winnerValue     = MT_ZERO;
        int winnerProcRank = 0;
        for (int proc = 0; proc < numProcs; proc++) {
          if (GlobalWeightForColId[proc] > winnerValue) {
            winnerValue    = GlobalWeightForColId[proc];
            winnerProcRank = proc;
          }
        }
 
        // winnerProcRank is the winner for handling globColId.
        // winnerProcRank is unique (even if two procs have the same weight for a column index)
 
        if (myRank != winnerProcRank) {
          // remove corresponding entry from permutedDiagCandidatesFiltered
          typename std::vector<std::pair<GlobalOrdinal, GlobalOrdinal> >::iterator p = permutedDiagCandidatesFiltered.begin();
          while (p != permutedDiagCandidatesFiltered.end()) {
            if ((*p).second == globColId)
              p = permutedDiagCandidatesFiltered.erase(p);
            else
              p++;
          }
        }
 
      }  // end if isNodeGlobalElement
    }    // end loop over global_procOverlappingWeights and eliminate wrong entries...
  }      // end if globalMultColRequests > 0
 
  // put together all pairs:
  // size_t sizeRowColPairs = keepDiagonalEntries.size() + permutedDiagCandidatesFiltered.size();
  std::vector<std::pair<GlobalOrdinal, GlobalOrdinal> > RowColPairs;
  RowColPairs.insert(RowColPairs.end(), keepDiagonalEntries.begin(), keepDiagonalEntries.end());
  RowColPairs.insert(RowColPairs.end(), permutedDiagCandidatesFiltered.begin(), permutedDiagCandidatesFiltered.end());
 
#ifdef DEBUG_OUTPUT
  //&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
  // plausibility check
  gColVec->putScalar(SC_ZERO);
  gDomVec->putScalar(SC_ZERO);
  typename std::vector<std::pair<GlobalOrdinal, GlobalOrdinal> >::iterator pl = RowColPairs.begin();
  while (pl != RowColPairs.end()) {
    // GlobalOrdinal ik = (*pl).first;
    GlobalOrdinal jk = (*pl).second;
 
    gColVec->sumIntoGlobalValue(jk, 1.0);
    pl++;
  }
  gDomVec->doExport(*gColVec, *exporter, Xpetra::ADD);
  for (size_t sz = 0; sz < gDomVec->getLocalLength(); ++sz) {
    Teuchos::ArrayRCP<const Scalar> arrDomVec = gDomVec->getData(0);
    if (arrDomVec[sz] > 1.0) {
      GetOStream(Runtime0) << "RowColPairs has multiple column [" << sz << "]=" << arrDomVec[sz] << std::endl;
    } else if (arrDomVec[sz] == SC_ZERO) {
      GetOStream(Runtime0) << "RowColPairs has empty column [" << sz << "]=" << arrDomVec[sz] << std::endl;
    }
  }
  //&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
#endif
 
  // assumption: on each processor RowColPairs now contains
  // a valid set of (row,column) pairs, where the row entries
  // are a subset of the processor's rows and the column entries
  // are unique throughout all processors.
  // Note: the RowColPairs are only defined for a subset of all rows,
  // so there might be rows without an entry in RowColPairs.
  // It can be, that some rows seem to be missing in RowColPairs, since
  // the entry in that row with maximum absolute value has been reserved
  // by another row already (e.g. as already diagonal dominant row outside
  // of perRowMap).
  // In fact, the RowColPairs vector only defines the (row,column) pairs
  // that will be definitely moved to the diagonal after permutation.
 
#ifdef DEBUG_OUTPUT
  //  for (typename std::vector<std::pair<GlobalOrdinal, GlobalOrdinal> >::const_iterator p = RowColPairs.begin(); p != RowColPairs.end(); ++p) {
  //    std::cout << "proc: " << myRank << " r/c: " << (*p).first << "/" << (*p).second << std::endl;
  //  }
  //    for (typename std::vector<std::pair<GlobalOrdinal, GlobalOrdinal> >::const_iterator p = RowColPairs.begin(); p != RowColPairs.end(); ++p)
  //    {
  //      std::cout << (*p).first +1 << " " << (*p).second+1 << std::endl;
  //    }
  //    std::cout << "\n";
#endif
 
  // vectors to store permutation information
  Teuchos::RCP<Vector> Pperm  = VectorFactory::Build(A->getRowMap());
  Teuchos::RCP<Vector> Qperm  = VectorFactory::Build(A->getDomainMap());  // global variant (based on domain map)
  Teuchos::RCP<Vector> lQperm = VectorFactory::Build(A->getColMap());     // local variant (based on column map)
 
  Pperm->putScalar(SC_ZERO);
  Qperm->putScalar(SC_ZERO);
  lQperm->putScalar(SC_ZERO);
 
  // setup exporter for Qperm
  Teuchos::RCP<Export> QpermExporter = ExportFactory::Build(lQperm->getMap(), Qperm->getMap());
 
  Teuchos::RCP<Vector> RowIdStatus  = VectorFactory::Build(A->getRowMap());
  Teuchos::RCP<Vector> ColIdStatus  = VectorFactory::Build(A->getDomainMap());  // global variant (based on domain map)
  Teuchos::RCP<Vector> lColIdStatus = VectorFactory::Build(A->getColMap());     // local variant (based on column map)
  Teuchos::RCP<Vector> ColIdUsed    = VectorFactory::Build(A->getDomainMap());  // mark column ids to be already in use
  RowIdStatus->putScalar(SC_ZERO);
  ColIdStatus->putScalar(SC_ZERO);
  lColIdStatus->putScalar(SC_ZERO);
  ColIdUsed->putScalar(SC_ZERO);  // no column ids are used
 
  // count wide-range permutations
  // a wide-range permutation is defined as a permutation of rows/columns which do not
  // belong to the same node
  LocalOrdinal lWideRangeRowPermutations  = 0;
  GlobalOrdinal gWideRangeRowPermutations = 0;
  LocalOrdinal lWideRangeColPermutations  = 0;
  GlobalOrdinal gWideRangeColPermutations = 0;
 
  // run 1: mark all "identity" permutations
  {
    Teuchos::ArrayRCP<Scalar> RowIdStatusArray  = RowIdStatus->getDataNonConst(0);
    Teuchos::ArrayRCP<Scalar> lColIdStatusArray = lColIdStatus->getDataNonConst(0);
 
    typename std::vector<std::pair<GlobalOrdinal, GlobalOrdinal> >::iterator p = RowColPairs.begin();
    while (p != RowColPairs.end()) {
      GlobalOrdinal ik = (*p).first;
      GlobalOrdinal jk = (*p).second;
 
      LocalOrdinal lik = A->getRowMap()->getLocalElement(ik);
      LocalOrdinal ljk = A->getColMap()->getLocalElement(jk);
 
      if (RowIdStatusArray[lik] == SC_ZERO) {
        RowIdStatusArray[lik]  = SC_ONE;  // use this row id
        lColIdStatusArray[ljk] = SC_ONE;  // use this column id
        Pperm->replaceLocalValue(lik, ik);
        lQperm->replaceLocalValue(ljk, ik);         // use column map
        ColIdUsed->replaceGlobalValue(ik, SC_ONE);  // ik is now used
        p = RowColPairs.erase(p);
 
        // detect wide range permutations
        if (floor(ik / nDofsPerNode) != floor(jk / nDofsPerNode)) {
          lWideRangeColPermutations++;
        }
      } else
        p++;
    }
 
    // communicate column map -> domain map
    Qperm->doExport(*lQperm, *QpermExporter, Xpetra::ABSMAX);
    ColIdStatus->doExport(*lColIdStatus, *QpermExporter, Xpetra::ABSMAX);
 
    // plausibility check
    if (RowColPairs.size() > 0) GetOStream(Warnings0) << "MueLu::PermutationFactory: There are Row/Col pairs left!!!" << std::endl;  // TODO fix me
 
    // close Pperm
 
    // count, how many row permutations are missing on current proc
    size_t cntFreeRowIdx = 0;
    std::queue<GlobalOrdinal> qFreeGRowIdx;  // store global row ids of "free" rows
    for (size_t lik = 0; lik < RowIdStatus->getLocalLength(); ++lik) {
      if (RowIdStatusArray[lik] == SC_ZERO) {
        cntFreeRowIdx++;
        qFreeGRowIdx.push(RowIdStatus->getMap()->getGlobalElement(lik));
      }
    }
 
    // fix Pperm
    for (size_t lik = 0; lik < RowIdStatus->getLocalLength(); ++lik) {
      if (RowIdStatusArray[lik] == SC_ZERO) {
        RowIdStatusArray[lik] = SC_ONE;  // use this row id
        Pperm->replaceLocalValue(lik, qFreeGRowIdx.front());
        // detect wide range permutations
        if (floor(qFreeGRowIdx.front() / nDofsPerNode) != floor(RowIdStatus->getMap()->getGlobalElement(lik) / nDofsPerNode)) {
          lWideRangeRowPermutations++;
        }
        qFreeGRowIdx.pop();
      }
    }
  }
 
  size_t cntUnusedColIdx = 0;
  std::queue<GlobalOrdinal> qUnusedGColIdx;  // store global column ids of "free" available columns
  size_t cntFreeColIdx = 0;
  std::queue<GlobalOrdinal> qFreeGColIdx;  // store global column ids of "free" available columns
  {
    Teuchos::ArrayRCP<Scalar> ColIdStatusArray = ColIdStatus->getDataNonConst(0);
    Teuchos::ArrayRCP<Scalar> ColIdUsedArray   = ColIdUsed->getDataNonConst(0);  // not sure about this
 
    // close Qperm (free permutation entries in Qperm)
    for (size_t ljk = 0; ljk < ColIdStatus->getLocalLength(); ++ljk) {
      if (ColIdStatusArray[ljk] == SC_ZERO) {
        cntFreeColIdx++;
        qFreeGColIdx.push(ColIdStatus->getMap()->getGlobalElement(ljk));
      }
    }
 
    for (size_t ljk = 0; ljk < ColIdUsed->getLocalLength(); ++ljk) {
      if (ColIdUsedArray[ljk] == SC_ZERO) {
        cntUnusedColIdx++;
        qUnusedGColIdx.push(ColIdUsed->getMap()->getGlobalElement(ljk));
      }
    }
 
    // fix Qperm with local entries
    for (size_t ljk = 0; ljk < ColIdStatus->getLocalLength(); ++ljk) {
      // stop if no (local) unused column idx are left
      if (cntUnusedColIdx == 0) break;
 
      if (ColIdStatusArray[ljk] == SC_ZERO) {
        ColIdStatusArray[ljk] = SC_ONE;                                 // use this row id
        Qperm->replaceLocalValue(ljk, qUnusedGColIdx.front());          // loop over ColIdStatus (lives on domain map)
        ColIdUsed->replaceGlobalValue(qUnusedGColIdx.front(), SC_ONE);  // ljk is now used, too
        // detect wide range permutations
        if (floor(qUnusedGColIdx.front() / nDofsPerNode) != floor(ColIdStatus->getMap()->getGlobalElement(ljk) / nDofsPerNode)) {
          lWideRangeColPermutations++;
        }
        qUnusedGColIdx.pop();
        cntUnusedColIdx--;
        cntFreeColIdx--;
      }
    }
 
    // Qperm->doExport(*lQperm,*QpermExporter,Xpetra::ABSMAX); // no export necessary, since changes only locally
    // ColIdStatus->doExport(*lColIdStatus,*QpermExporter,Xpetra::ABSMAX);
 
    // count, how many unused column idx are needed on current processor
    // to complete Qperm
    cntFreeColIdx = 0;
    for (size_t ljk = 0; ljk < ColIdStatus->getLocalLength(); ++ljk) {  // TODO avoid this loop
      if (ColIdStatusArray[ljk] == SC_ZERO) {
        cntFreeColIdx++;
      }
    }
  }
 
  GlobalOrdinal global_cntFreeColIdx = 0;
  LocalOrdinal local_cntFreeColIdx   = cntFreeColIdx;
  MueLu_sumAll(comm, Teuchos::as<GlobalOrdinal>(local_cntFreeColIdx), global_cntFreeColIdx);
#ifdef DEBUG_OUTPUT
  std::cout << "global # of empty column idx entries in Qperm: " << global_cntFreeColIdx << std::endl;
#endif
 
  // avoid global communication if possible
  if (global_cntFreeColIdx > 0) {
    // 1) count how many unused column ids are left
    GlobalOrdinal global_cntUnusedColIdx = 0;
    LocalOrdinal local_cntUnusedColIdx   = cntUnusedColIdx;
    MueLu_sumAll(comm, Teuchos::as<GlobalOrdinal>(local_cntUnusedColIdx), global_cntUnusedColIdx);
#ifdef DEBUG_OUTPUT
    std::cout << "global # of unused column idx: " << global_cntUnusedColIdx << std::endl;
#endif
 
    // 2) communicate how many unused column ids are available on procs
    std::vector<LocalOrdinal> local_UnusedColIdxOnProc(numProcs);
    std::vector<LocalOrdinal> global_UnusedColIdxOnProc(numProcs);
    local_UnusedColIdxOnProc[myRank] = local_cntUnusedColIdx;
    Teuchos::reduceAll(*comm, Teuchos::REDUCE_MAX, numProcs,
                       local_UnusedColIdxOnProc.data(),
                       global_UnusedColIdxOnProc.data());
 
#ifdef DEBUG_OUTPUT
    std::cout << "PROC " << myRank << " global num unused indices per proc: ";
    for (size_t ljk = 0; ljk < global_UnusedColIdxOnProc.size(); ++ljk) {
      std::cout << " " << global_UnusedColIdxOnProc[ljk];
    }
    std::cout << std::endl;
#endif
 
    // 3) build array of length global_cntUnusedColIdx to globally replicate unused column idx
    std::vector<GlobalOrdinal> local_UnusedColIdxVector(Teuchos::as<size_t>(global_cntUnusedColIdx));
    std::vector<GlobalOrdinal> global_UnusedColIdxVector(Teuchos::as<size_t>(global_cntUnusedColIdx));
    GlobalOrdinal global_cntUnusedColIdxStartIter = 0;
    for (int proc = 0; proc < myRank; proc++) {
      global_cntUnusedColIdxStartIter += global_UnusedColIdxOnProc[proc];
    }
    for (GlobalOrdinal k = global_cntUnusedColIdxStartIter; k < global_cntUnusedColIdxStartIter + local_cntUnusedColIdx; k++) {
      local_UnusedColIdxVector[k] = qUnusedGColIdx.front();
      qUnusedGColIdx.pop();
    }
    Teuchos::reduceAll(*comm, Teuchos::REDUCE_MAX, Teuchos::as<int>(global_cntUnusedColIdx),
                       local_UnusedColIdxVector.data(),
                       global_UnusedColIdxVector.data());
#ifdef DEBUG_OUTPUT
    std::cout << "PROC " << myRank << " global UnusedGColIdx: ";
    for (size_t ljk = 0; ljk < global_UnusedColIdxVector.size(); ++ljk) {
      std::cout << " " << global_UnusedColIdxVector[ljk];
    }
    std::cout << std::endl;
#endif
 
    // 4) communicate, how many column idx are needed on each processor
    //    to complete Qperm
    std::vector<LocalOrdinal> local_EmptyColIdxOnProc(numProcs);
    std::vector<LocalOrdinal> global_EmptyColIdxOnProc(numProcs);
    local_EmptyColIdxOnProc[myRank] = local_cntFreeColIdx;
    Teuchos::reduceAll(*comm, Teuchos::REDUCE_MAX, numProcs,
                       local_EmptyColIdxOnProc.data(),
                       global_EmptyColIdxOnProc.data());
 
#ifdef DEBUG_OUTPUT
    std::cout << "PROC " << myRank << " global num of needed column indices: ";
    for (size_t ljk = 0; ljk < global_EmptyColIdxOnProc.size(); ++ljk) {
      std::cout << " " << global_EmptyColIdxOnProc[ljk];
    }
    std::cout << std::endl;
#endif
 
    // 5) determine first index in global_UnusedColIdxVector for unused column indices,
    //    that are marked to be used by this processor
    GlobalOrdinal global_UnusedColStartIdx = 0;
    for (int proc = 0; proc < myRank; proc++) {
      global_UnusedColStartIdx += global_EmptyColIdxOnProc[proc];
    }
 
#ifdef DEBUG_OUTPUT
    GetOStream(Statistics0) << "PROC " << myRank << " is allowd to use the following column gids: ";
    for (GlobalOrdinal k = global_UnusedColStartIdx; k < global_UnusedColStartIdx + Teuchos::as<GlobalOrdinal>(cntFreeColIdx); k++) {
      GetOStream(Statistics0) << global_UnusedColIdxVector[k] << " ";
    }
    GetOStream(Statistics0) << std::endl;
#endif
 
    // 6.) fix Qperm with global entries
    {
      Teuchos::ArrayRCP<Scalar> ColIdStatusArray = ColIdStatus->getDataNonConst(0);
      GlobalOrdinal array_iter                   = 0;
      for (size_t ljk = 0; ljk < ColIdStatus->getLocalLength(); ++ljk) {
        if (ColIdStatusArray[ljk] == SC_ZERO) {
          ColIdStatusArray[ljk] = SC_ONE;  // use this row id
          Qperm->replaceLocalValue(ljk, global_UnusedColIdxVector[global_UnusedColStartIdx + array_iter]);
          ColIdUsed->replaceGlobalValue(global_UnusedColIdxVector[global_UnusedColStartIdx + array_iter], SC_ONE);
          // detect wide range permutations
          if (floor(global_UnusedColIdxVector[global_UnusedColStartIdx + array_iter] / nDofsPerNode) != floor(ColIdStatus->getMap()->getGlobalElement(ljk) / nDofsPerNode)) {
            lWideRangeColPermutations++;
          }
          array_iter++;
          // cntUnusedColIdx--; // check me
        }
      }
    }
  }  // end if global_cntFreeColIdx > 0
 
  // create new empty Matrix
  Teuchos::RCP<CrsMatrixWrap> permPTmatrix = Teuchos::rcp(new CrsMatrixWrap(A->getRowMap(), 1));
  Teuchos::RCP<CrsMatrixWrap> permQTmatrix = Teuchos::rcp(new CrsMatrixWrap(A->getRowMap(), 1));
 
  {
    Teuchos::ArrayRCP<Scalar> PpermData = Pperm->getDataNonConst(0);
    Teuchos::ArrayRCP<Scalar> QpermData = Qperm->getDataNonConst(0);
 
    for (size_t row = 0; row < A->getLocalNumRows(); row++) {
      // FIXME (mfh 30 Oct 2015): Teuchos::as doesn't know how to
      // convert from complex Scalar to GO, so we have to take the real
      // part first.  I think that's the right thing to do in this case.
      Teuchos::ArrayRCP<GlobalOrdinal> indoutP(1, Teuchos::as<GO>(Teuchos::ScalarTraits<Scalar>::real(PpermData[row])));  // column idx for Perm^T
      Teuchos::ArrayRCP<GlobalOrdinal> indoutQ(1, Teuchos::as<GO>(Teuchos::ScalarTraits<Scalar>::real(QpermData[row])));  // column idx for Qperm
      Teuchos::ArrayRCP<Scalar> valout(1, SC_ONE);
      permPTmatrix->insertGlobalValues(A->getRowMap()->getGlobalElement(row), indoutP.view(0, indoutP.size()), valout.view(0, valout.size()));
      permQTmatrix->insertGlobalValues(A->getRowMap()->getGlobalElement(row), indoutQ.view(0, indoutQ.size()), valout.view(0, valout.size()));
    }
  }
 
  permPTmatrix->fillComplete();
  permQTmatrix->fillComplete();
 
  Teuchos::RCP<Matrix> permPmatrix = Utilities::Transpose(*permPTmatrix, true);
 
  for (size_t row = 0; row < permPTmatrix->getLocalNumRows(); row++) {
    if (permPTmatrix->getNumEntriesInLocalRow(row) != 1)
      GetOStream(Warnings0) << "#entries in row " << row << " of permPTmatrix is " << permPTmatrix->getNumEntriesInLocalRow(row) << std::endl;
    if (permPmatrix->getNumEntriesInLocalRow(row) != 1)
      GetOStream(Warnings0) << "#entries in row " << row << " of permPmatrix is " << permPmatrix->getNumEntriesInLocalRow(row) << std::endl;
    if (permQTmatrix->getNumEntriesInLocalRow(row) != 1)
      GetOStream(Warnings0) << "#entries in row " << row << " of permQmatrix is " << permQTmatrix->getNumEntriesInLocalRow(row) << std::endl;
  }
 
  // build permP * A * permQT
  Teuchos::RCP<Matrix> ApermQt      = Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Multiply(*A, false, *permQTmatrix, false, GetOStream(Statistics2));
  Teuchos::RCP<Matrix> permPApermQt = Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Multiply(*permPmatrix, false, *ApermQt, false, GetOStream(Statistics2));
 
  /*
    MueLu::Utils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Write("A.mat", *A);
    MueLu::Utils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Write("permP.mat", *permPmatrix);
    MueLu::Utils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Write("permQt.mat", *permQTmatrix);
    MueLu::Utils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Write("permPApermQt.mat", *permPApermQt);
  */
  // build scaling matrix
  Teuchos::RCP<Vector> diagVec              = VectorFactory::Build(permPApermQt->getRowMap(), true);
  Teuchos::RCP<Vector> invDiagVec           = VectorFactory::Build(permPApermQt->getRowMap(), true);
  Teuchos::RCP<CrsMatrixWrap> diagScalingOp = Teuchos::rcp(new CrsMatrixWrap(permPApermQt->getRowMap(), 1));
  {
    permPApermQt->getLocalDiagCopy(*diagVec);
    Teuchos::ArrayRCP<const Scalar> diagVecData = diagVec->getData(0);
    Teuchos::ArrayRCP<Scalar> invDiagVecData    = invDiagVec->getDataNonConst(0);
    for (size_t i = 0; i < diagVec->getMap()->getLocalNumElements(); ++i) {
      if (diagVecData[i] != SC_ZERO)
        invDiagVecData[i] = SC_ONE / diagVecData[i];
      else {
        invDiagVecData[i] = SC_ONE;
        GetOStream(Statistics0) << "MueLu::PermutationFactory: found zero on diagonal in row " << i << std::endl;
      }
    }
 
    for (size_t row = 0; row < A->getLocalNumRows(); row++) {
      Teuchos::ArrayRCP<GlobalOrdinal> indout(1, permPApermQt->getRowMap()->getGlobalElement(row));  // column idx for Perm^T
      Teuchos::ArrayRCP<Scalar> valout(1, invDiagVecData[row]);
      diagScalingOp->insertGlobalValues(A->getRowMap()->getGlobalElement(row), indout.view(0, indout.size()), valout.view(0, valout.size()));
    }
  }
  diagScalingOp->fillComplete();
 
  Teuchos::RCP<Matrix> scaledA = Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Multiply(*diagScalingOp, false, *permPApermQt, false, GetOStream(Statistics2));
  currentLevel.Set("A", Teuchos::rcp_dynamic_cast<Matrix>(scaledA), genFactory /*this*/);
 
  currentLevel.Set("permA", Teuchos::rcp_dynamic_cast<Matrix>(permPApermQt), genFactory /*this*/);  // TODO careful with this!!!
  currentLevel.Set("permP", Teuchos::rcp_dynamic_cast<Matrix>(permPmatrix), genFactory /*this*/);
  currentLevel.Set("permQT", Teuchos::rcp_dynamic_cast<Matrix>(permQTmatrix), genFactory /*this*/);
  currentLevel.Set("permScaling", Teuchos::rcp_dynamic_cast<Matrix>(diagScalingOp), genFactory /*this*/);
 
  // count zeros on diagonal in P -> number of row permutations
  Teuchos::RCP<Vector> diagPVec = VectorFactory::Build(permPmatrix->getRowMap(), true);
  permPmatrix->getLocalDiagCopy(*diagPVec);
  LocalOrdinal lNumRowPermutations  = 0;
  GlobalOrdinal gNumRowPermutations = 0;
  {
    Teuchos::ArrayRCP<const Scalar> diagPVecData = diagPVec->getData(0);
    for (size_t i = 0; i < diagPVec->getMap()->getLocalNumElements(); ++i) {
      if (diagPVecData[i] == SC_ZERO) {
        lNumRowPermutations++;
      }
    }
  }
 
  // sum up all entries in multipleColRequests over all processors
  MueLu_sumAll(diagPVec->getMap()->getComm(), Teuchos::as<GlobalOrdinal>(lNumRowPermutations), gNumRowPermutations);
 
  // count zeros on diagonal in Q^T -> number of column permutations
  Teuchos::RCP<Vector> diagQTVec = VectorFactory::Build(permQTmatrix->getRowMap(), true);
  permQTmatrix->getLocalDiagCopy(*diagQTVec);
  LocalOrdinal lNumColPermutations  = 0;
  GlobalOrdinal gNumColPermutations = 0;
  {
    Teuchos::ArrayRCP<const Scalar> diagQTVecData = diagQTVec->getData(0);
    for (size_t i = 0; i < diagQTVec->getMap()->getLocalNumElements(); ++i) {
      if (diagQTVecData[i] == SC_ZERO) {
        lNumColPermutations++;
      }
    }
  }
 
  // sum up all entries in multipleColRequests over all processors
  MueLu_sumAll(diagQTVec->getMap()->getComm(), Teuchos::as<GlobalOrdinal>(lNumColPermutations), gNumColPermutations);
 
  currentLevel.Set("#RowPermutations", gNumRowPermutations, genFactory /*this*/);
  currentLevel.Set("#ColPermutations", gNumColPermutations, genFactory /*this*/);
  currentLevel.Set("#WideRangeRowPermutations", gWideRangeRowPermutations, genFactory /*this*/);
  currentLevel.Set("#WideRangeColPermutations", gWideRangeColPermutations, genFactory /*this*/);
 
  GetOStream(Statistics0) << "#Row    permutations/max possible permutations: " << gNumRowPermutations << "/" << diagPVec->getMap()->getGlobalNumElements() << std::endl;
  GetOStream(Statistics0) << "#Column permutations/max possible permutations: " << gNumColPermutations << "/" << diagQTVec->getMap()->getGlobalNumElements() << std::endl;
  GetOStream(Runtime1) << "#wide range row permutations: " << gWideRangeRowPermutations << " #wide range column permutations: " << gWideRangeColPermutations << std::endl;
}
 
}  // namespace MueLu
 
#endif /* MUELU_ALGEBRAICPERMUTATIONSTRATEGY_DEF_HPP_ */