MueLu_MatrixAnalysisFactory_def.hpp

Go to the documentation of this file.

// @HEADER
//
// ***********************************************************************
//
//        MueLu: A package for multigrid based preconditioning
//                  Copyright 2012 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
//                    Jonathan Hu       (jhu@sandia.gov)
//                    Andrey Prokopenko (aprokop@sandia.gov)
//                    Ray Tuminaro      (rstumin@sandia.gov)
//
// ***********************************************************************
//
// @HEADER

#ifndef MUELU_MATRIXANALYSISFACTORY_DEF_HPP_
#define MUELU_MATRIXANALYSISFACTORY_DEF_HPP_

#include "MueLu_MatrixAnalysisFactory_decl.hpp"

#include <Xpetra_Map.hpp>
#include <Xpetra_StridedMap.hpp>    // for nDofsPerNode...
#include <Xpetra_Vector.hpp>
#include <Xpetra_VectorFactory.hpp>
#include <Xpetra_Matrix.hpp>

#include "MueLu_Level.hpp"
#include "MueLu_Utilities.hpp"
#include "MueLu_Monitor.hpp"

namespace MueLu {
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
MatrixAnalysisFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::MatrixAnalysisFactory()
  { }

template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
MatrixAnalysisFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::~MatrixAnalysisFactory() {}

template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
RCP<const ParameterList> MatrixAnalysisFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::GetValidParameterList() const {
  RCP<ParameterList> validParamList = rcp(new ParameterList());

  validParamList->set< RCP<const FactoryBase> >("A", Teuchos::null, "Generating factory of the matrix A to be permuted.");

  return validParamList;
}

template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void MatrixAnalysisFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::DeclareInput(Level &fineLevel, Level &/* coarseLevel */) const {
  Input(fineLevel, "A");
}

template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void MatrixAnalysisFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(Level &fineLevel, Level &/* coarseLevel */) const {
  FactoryMonitor m(*this, "MatrixAnalysis Factory ", fineLevel);

  Teuchos::RCP<Matrix> A = Get< Teuchos::RCP<Matrix> > (fineLevel, "A");
  Teuchos::RCP<const Teuchos::Comm<int> > comm = A->getRangeMap()->getComm();

  //const ParameterList & pL = GetParameterList();

  // General information
  {
    GetOStream(Runtime0) << "~~~~~~~~ GENERAL INFORMATION ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ (Level " << fineLevel.GetLevelID() << ")" << std::endl;
    GetOStream(Runtime0) << "A is a " << A->getRangeMap()->getGlobalNumElements() << " x " << A->getDomainMap()->getGlobalNumElements() << " matrix." << std::endl;

    if (A->IsView("stridedMaps") && Teuchos::rcp_dynamic_cast<const StridedMap>(A->getRowMap("stridedMaps")) != Teuchos::null) {
      Teuchos::RCP<const StridedMap> strRowMap = Teuchos::rcp_dynamic_cast<const StridedMap>(A->getRowMap("stridedMaps"));
      LocalOrdinal blockid       = strRowMap->getStridedBlockId();
      if (blockid > -1) {
        std::vector<size_t> stridingInfo = strRowMap->getStridingData();
        LO dofsPerNode = Teuchos::as<LocalOrdinal>(stridingInfo[blockid]);
        GetOStream(Runtime0) << "Strided row: " << dofsPerNode << " dofs per node. Strided block id = " << blockid << std::endl;
      } else {
        GetOStream(Runtime0) << "Row: " << strRowMap->getFixedBlockSize() << " dofs per node" << std::endl;
      }
    }

    if (A->IsView("stridedMaps") && Teuchos::rcp_dynamic_cast<const StridedMap>(A->getColMap("stridedMaps")) != Teuchos::null) {
      Teuchos::RCP<const StridedMap> strDomMap = Teuchos::rcp_dynamic_cast<const StridedMap>(A->getColMap("stridedMaps"));
      LocalOrdinal blockid       = strDomMap->getStridedBlockId();
      if (blockid > -1) {
        std::vector<size_t> stridingInfo = strDomMap->getStridingData();
        LO dofsPerNode = Teuchos::as<LocalOrdinal>(stridingInfo[blockid]);
        GetOStream(Runtime0) << "Strided column: " << dofsPerNode << " dofs per node. Strided block id = " << blockid << std::endl;
      } else {
        GetOStream(Runtime0) << "Column: " << strDomMap->getFixedBlockSize() << " dofs per node" << std::endl;
      }
    }


    GetOStream(Runtime0) << "A is distributed over " << comm->getSize() << " processors" << std::endl;

    // empty processors
    std::vector<LO> lelePerProc(comm->getSize(),0);
    std::vector<LO> gelePerProc(comm->getSize(),0);
    lelePerProc[comm->getRank()] = A->getNodeNumEntries ();
    Teuchos::reduceAll(*comm,Teuchos::REDUCE_MAX,comm->getSize(),&lelePerProc[0],&gelePerProc[0]);
    if(comm->getRank() == 0) {
      for(int i=0; i<comm->getSize(); i++) {
        if(gelePerProc[i] == 0) {
          GetOStream(Runtime0) << "Proc " << i << " is empty." << std::endl;
        }
      }
    }
  }

  // Matrix diagonal
  {
    GetOStream(Runtime0) << "~~~~~~~~ MATRIX DIAGONAL ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ (Level " << fineLevel.GetLevelID() << ")" << std::endl;
    Teuchos::RCP<Vector> diagAVec = VectorFactory::Build(A->getRowMap(),true);
    A->getLocalDiagCopy(*diagAVec);
    Teuchos::ArrayRCP< const Scalar > diagAVecData = diagAVec->getData(0);
    GlobalOrdinal lzerosOnDiagonal = 0;
    GlobalOrdinal gzerosOnDiagonal = 0;
    GlobalOrdinal lnearlyzerosOnDiagonal = 0;
    GlobalOrdinal gnearlyzerosOnDiagonal = 0;
    GlobalOrdinal lnansOnDiagonal = 0;
    GlobalOrdinal gnansOnDiagonal = 0;

    for(size_t i = 0; i<diagAVec->getMap()->getNodeNumElements(); ++i) {
      if(diagAVecData[i] == Teuchos::ScalarTraits<Scalar>::zero()) {
        lzerosOnDiagonal++;
      } else if(Teuchos::ScalarTraits<Scalar>::magnitude(diagAVecData[i]) < 1e-6) {
        lnearlyzerosOnDiagonal++;
      } else if(Teuchos::ScalarTraits<Scalar>::isnaninf(diagAVecData[i])) {
        lnansOnDiagonal++;
      }
    }

    // sum up all entries in multipleColRequests over all processors
    MueLu_sumAll(comm, lzerosOnDiagonal, gzerosOnDiagonal);
    MueLu_sumAll(comm, lnearlyzerosOnDiagonal, gnearlyzerosOnDiagonal);
    MueLu_sumAll(comm, lnansOnDiagonal, gnansOnDiagonal);

    if(gzerosOnDiagonal > 0) GetOStream(Runtime0) << "Found " << gzerosOnDiagonal << " zeros on diagonal of A" << std::endl;
    if(gnearlyzerosOnDiagonal > 0) GetOStream(Runtime0) << "Found " << gnearlyzerosOnDiagonal << " entries with absolute value < 1.0e-6 on diagonal of A" << std::endl;
    if(gnansOnDiagonal > 0) GetOStream(Runtime0) << "Found " << gnansOnDiagonal << " entries with NAN or INF on diagonal of A" << std::endl;
  }

  // Diagonal dominance?
  {
    GetOStream(Runtime0) << "~~~~~~~~ DIAGONAL DOMINANCE ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ (Level " << fineLevel.GetLevelID() << ")" << std::endl;
    // loop over all local rows in matrix A and keep diagonal entries if corresponding
    // matrix rows are not contained in permRowMap
    GlobalOrdinal lnumWeakDiagDomRows = 0;
    GlobalOrdinal gnumWeakDiagDomRows = 0;
    GlobalOrdinal lnumStrictDiagDomRows = 0;
    GlobalOrdinal gnumStrictDiagDomRows = 0;
    double worstRatio = 99999999.9;
    for (size_t row = 0; row < A->getRowMap()->getNodeNumElements(); row++) {
      GlobalOrdinal grow = A->getRowMap()->getGlobalElement(row);

      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::MatrixAnalysisFactory::Build: number of nonzeros not equal to number of indices? Error.");

      // find column entry with max absolute value
      double norm1 = 0.0;
      double normdiag = 0.0;
      for (size_t j = 0; j < Teuchos::as<size_t>(indices.size()); j++) {
        GO gcol = A->getColMap()->getGlobalElement(indices[j]);
        if (gcol==grow) normdiag = Teuchos::as<double>(Teuchos::ScalarTraits<Scalar>::magnitude(vals[j]));
        else norm1 +=  Teuchos::as<double>(Teuchos::ScalarTraits<Scalar>::magnitude(vals[j]));
      }

      if      (normdiag >= norm1) lnumWeakDiagDomRows++;
      else if (normdiag >  norm1) lnumStrictDiagDomRows++;

      if (norm1 != 0.0) {
        if (normdiag / norm1 < worstRatio) worstRatio = normdiag / norm1;
      }
    }

    MueLu_sumAll(comm, lnumWeakDiagDomRows, gnumWeakDiagDomRows);
    MueLu_sumAll(comm, lnumStrictDiagDomRows, gnumStrictDiagDomRows);

    GetOStream(Runtime0) << "A has " << gnumWeakDiagDomRows << "/" << A->getRangeMap()->getGlobalNumElements() << " weakly diagonal dominant rows. (" << Teuchos::as<Scalar>(gnumWeakDiagDomRows/A->getRangeMap()->getGlobalNumElements()*100) << "%)" << std::endl;
    GetOStream(Runtime0) << "A has " << gnumStrictDiagDomRows << "/" << A->getRangeMap()->getGlobalNumElements() << " strictly diagonal dominant rows. (" << Teuchos::as<Scalar>(gnumStrictDiagDomRows/A->getRangeMap()->getGlobalNumElements()*100) << "%)" << std::endl;

    double gworstRatio;
    Teuchos::reduceAll(*comm,Teuchos::REDUCE_MIN,comm->getSize(),&worstRatio,&gworstRatio);
    GetOStream(Runtime0) << "The minimum of the ratio of diagonal element/ sum of off-diagonal elements is " << gworstRatio << ". Values about 1.0 are ok." << std::endl;
  }


  SC one = Teuchos::ScalarTraits<Scalar>::one(), zero = Teuchos::ScalarTraits<Scalar>::zero();

  // multiply with one vector
  {
    GetOStream(Runtime0) << "~~~~~~~~ Av with one vector ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ (Level " << fineLevel.GetLevelID() << ")" << std::endl;
    Teuchos::RCP<Vector> ones = VectorFactory::Build(A->getDomainMap(), 1);
    ones->putScalar(one);

    Teuchos::RCP<Vector> res1 = VectorFactory::Build(A->getRangeMap(), false);
    A->apply(*ones, *res1, Teuchos::NO_TRANS, one, zero);
    checkVectorEntries(res1,std::string("after applying the one vector to A"));
  }

  {
    GetOStream(Runtime0) << "~~~~~~~~ Av with random vector ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ (Level " << fineLevel.GetLevelID() << ")" << std::endl;
    Teuchos::RCP<Vector> randvec = VectorFactory::Build(A->getDomainMap(), 1);
    randvec->randomize();

    Teuchos::RCP<Vector> resrand = VectorFactory::Build(A->getRangeMap(), false);
    A->apply(*randvec, *resrand, Teuchos::NO_TRANS, one, zero);
    checkVectorEntries(resrand,std::string("after applying random vector to A"));
  }

  // apply Jacobi sweep
  {
    GetOStream(Runtime0) << "~~~~~~~~ Damped Jacobi sweep (one vector) ~~~~~~~~~~~~~~~~~~~~~~~ (Level " << fineLevel.GetLevelID() << ")" << std::endl;
    Teuchos::RCP<Vector> ones = VectorFactory::Build(A->getDomainMap(), 1);
    ones->putScalar(one);

    Teuchos::RCP<Vector> res1 = VectorFactory::Build(A->getRangeMap(), false);
    A->apply(*ones, *res1, Teuchos::NO_TRANS, one, zero);
    checkVectorEntries(res1,std::string("after applying one vector to A"));

    Teuchos::RCP<Vector> invDiag = Utilities::GetMatrixDiagonalInverse(*A);
    checkVectorEntries(invDiag,std::string("in invDiag"));

    Teuchos::RCP<Vector> res2 = VectorFactory::Build(A->getRangeMap(), false);
    res2->elementWiseMultiply (0.8, *invDiag, *res1, 0.0);
    checkVectorEntries(res2,std::string("after scaling Av with invDiag (with v the one vector)"));
    res2->update(1.0, *ones, -1.0);

    checkVectorEntries(res2,std::string("after applying one damped Jacobi sweep (with one vector)"));
  }

  // apply Jacobi sweep
  {
    GetOStream(Runtime0) << "~~~~~~~~ Damped Jacobi sweep (random vector) ~~~~~~~~~~~~~~~~~~~~ (Level " << fineLevel.GetLevelID() << ")" << std::endl;
    Teuchos::RCP<Vector> ones = VectorFactory::Build(A->getDomainMap(), 1);
    ones->randomize();

    Teuchos::RCP<Vector> res1 = VectorFactory::Build(A->getRangeMap(), false);
    A->apply(*ones, *res1, Teuchos::NO_TRANS, one, zero);
    checkVectorEntries(res1,std::string("after applying a random vector to A"));

    Teuchos::RCP<Vector> invDiag = Utilities::GetMatrixDiagonalInverse(*A);
    checkVectorEntries(invDiag,std::string("in invDiag"));

    Teuchos::RCP<Vector> res2 = VectorFactory::Build(A->getRangeMap(), false);
    res2->elementWiseMultiply (0.8, *invDiag, *res1, 0.0);
    checkVectorEntries(res2,std::string("after scaling Av with invDiag (with v a random vector)"));
    res2->update(1.0, *ones, -1.0);

    checkVectorEntries(res2,std::string("after applying one damped Jacobi sweep (with v a random vector)"));
  }

  GetOStream(Runtime0) << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ (Level " << fineLevel.GetLevelID() << ")" << std::endl;
}

template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void MatrixAnalysisFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::checkVectorEntries(const Teuchos::RCP<Vector>& vec, std::string info) const {
  SC zero = Teuchos::ScalarTraits<Scalar>::zero();
  Teuchos::RCP<const Teuchos::Comm<int> > comm = vec->getMap()->getComm();
  Teuchos::ArrayRCP< const Scalar > vecData = vec->getData(0);
  GO lzeros = 0;
  GO gzeros = 0;
  GO lnans = 0;
  GO gnans = 0;

  for(size_t i = 0; i<vec->getMap()->getNodeNumElements(); ++i) {
    if(vecData[i] == zero) {
      lzeros++;
    } else if(Teuchos::ScalarTraits<Scalar>::isnaninf(vecData[i])) {
      lnans++;
    }
  }

  // sum up all entries in multipleColRequests over all processors
  MueLu_sumAll(comm, lzeros, gzeros);
  MueLu_sumAll(comm, lnans, gnans);

  if(gzeros > 0) GetOStream(Runtime0) << "Found " << gzeros << " zeros " << info << std::endl;
  if(gnans > 0) GetOStream(Runtime0) << "Found " << gnans << " entries " << info << std::endl;
}

} // namespace MueLu


#endif /* MUELU_MATRIXANALYSISFACTORY_DEF_HPP_ */