Teko_DiagonalPreconditionerFactory.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_DiagonalPreconditionerFactory.hpp"
#include "Teko_DiagonalPreconditionerOp.hpp"
#include "Thyra_get_Epetra_Operator.hpp"
#include "Epetra_CrsMatrix.h"
#include "EpetraExt_PointToBlockDiagPermute.h"

#include "Teko_TpetraHelpers.hpp"
#include "Thyra_EpetraLinearOp.hpp"
#include "Thyra_TpetraLinearOp.hpp"

using Teuchos::rcp;
using Teuchos::RCP;

namespace Teko {

DiagonalPrecondState::DiagonalPrecondState() {}

/*****************************************************/

DiagonalPreconditionerFactory::DiagonalPreconditionerFactory() {}

RCP<PreconditionerState> DiagonalPreconditionerFactory::buildPreconditionerState() const {
  return Teuchos::rcp(new DiagonalPrecondState());
}

LinearOp DiagonalPreconditionerFactory::buildPreconditionerOperator(
    LinearOp& lo, PreconditionerState& state) const {
  if (diagonalType_ == BlkDiag) {
    // Sanity check the state
    DiagonalPrecondState& MyState = Teuchos::dyn_cast<DiagonalPrecondState>(state);

    TEUCHOS_TEST_FOR_EXCEPTION(TpetraHelpers::isTpetraLinearOp(lo), std::runtime_error,
                               "BlkDiag not implemented for Tpetra operators");

    // Get the underlying Epetra_CrsMatrix, if we have one
    Teuchos::RCP<const Epetra_Operator> eo = Thyra::get_Epetra_Operator(*lo);
    TEUCHOS_ASSERT(eo != Teuchos::null);
    Teuchos::RCP<const Epetra_CrsMatrix> MAT =
        Teuchos::rcp_dynamic_cast<const Epetra_CrsMatrix>(eo);
    TEUCHOS_ASSERT(MAT != Teuchos::null);

    // Create a new EpetraExt_PointToBlockDiagPermute for the state object, if we don't have one
    Teuchos::RCP<EpetraExt_PointToBlockDiagPermute> BDP;
    if (MyState.BDP_ == Teuchos::null) {
      BDP = Teuchos::rcp(new EpetraExt_PointToBlockDiagPermute(*MAT));
      BDP->SetParameters(List_);
      BDP->Compute();
      MyState.BDP_ = BDP;
    }

    RCP<Epetra_FECrsMatrix> Hcrs = rcp(MyState.BDP_->CreateFECrsMatrix());
    return Thyra::epetraLinearOp(Hcrs);

    // Build the LinearOp object  (NTS: swapping the range and domain)
    // LinearOp MyOp = Teuchos::rcp(new
    // DiagonalPreconditionerOp(MyState.BDP_,lo->domain(),lo->range()));
  }

  return getInvDiagonalOp(lo, diagonalType_);
}

void DiagonalPreconditionerFactory::initializeFromParameterList(const Teuchos::ParameterList& pl) {
  List_ = pl;

  diagonalType_ = BlkDiag;
  if (pl.isParameter("Diagonal Type")) {
    diagonalType_ = getDiagonalType(pl.get<std::string>("Diagonal Type"));
    TEUCHOS_TEST_FOR_EXCEPT(diagonalType_ == NotDiag);
  }

  if (diagonalType_ == BlkDiag) {
    // Reset default to invert mode if the user hasn't specified something else
    Teuchos::ParameterList& SubList = List_.sublist("blockdiagmatrix: list");
    SubList.set("apply mode", SubList.get("apply mode", "invert"));
  }
}

}  // end namespace Teko