Teko_LU2x2DiagonalStrategy.cpp

/*
// @HEADER
//
// ***********************************************************************
//
//      Teko: A package for block and physics based preconditioning
//                  Copyright 2010 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 Eric C. Cyr (eccyr@sandia.gov)
//
// ***********************************************************************
//
// @HEADER

*/

#include "Teko_LU2x2DiagonalStrategy.hpp"

#include "Teuchos_TimeMonitor.hpp"

namespace Teko {

using Teuchos::TimeMonitor;

Teuchos::RCP<Teuchos::Time> LU2x2DiagonalStrategy::initTimer_;
Teuchos::RCP<Teuchos::Time> LU2x2DiagonalStrategy::invSTimer_;
Teuchos::RCP<Teuchos::Time> LU2x2DiagonalStrategy::invA00Timer_;
Teuchos::RCP<Teuchos::Time> LU2x2DiagonalStrategy::opsTimer_;

void LU2x2DiagonalStrategy::buildTimers() {
  if (initTimer_ == Teuchos::null)
    initTimer_ = TimeMonitor::getNewTimer("LU2x2DiagonalStrategy::initializePrec");

  if (invSTimer_ == Teuchos::null)
    invSTimer_ = TimeMonitor::getNewTimer("LU2x2DiagonalStrategy::initializePrec invS");

  if (invA00Timer_ == Teuchos::null)
    invA00Timer_ = TimeMonitor::getNewTimer("LU2x2DiagonalStrategy::initializePrec invA00");

  if (opsTimer_ == Teuchos::null)
    opsTimer_ = TimeMonitor::getNewTimer("LU2x2DiagonalStrategy::initializePrec buildOps");
}

LU2x2DiagonalStrategy::LU2x2DiagonalStrategy() : a00InverseType_(Diagonal) { buildTimers(); }

LU2x2DiagonalStrategy::LU2x2DiagonalStrategy(const Teuchos::RCP<InverseFactory>& invFA,
                                             const Teuchos::RCP<InverseFactory>& invS)
    : invFactoryA00_(invFA), invFactoryS_(invS), a00InverseType_(Diagonal) {
  buildTimers();
}

const Teko::LinearOp LU2x2DiagonalStrategy::getHatInvA00(const Teko::BlockedLinearOp& A,
                                                         BlockPreconditionerState& state) const {
  initializeState(A, state);

  return state.getModifiableOp("invA00");
}

const Teko::LinearOp LU2x2DiagonalStrategy::getTildeInvA00(const Teko::BlockedLinearOp& A,
                                                           BlockPreconditionerState& state) const {
  initializeState(A, state);

  return state.getModifiableOp("invA00");
}

const Teko::LinearOp LU2x2DiagonalStrategy::getInvS(const Teko::BlockedLinearOp& A,
                                                    BlockPreconditionerState& state) const {
  initializeState(A, state);

  return state.getModifiableOp("invS");
}

void LU2x2DiagonalStrategy::initializeState(const Teko::BlockedLinearOp& A,
                                            BlockPreconditionerState& state) const {
  Teko_DEBUG_SCOPE("LU2x2DiagonalStrategy::initializeState", 10);

  // no work to be done
  if (state.isInitialized()) return;

  Teuchos::TimeMonitor timer(*initTimer_, true);

  // extract sub blocks
  LinearOp A00 = Teko::getBlock(0, 0, A);
  LinearOp A01 = Teko::getBlock(0, 1, A);
  LinearOp A10 = Teko::getBlock(1, 0, A);
  LinearOp A11 = Teko::getBlock(1, 1, A);

  // build the Schur complement
  ModifiableLinearOp& S = state.getModifiableOp("S");
  {
    Teko_DEBUG_SCOPE("Building S", 5);
    Teuchos::TimeMonitor timerS(*opsTimer_, true);

    LinearOp diagA00 = getInvDiagonalOp(A00, a00InverseType_);

    // build Schur-complement
    ModifiableLinearOp& triple = state.getModifiableOp("triple");
    triple                     = explicitMultiply(A10, diagA00, A01, triple);
    S                          = explicitAdd(scale(-1.0, A11), triple, S);
  }

  // build inverse S
  {
    Teko_DEBUG_SCOPE("Building inverse(S)", 5);
    Teuchos::TimeMonitor timerInvS(*invSTimer_, true);

    ModifiableLinearOp& invS = state.getModifiableOp("invS");
    if (invS == Teuchos::null)
      invS = buildInverse(*invFactoryS_, S);
    else
      rebuildInverse(*invFactoryS_, S, invS);
  }

  // build inverse A00
  {
    Teko_DEBUG_SCOPE("Building inverse(A00)", 5);
    Teuchos::TimeMonitor timerInvA00(*invA00Timer_, true);

    ModifiableLinearOp& invA00 = state.getModifiableOp("invA00");
    *getOutputStream() << "(LU2x2) invA00 pointer = " << invA00 << std::endl;
    if (invA00 == Teuchos::null)
      invA00 = buildInverse(*invFactoryA00_, A00);
    else
      rebuildInverse(*invFactoryA00_, A00, invA00);
  }

  // mark state as initialized
  state.setInitialized(true);
}

void LU2x2DiagonalStrategy::initializeFromParameterList(const Teuchos::ParameterList& pl,
                                                        const InverseLibrary& invLib) {
  Teko_DEBUG_SCOPE("LU2x2DiagonalStrategy::initializeFromParameterList", 10);

  std::string invStr = "", invA00Str = "", invSStr = "";
#if defined(Teko_ENABLE_Amesos)
  invStr = "Amesos";
#elif defined(Teko_ENABLE_Amesos2)
  invStr = "Amesos2";
#endif

  // "parse" the parameter list
  if (pl.isParameter("Inverse Type")) invStr = pl.get<std::string>("Inverse Type");
  if (pl.isParameter("Inverse A00 Type")) invA00Str = pl.get<std::string>("Inverse A00 Type");
  if (pl.isParameter("Inverse Schur Type")) invSStr = pl.get<std::string>("Inverse Schur Type");
  if (pl.isParameter("Diagonal Type")) {
    std::string massInverseStr = pl.get<std::string>("Diagonal Type");

    // build inverse types
    a00InverseType_ = getDiagonalType(massInverseStr);
  }

  // set defaults as needed
  if (invA00Str == "") invA00Str = invStr;
  if (invSStr == "") invSStr = invStr;

  Teko_DEBUG_MSG_BEGIN(5) DEBUG_STREAM << "LU2x2 Diagonal Strategy Parameters: " << std::endl;
  DEBUG_STREAM << "   inv type   = \"" << invStr << "\"" << std::endl;
  DEBUG_STREAM << "   inv A00 type = \"" << invA00Str << "\"" << std::endl;
  DEBUG_STREAM << "   inv S type = \"" << invSStr << "\"" << std::endl;
  DEBUG_STREAM << "LU2x2 Diagonal Strategy Parameter list: " << std::endl;
  pl.print(DEBUG_STREAM);
  Teko_DEBUG_MSG_END()

      // build velocity inverse factory
      invFactoryA00_ = invLib.getInverseFactory(invA00Str);

  if (invA00Str == invSStr)
    invFactoryS_ = invFactoryA00_;
  else
    invFactoryS_ = invLib.getInverseFactory(invSStr);
}
}  // end namespace Teko