MueLu_RAPShiftFactory_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_RAPSHIFTFACTORY_DEF_HPP
#define MUELU_RAPSHIFTFACTORY_DEF_HPP

#include <sstream>

#include <Xpetra_Matrix.hpp>
#include <Xpetra_MatrixMatrix.hpp>
#include <Xpetra_MatrixUtils.hpp>
#include <Xpetra_Vector.hpp>
#include <Xpetra_VectorFactory.hpp>


#include "MueLu_RAPShiftFactory_decl.hpp"
#include "MueLu_MasterList.hpp"
#include "MueLu_Monitor.hpp"
#include "MueLu_PerfUtils.hpp"
#include "MueLu_Utilities.hpp"

namespace MueLu {

  /*********************************************************************************************************/
  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
  RAPShiftFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::RAPShiftFactory()
    : implicitTranspose_(false)  { }


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

#define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name))
    SET_VALID_ENTRY("transpose: use implicit");
    SET_VALID_ENTRY("rap: fix zero diagonals");
    SET_VALID_ENTRY("rap: shift");
    SET_VALID_ENTRY("rap: shift array");
    SET_VALID_ENTRY("rap: cfl array");
    SET_VALID_ENTRY("rap: shift diagonal M");
    SET_VALID_ENTRY("rap: shift low storage");
    SET_VALID_ENTRY("rap: relative diagonal floor");
#undef  SET_VALID_ENTRY

    validParamList->set< RCP<const FactoryBase> >("A",              Teuchos::null, "Generating factory of the matrix A used during the prolongator smoothing process");
    validParamList->set< RCP<const FactoryBase> >("M",              Teuchos::null, "Generating factory of the matrix M used during the non-Galerkin RAP");
    validParamList->set< RCP<const FactoryBase> >("Mdiag",          Teuchos::null, "Generating factory of the matrix Mdiag used during the non-Galerkin RAP");
    validParamList->set< RCP<const FactoryBase> >("K",              Teuchos::null, "Generating factory of the matrix K used during the non-Galerkin RAP");
    validParamList->set< RCP<const FactoryBase> >("P",              Teuchos::null, "Prolongator factory");
    validParamList->set< RCP<const FactoryBase> >("R",              Teuchos::null, "Restrictor factory");

    validParamList->set< bool >                  ("CheckMainDiagonal",  false, "Check main diagonal for zeros");
    validParamList->set< bool >                  ("RepairMainDiagonal", false, "Repair zeros on main diagonal");

    validParamList->set<RCP<const FactoryBase> > ("deltaT", Teuchos::null, "user deltaT");
    validParamList->set<RCP<const FactoryBase> > ("cfl", Teuchos::null, "user cfl");
    validParamList->set<RCP<const FactoryBase> > ("cfl-based shift array", Teuchos::null, "MueLu-generated shift array for CFL-based shifting");

    // Make sure we don't recursively validate options for the matrixmatrix kernels
    ParameterList norecurse;
    norecurse.disableRecursiveValidation();
    validParamList->set<ParameterList> ("matrixmatrix: kernel params", norecurse, "MatrixMatrix kernel parameters");

