MueLu_SchurComplementFactory_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_SCHURCOMPLEMENTFACTORY_DEF_HPP_
#define MUELU_SCHURCOMPLEMENTFACTORY_DEF_HPP_

#include <Xpetra_BlockedCrsMatrix.hpp>
#include <Xpetra_MultiVectorFactory.hpp>
#include <Xpetra_VectorFactory.hpp>
#include <Xpetra_MatrixFactory.hpp>
#include <Xpetra_Matrix.hpp>
#include <Xpetra_MatrixMatrix.hpp>
#include <Xpetra_TripleMatrixMultiply.hpp>
#include <Xpetra_CrsMatrixWrap.hpp>
#include <Xpetra_BlockedCrsMatrix.hpp>
#include <Xpetra_CrsMatrix.hpp>

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

namespace MueLu {

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

    const SC one = Teuchos::ScalarTraits<SC>::one();

    validParamList->set<RCP<const FactoryBase> >("A"    , NoFactory::getRCP(), "Generating factory of the matrix A used for building Schur complement (must be a 2x2 blocked operator)");
    validParamList->set<RCP<const FactoryBase> >("Ainv" , Teuchos::null,       "Generating factory of the inverse matrix used in the Schur complement");

    validParamList->set<SC>                     ("omega", one, "Scaling parameter in S = A(1,1) - 1/omega A(1,0) Ainv A(0,1)");

