Teko_DiagnosticLinearOp.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_DiagnosticLinearOp.hpp"

#include "Teuchos_TimeMonitor.hpp"

#include "Thyra_MultiVectorStdOps.hpp"

namespace Teko {

DiagnosticLinearOp::DiagnosticLinearOp(const Teuchos::RCP<std::ostream>& ostrm,
                                       const ModifiableLinearOp& A,
                                       const std::string& diagnosticString)
    : outputStream_(ostrm),
      wrapOpA_(A),
      wrapOpA_lo_(A),
      diagString_(diagnosticString),
      timer_(diagnosticString) {}

DiagnosticLinearOp::DiagnosticLinearOp(const Teuchos::RCP<std::ostream>& ostrm, const LinearOp& A,
                                       const std::string& diagnosticString)
    : outputStream_(ostrm),
      wrapOpA_(Teuchos::null),
      wrapOpA_lo_(A),
      diagString_(diagnosticString),
      timer_(diagnosticString) {}

DiagnosticLinearOp::DiagnosticLinearOp(const Teuchos::RCP<std::ostream>& ostrm,
                                       const LinearOp& fwdOp, const ModifiableLinearOp& A,
                                       const std::string& diagnosticString)
    : outputStream_(ostrm),
      wrapOpA_(A),
      wrapOpA_lo_(A),
      fwdOp_(fwdOp),
      diagString_(diagnosticString),
      timer_(diagnosticString) {}

DiagnosticLinearOp::~DiagnosticLinearOp() {
  double elapsedTime = totalTime();
  int applications   = numApplications();

  (*outputStream_) << "DiagnosticLinearOp \"" << diagString_ << "\": "
                   << "elapsed = " << elapsedTime << ", "
                   << "applications = " << applications << ", ";
  if (applications > 0)
    (*outputStream_) << "timer/app = " << elapsedTime / double(applications) << std::endl;
  else
    (*outputStream_) << "timer/app = "
                     << "none" << std::endl;
}

void DiagnosticLinearOp::implicitApply(const MultiVector& x, MultiVector& y, const double alpha,
                                       const double beta) const {
  Teko_DEBUG_SCOPE("DiagnosticLinearOp::implicityApply", 10);

  // start timer on construction, end on destruction
  Teuchos::TimeMonitor monitor(timer_, false);

  MultiVector z;  // for temporary storage dealing with nozero beta
  if (beta != 0.0) z = deepcopy(y);

  wrapOpA_lo_->apply(Thyra::NOTRANS, *x, y.ptr(), alpha, beta);

  // print residual if there is a fwd Op
  bool printResidual = (fwdOp_ != Teuchos::null);
  if (printResidual) {
    // compute residual
    MultiVector residual = Teko::deepcopy(x);
    // fwdOp_->apply(Thyra::NOTRANS,*y,residual.ptr(),-1.0,1.0);

    fwdOp_->apply(Thyra::NOTRANS, *y, residual.ptr(), -1.0, alpha);
    if (beta != 0.0) fwdOp_->apply(Thyra::NOTRANS, *z, residual.ptr(), beta, 1.0);

    // calculate norms
    std::vector<double> norms(y->domain()->dim());     // size of column count
    std::vector<double> rhsNorms(x->domain()->dim());  // size of column count
    Thyra::norms_2<double>(*residual, Teuchos::arrayViewFromVector(norms));
    Thyra::norms_2<double>(*x, Teuchos::arrayViewFromVector(rhsNorms));

    // print out residual norms
    (*outputStream_) << "DiagnosticLinearOp \"" << diagString_ << "\": residual = [";
    for (std::size_t i = 0; i < norms.size(); ++i)
      (*outputStream_) << " " << std::scientific << std::setprecision(4)
                       << norms[i] / rhsNorms[i];  // << " (" <<rhsNorms[i]<<") ";
    (*outputStream_) << " ]" << std::endl;

    residualNorm_ = norms[0];
  }
}

}  // end namespace Teko