    return validParamList;
  }

  /*********************************************************************************************************/
  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
  void RAPShiftFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::DeclareInput(Level &fineLevel, Level &coarseLevel) const {
    const Teuchos::ParameterList& pL = GetParameterList();

    bool use_mdiag = false;
    if(pL.isParameter("rap: shift diagonal M")) 
      use_mdiag = pL.get<bool>("rap: shift diagonal M");

    // The low storage version requires mdiag
    bool use_low_storage = false;
    if(pL.isParameter("rap: shift low storage")) {
      use_low_storage = pL.get<bool>("rap: shift low storage");
      use_mdiag = use_low_storage ? true : use_mdiag;
    }

    if (implicitTranspose_ == false) {
      Input(coarseLevel, "R");
    }

    if(!use_low_storage) Input(fineLevel,   "K");
    else Input(fineLevel,   "A");
    Input(coarseLevel, "P");

    if(!use_mdiag) Input(fineLevel, "M");
    else Input(fineLevel, "Mdiag");

    // CFL array stuff
    if(pL.isParameter("rap: cfl array") && pL.get<Teuchos::Array<double> >("rap: cfl array").size() > 0) {
      if(fineLevel.GetLevelID() == 0) {
        if(fineLevel.IsAvailable("deltaT", NoFactory::get())) {
          fineLevel.DeclareInput("deltaT", NoFactory::get(), this);
        } else {
          TEUCHOS_TEST_FOR_EXCEPTION(fineLevel.IsAvailable("fine deltaT", NoFactory::get()),
                                     Exceptions::RuntimeError,
                                     "deltaT was not provided by the user on level0!");
        }
        
        if(fineLevel.IsAvailable("cfl", NoFactory::get())) {
          fineLevel.DeclareInput("cfl", NoFactory::get(), this);
        } else {
          TEUCHOS_TEST_FOR_EXCEPTION(fineLevel.IsAvailable("fine cfl", NoFactory::get()),
                                   Exceptions::RuntimeError,
                                     "cfl was not provided by the user on level0!");
        }      
      }
      else {
        Input(fineLevel,"cfl-based shift array");
      }
    }

    // call DeclareInput of all user-given transfer factories
    for(std::vector<RCP<const FactoryBase> >::const_iterator it = transferFacts_.begin(); it!=transferFacts_.end(); ++it) {
      (*it)->CallDeclareInput(coarseLevel);
    }
  }

  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
  void RAPShiftFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(Level &fineLevel, Level &coarseLevel) const { // FIXME make fineLevel const
    {
      FactoryMonitor m(*this, "Computing Ac", coarseLevel);      
      const Teuchos::ParameterList& pL = GetParameterList();

      bool M_is_diagonal = false; 
      if(pL.isParameter("rap: shift diagonal M"))
         M_is_diagonal = pL.get<bool>("rap: shift diagonal M");

      // The low storage version requires mdiag
      bool use_low_storage = false;
      if(pL.isParameter("rap: shift low storage")) {
        use_low_storage = pL.get<bool>("rap: shift low storage");
        M_is_diagonal = use_low_storage ? true : M_is_diagonal;
      }

      Teuchos::ArrayView<const double> doubleShifts;
      Teuchos::ArrayRCP<double> myshifts;
      if(pL.isParameter("rap: shift array") && pL.get<Teuchos::Array<double> >("rap: shift array").size() > 0 ) {
        // Do we have an array of shifts?  If so, we set doubleShifts_
        doubleShifts = pL.get<Teuchos::Array<double> >("rap: shift array")();
      }
      if(pL.isParameter("rap: cfl array") && pL.get<Teuchos::Array<double> >("rap: cfl array").size() > 0) {
        // Do we have an array of CFLs?  If so, we calculated the shifts from them.        
        Teuchos::ArrayView<const double> CFLs = pL.get<Teuchos::Array<double> >("rap: cfl array")();
        if(fineLevel.GetLevelID() == 0) {
          double dt  = Get<double>(fineLevel,"deltaT");
          double cfl = Get<double>(fineLevel,"cfl");
          double ts_at_cfl1 = dt / cfl;
          myshifts.resize(CFLs.size());
          Teuchos::Array<double> myCFLs(CFLs.size());
          myCFLs[0] = cfl;

          // Never make the CFL bigger
          for(int i=1; i<(int)CFLs.size(); i++)
            myCFLs[i] = (CFLs[i]> cfl) ? cfl : CFLs[i];
          
          {
            std::ostringstream ofs;
            ofs<<"RAPShiftFactory: CFL schedule = ";
            for(int i=0; i<(int)CFLs.size(); i++)
              ofs<<" "<<myCFLs[i];
            GetOStream(Statistics0) <<ofs.str() << std::endl;
          }
          GetOStream(Statistics0)<< "RAPShiftFactory: Timestep at CFL=1 is "<< ts_at_cfl1 << "    " <<std::endl;

          // The shift array needs to be 1/dt
          for(int i=0; i<(int)myshifts.size(); i++)
            myshifts[i] =  1.0 / (ts_at_cfl1*myCFLs[i]);
          doubleShifts = myshifts();
          
          {
            std::ostringstream ofs;
            ofs<<"RAPShiftFactory: shift schedule = ";
            for(int i=0; i<(int)doubleShifts.size(); i++)
              ofs<<" "<<doubleShifts[i];
            GetOStream(Statistics0) <<ofs.str() << std::endl;
          }
          Set(coarseLevel,"cfl-based shift array",myshifts);
        }
        else {
          myshifts = Get<Teuchos::ArrayRCP<double> > (fineLevel,"cfl-based shift array");
          doubleShifts = myshifts();
          Set(coarseLevel,"cfl-based shift array",myshifts);
          // NOTE: If we're not on level zero, then we should have a shift array
        }
      }
      
      // Inputs: K, M, P
      // Note: In the low-storage case we do not keep a separate "K", we just use A
      RCP<Matrix> K;
      RCP<Matrix> M;
      RCP<Vector> Mdiag;

      if(use_low_storage)  K = Get< RCP<Matrix> >(fineLevel, "A");
      else K = Get< RCP<Matrix> >(fineLevel, "K");
      if(!M_is_diagonal) M = Get< RCP<Matrix> >(fineLevel, "M");
      else Mdiag = Get< RCP<Vector> >(fineLevel, "Mdiag");

      RCP<Matrix> P = Get< RCP<Matrix> >(coarseLevel, "P");

      // Build Kc = RKP, Mc = RMP
      RCP<Matrix> KP, MP;

      // Reuse pattern if available (multiple solve)
      // FIXME: Old style reuse doesn't work any more
      //      if (IsAvailable(coarseLevel, "AP Pattern")) {
      //        KP = Get< RCP<Matrix> >(coarseLevel, "AP Pattern");
      //        MP = Get< RCP<Matrix> >(coarseLevel, "AP Pattern");
      //      }

      {
        SubFactoryMonitor subM(*this, "MxM: K x P", coarseLevel);
        KP = Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Multiply(*K, false, *P, false, KP, GetOStream(Statistics2));
        if(!M_is_diagonal) {
          MP = Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Multiply(*M, false, *P, false, MP, GetOStream(Statistics2));
        }
        else {
          MP = Xpetra::MatrixFactory2<Scalar, LocalOrdinal, GlobalOrdinal, Node>::BuildCopy(P);
          MP->leftScale(*Mdiag);
        }

        Set(coarseLevel, "AP Pattern", KP);
      }

      bool doOptimizedStorage = true;

      RCP<Matrix> Ac, Kc, Mc;

      // Reuse pattern if available (multiple solve)
      //     if (IsAvailable(coarseLevel, "RAP Pattern"))
      // Ac = Get< RCP<Matrix> >(coarseLevel, "RAP Pattern");

      bool doFillComplete=true;
      if (implicitTranspose_) {
        SubFactoryMonitor m2(*this, "MxM: P' x (KP) (implicit)", coarseLevel);
        Kc = Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Multiply(*P, true, *KP, false, Kc, GetOStream(Statistics2), doFillComplete, doOptimizedStorage);
        Mc = Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Multiply(*P, true, *MP, false, Mc, GetOStream(Statistics2), doFillComplete, doOptimizedStorage);
      }
      else {
        RCP<Matrix> R = Get< RCP<Matrix> >(coarseLevel, "R");
        SubFactoryMonitor m2(*this, "MxM: R x (KP) (explicit)", coarseLevel);
        Kc = Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Multiply(*R, false, *KP, false, Kc, GetOStream(Statistics2), doFillComplete, doOptimizedStorage);
        Mc = Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Multiply(*R, false, *MP, false, Mc, GetOStream(Statistics2), doFillComplete, doOptimizedStorage);
      }

      // Get the shift
      // FIXME - We should really get rid of the shifts array and drive this the same way everything else works
      // If we're using the recursive "low storage" version, we need to shift by ( \prod_{i=1}^k shift[i] - \prod_{i=1}^{k-1} shift[i]) to
      // get the recursive relationships correct
      int level      = coarseLevel.GetLevelID();
      Scalar shift   = Teuchos::ScalarTraits<Scalar>::zero();
      if(!use_low_storage) {
        // High Storage version
        if(level < (int)shifts_.size()) shift = shifts_[level];
        else shift  = Teuchos::as<Scalar>(pL.get<double>("rap: shift"));
      }
      else {
        // Low Storage Version
        if(level < (int)shifts_.size()) {
          if(level==1) shift = shifts_[level];
          else {
            Scalar prod1 = Teuchos::ScalarTraits<Scalar>::one();
            for(int i=1; i < level-1; i++) {
              prod1 *= shifts_[i];
            } 
            shift = (prod1 * shifts_[level] - prod1);
          }
        }
        else if(doubleShifts.size() != 0) {
          double d_shift = 0.0;
          if(level < doubleShifts.size())
            d_shift = doubleShifts[level] - doubleShifts[level-1];

          if(d_shift < 0.0) 
            GetOStream(Warnings1) << "WARNING: RAPShiftFactory has detected a negative shift... This implies a less stable coarse grid."<<std::endl;
          shift = Teuchos::as<Scalar>(d_shift);
        }
        else {
          double base_shift = pL.get<double>("rap: shift");
          if(level == 1) shift = Teuchos::as<Scalar>(base_shift);
          else shift = Teuchos::as<Scalar>(pow(base_shift,level) - pow(base_shift,level-1));
        }
      }
      GetOStream(Runtime0) << "RAPShiftFactory: Using shift " << shift << std::endl;

        
      // recombine to get K+shift*M
      {
  SubFactoryMonitor m2(*this, "Add: RKP + s*RMP", coarseLevel);
  Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::TwoMatrixAdd(*Kc, false, Teuchos::ScalarTraits<Scalar>::one(), *Mc, false, shift, Ac, GetOStream(Statistics2));
  Ac->fillComplete();
      }

      Teuchos::ArrayView<const double> relativeFloor = pL.get<Teuchos::Array<double> >("rap: relative diagonal floor")();
      if(relativeFloor.size() > 0) 
        Xpetra::MatrixUtils<SC,LO,GO,NO>::RelativeDiagonalBoost(Ac, relativeFloor,GetOStream(Statistics2));
   

      bool repairZeroDiagonals = pL.get<bool>("RepairMainDiagonal") || pL.get<bool>("rap: fix zero diagonals");
      bool checkAc             = pL.get<bool>("CheckMainDiagonal")|| pL.get<bool>("rap: fix zero diagonals"); ;
      if (checkAc || repairZeroDiagonals)
        Xpetra::MatrixUtils<SC,LO,GO,NO>::CheckRepairMainDiagonal(Ac, repairZeroDiagonals, GetOStream(Warnings1));

      RCP<ParameterList> params = rcp(new ParameterList());;
      params->set("printLoadBalancingInfo", true);
      GetOStream(Statistics0) << PerfUtils::PrintMatrixInfo(*Ac, "Ac", params);

      Set(coarseLevel, "A", Ac);
      // We only need K in the 'high storage' mode
      if(!use_low_storage)
        Set(coarseLevel, "K", Kc);

      if(!M_is_diagonal) {
        Set(coarseLevel, "M", Mc);
      }
      else {
        // If M is diagonal, then we only pass that part down the hierarchy
        // NOTE: Should we be doing some kind of rowsum instead?
        RCP<Vector> Mcv = Xpetra::VectorFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(Mc->getRowMap(),false);
        Mc->getLocalDiagCopy(*Mcv);
        Set(coarseLevel, "Mdiag", Mcv);
      }

      //      Set(coarseLevel, "RAP Pattern", Ac);
    }

    if (transferFacts_.begin() != transferFacts_.end()) {
      SubFactoryMonitor m(*this, "Projections", coarseLevel);

      // call Build of all user-given transfer factories
      for (std::vector<RCP<const FactoryBase> >::const_iterator it = transferFacts_.begin(); it != transferFacts_.end(); ++it) {
        RCP<const FactoryBase> fac = *it;
        GetOStream(Runtime0) << "RAPShiftFactory: call transfer factory: " << fac->description() << std::endl;
        fac->CallBuild(coarseLevel);
        // AP (11/11/13): I am not sure exactly why we need to call Release, but we do need it to get rid
        // of dangling data for CoordinatesTransferFactory
        coarseLevel.Release(*fac);
      }
    }
  }

  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
  void RAPShiftFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::AddTransferFactory(const RCP<const FactoryBase>& factory) {
    // check if it's a TwoLevelFactoryBase based transfer factory
    TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::rcp_dynamic_cast<const TwoLevelFactoryBase>(factory) == Teuchos::null, Exceptions::BadCast, "MueLu::RAPShiftFactory::AddTransferFactory: Transfer factory is not derived from TwoLevelFactoryBase. This is very strange. (Note: you can remove this exception if there's a good reason for)");
    transferFacts_.push_back(factory);
  }

} //namespace MueLu

#define MUELU_RAPSHIFTFACTORY_SHORT
#endif // MUELU_RAPSHIFTFACTORY_DEF_HPP