    return validParamList;
  }

  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
  void SchurComplementFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::DeclareInput(Level& currentLevel) const {
    Input(currentLevel, "A");
    
    // Get default or user-given inverse approximation factory
    RCP<const FactoryBase> AinvFact = GetFactory("Ainv");
    currentLevel.DeclareInput("Ainv", AinvFact.get(), this);
  }

  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
  void SchurComplementFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(Level& currentLevel) const {
    FactoryMonitor m(*this, "Build", currentLevel);

    RCP<Matrix> A = Get<RCP<Matrix> >(currentLevel, "A");
    RCP<BlockedCrsMatrix> bA = rcp_dynamic_cast<BlockedCrsMatrix>(A);

    TEUCHOS_TEST_FOR_EXCEPTION(bA.is_null(), Exceptions::BadCast,
                               "MueLu::SchurComplementFactory::Build: input matrix A is not of type BlockedCrsMatrix!");
    TEUCHOS_TEST_FOR_EXCEPTION(bA->Rows() != 2 || bA->Cols() != 2, Exceptions::RuntimeError,
                               "MueLu::SchurComplementFactory::Build: input matrix A is a " << bA->Rows() << "x" << bA->Cols() << " block matrix. We expect a 2x2 blocked operator.");

    // Calculate Schur Complement
    RCP<Matrix> Ainv = currentLevel.Get<RCP<Matrix> >("Ainv", this->GetFactory("Ainv").get());
    RCP<Matrix> S = ComputeSchurComplement(bA, Ainv);

    GetOStream(Statistics1) << "S has " << S->getGlobalNumRows() << "x" << S->getGlobalNumCols() << " rows and columns." << std::endl;

    // NOTE: "A" generated by this factory is actually the Schur complement
    // matrix, but it is required as all smoothers expect "A"
    Set(currentLevel, "A", S);
  }

  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
  RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
  SchurComplementFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::ComputeSchurComplement(RCP<BlockedCrsMatrix>& bA, RCP<Matrix>& Ainv) const {

    using STS = Teuchos::ScalarTraits<SC>;
    const SC zero = STS::zero(), one = STS::one();

    RCP<Matrix> A01 = bA->getMatrix(0,1);
    RCP<Matrix> A10 = bA->getMatrix(1,0);
    RCP<Matrix> A11 = bA->getMatrix(1,1);

    RCP<BlockedCrsMatrix> bA01 = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(A01);
    const bool isBlocked = (bA01 == Teuchos::null ? false : true);

    const ParameterList& pL = GetParameterList();
    const SC omega = pL.get<Scalar>("omega");

    TEUCHOS_TEST_FOR_EXCEPTION(omega == zero, Exceptions::RuntimeError,
                               "MueLu::SchurComplementFactory::Build: Scaling parameter omega must not be zero to avoid division by zero.");

    RCP<Matrix> S = Teuchos::null; // Schur complement
    RCP<Matrix> D = Teuchos::null; // temporary result for A10*Ainv*A01

    // only if the off-diagonal blocks A10 and A01 are non-zero we have to do the MM multiplication
    if(A01.is_null() == false && A10.is_null() == false) {
      // scale with -1/omega
      Ainv->scale(Teuchos::as<Scalar>(-one/omega));

      // build Schur complement operator
      if (!isBlocked) {
        RCP<ParameterList> myparams = rcp(new ParameterList);
        myparams->set("compute global constants", true);

        // -1/omega*Ainv*A01
        TEUCHOS_TEST_FOR_EXCEPTION(A01->getRangeMap()->isSameAs(*(Ainv->getDomainMap())) == false, Exceptions::RuntimeError,
                                   "MueLu::SchurComplementFactory::Build: RangeMap of A01 and domain map of Ainv are not the same.");
        RCP<Matrix> C = MatrixMatrix::Multiply(*Ainv, false, *A01, false, GetOStream(Statistics2), true, true, std::string("SchurComplementFactory"), myparams);

        // -1/omega*A10*Ainv*A01
        TEUCHOS_TEST_FOR_EXCEPTION(A01->getRangeMap()->isSameAs(*(A10->getDomainMap())) == false, Exceptions::RuntimeError,
                                   "MueLu::SchurComplementFactory::Build: RangeMap of A10 and domain map A01 are not the same.");
        D = MatrixMatrix::Multiply(*A10, false, *C, false, GetOStream(Statistics2), true, true, std::string("SchurComplementFactory"), myparams);
      }
      else {
        // nested blocking
        auto bA10  = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(A10);
        auto bAinv = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(Ainv);
        TEUCHOS_TEST_FOR_EXCEPTION(bAinv == Teuchos::null, Exceptions::RuntimeError,
                                   "MueLu::SchurComplementFactory::Build: Casting Ainv to BlockedCrsMatrix not possible.");

        // -1/omega*bAinv*bA01
        TEUCHOS_TEST_FOR_EXCEPTION(bA01->Rows() != bAinv->Cols(), Exceptions::RuntimeError,
                                   "MueLu::SchurComplementFactory::Build: Block rows and cols of bA01 and bAinv are not compatible.");
        RCP<BlockedCrsMatrix> C = MatrixMatrix::TwoMatrixMultiplyBlock(*bAinv, false, *bA01, false, GetOStream(Statistics2));

        // -1/omega*A10*Ainv*A01
        TEUCHOS_TEST_FOR_EXCEPTION(bA10->Rows() != bA01->Cols(), Exceptions::RuntimeError,
                                   "MueLu::SchurComplementFactory::Build: Block rows and cols of bA10 and bA01 are not compatible.");
        D = MatrixMatrix::TwoMatrixMultiplyBlock(*bA10, false, *C, false, GetOStream(Statistics2));
      }
      if (!A11.is_null()) {
        MatrixMatrix::TwoMatrixAdd(*A11, false, one, *D, false, one, S, GetOStream(Statistics2));
        S->fillComplete();

        TEUCHOS_TEST_FOR_EXCEPTION(A11->getRangeMap()->isSameAs(*(S->getRangeMap())) == false, Exceptions::RuntimeError,
                                   "MueLu::SchurComplementFactory::Build: RangeMap of A11 and S are not the same.");
        TEUCHOS_TEST_FOR_EXCEPTION(A11->getDomainMap()->isSameAs(*(S->getDomainMap())) == false, Exceptions::RuntimeError,
                                   "MueLu::SchurComplementFactory::Build: DomainMap of A11 and S are not the same.");
      }
      else {
        S = MatrixFactory::BuildCopy(D);
      }
    }
    else {
      if (!A11.is_null()) {
        S = MatrixFactory::BuildCopy(A11);
      } else {
        S = MatrixFactory::Build(A11->getRowMap(), 10 /*A11->getLocalMaxNumRowEntries()*/);
        S->fillComplete(A11->getDomainMap(),A11->getRangeMap());
      }
    }

    // Check whether Schur complement operator is a 1x1 block matrix.
    // If so, unwrap it and return the CrsMatrix based Matrix object
    // We need this, as single-block smoothers expect it this way.
    // In case of Thyra GIDs we obtain a Schur complement operator in Thyra GIDs
    // This may make some special handling in feeding the SchurComplement solver Apply routine
    // necessary!
    if (isBlocked) {
      RCP<BlockedCrsMatrix> bS = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(S);

      if (bS != Teuchos::null && bS->Rows() == 1 && bS->Cols() == 1) {
        RCP<Matrix> temp = bS->getCrsMatrix();
        S.swap(temp);
      }
    }

    return S;
  }

} // namespace MueLu

#endif /* MUELU_SCHURCOMPLEMENTFACTORY_DEF_HPP_ */