Teko_BlockLowerTriInverseOp.cpp

// @HEADER
// *****************************************************************************
//      Teko: A package for block and physics based preconditioning
//
// Copyright 2010 NTESS and the Teko contributors.
// SPDX-License-Identifier: BSD-3-Clause
// *****************************************************************************
// @HEADER
 
#include "Teko_BlockLowerTriInverseOp.hpp"
 
#include "Teuchos_Utils.hpp"
 
namespace Teko {
 
using Teuchos::RCP;
 
BlockLowerTriInverseOp::BlockLowerTriInverseOp(BlockedLinearOp& L,
                                               const std::vector<LinearOp>& invDiag)
    : L_(L) {
  invDiag_ = invDiag;
 
  // sanity check
  int blocks = blockRowCount(L_);
  TEUCHOS_ASSERT(blocks > 0);
  TEUCHOS_ASSERT(blocks == blockColCount(L_));
  TEUCHOS_ASSERT(blocks == (int)invDiag_.size());
 
  // create the range and product space
 
  // just flip flop them!
  productRange_  = L_->productDomain();
  productDomain_ = L_->productRange();
}
 
namespace {
bool compatibleMultiVectors(const BlockedMultiVector& a, const BlockedMultiVector& b) {
  if (!a || !b) return false;
 
  return a->domain()->dim() == b->domain()->dim();
}
}  // namespace
 
void BlockLowerTriInverseOp::implicitApply(const BlockedMultiVector& src, BlockedMultiVector& dst,
                                           const double alpha, const double beta) const {
  int blocks = blockCount(src);
 
  TEUCHOS_ASSERT(blocks == blockRowCount(L_));
  TEUCHOS_ASSERT(blocks == blockCount(dst));
 
  bool needsAlloc = !allocated;
  needsAlloc |= !compatibleMultiVectors(srcScrap_, src);
  needsAlloc |= !compatibleMultiVectors(dstScrap_, dst);
  if (needsAlloc) {
    srcScrap_ = deepcopy(src);
    dstScrap_ = deepcopy(dst);
    allocated = true;
  }
 
  // build a scrap vector for storing work
  Thyra::assign<double>(srcScrap_.ptr(), *src);
  BlockedMultiVector dstCopy;
  if (beta != 0.0) {
    Thyra::assign<double>(dstScrap_.ptr(), *dst);
    dstCopy = dstScrap_;
  } else
    dstCopy = dst;  // shallow copy
 
  // extract the blocks componets from
  // the source and destination vectors
  std::vector<MultiVector> dstVec;
  std::vector<MultiVector> scrapVec;
  for (int b = 0; b < blocks; b++) {
    dstVec.push_back(getBlock(b, dstCopy));
    scrapVec.push_back(getBlock(b, srcScrap_));
  }
 
  // run forward-substituion: run over each column
  //    From Heath pg. 65
  for (int b = 0; b < blocks; b++) {
    applyOp(invDiag_[b], scrapVec[b], dstVec[b]);
 
    // loop over each row
    for (int i = b + 1; i < blocks; i++) {
      LinearOp u_ib = getBlock(i, b, L_);
      if (u_ib != Teuchos::null) {
        applyOp(u_ib, dstVec[b], scrapVec[i], -1.0, 1.0);
      }
    }
  }
 
  // scale result by alpha
  if (beta != 0)
    update(alpha, dstCopy, beta, dst);  // dst = alpha * dstCopy + beta * dst
  else if (alpha != 1.0)
    scale(alpha, dst);  // dst = alpha * dst
}
 
void BlockLowerTriInverseOp::describe(Teuchos::FancyOStream& out_arg,
                                      const Teuchos::EVerbosityLevel verbLevel) const {
  using Teuchos::OSTab;
 
  RCP<Teuchos::FancyOStream> out = rcp(&out_arg, false);
  OSTab tab(out);
  switch (verbLevel) {
    case Teuchos::VERB_DEFAULT:
    case Teuchos::VERB_LOW: *out << this->description() << std::endl; break;
    case Teuchos::VERB_MEDIUM:
    case Teuchos::VERB_HIGH:
    case Teuchos::VERB_EXTREME: {
      *out << Teuchos::Describable::description() << "{"
           << "rangeDim=" << this->range()->dim() << ",domainDim=" << this->domain()->dim()
           << ",rows=" << blockRowCount(L_) << ",cols=" << blockColCount(L_) << "}\n";
      {
        OSTab tab2(out);
        *out << "[L Operator] = ";
        *out << Teuchos::describe(*L_, verbLevel);
      }
      {
        OSTab tab2(out);
        *out << "[invDiag Operators]:\n";
        tab.incrTab();
        for (int i = 0; i < blockRowCount(L_); i++) {
          *out << "[invD(" << i << ")] = ";
          *out << Teuchos::describe(*invDiag_[i], verbLevel);
        }
      }
      break;
    }
    default: TEUCHOS_TEST_FOR_EXCEPT(true);  // Should never get here!
  }
}
 
}  // end namespace Teko