schursolveu.h

Go to the documentation of this file.

// ========================================================================= //
// This file is part of MyraMath, copyright (c) 2014-2019 by Ryan A Chilton  //
// and distributed by MyraCore, LLC. See LICENSE.txt for license terms.      //
// ========================================================================= //

#ifndef MYRAMATH_MULTIFRONTAL_LU_SCHURSOLVEU_H
#define MYRAMATH_MULTIFRONTAL_LU_SCHURSOLVEU_H

// Note that this actually solves by U'

#include <myramath/multifrontal/detail/schursolve.h>

namespace myra {
namespace multifrontal {
namespace lu {
namespace schursolveu {

template<class Number> class JobGraphBase : public ::myra::multifrontal::detail::schursolve::JobGraphBase2<Number>
  {
  public:

    // Structural typedefs for base class.
    typedef LUKernel<Number> Kernel;
    typedef ::myra::multifrontal::detail::schursolve::JobGraphBase2<Number> Base;

    // Typedefs related to numeric storage.
    typedef ::myra::multifrontal::detail::lu::LUContainer<Kernel> LUContainer;
    typedef ::myra::multifrontal::detail::XContainer<Number> XContainer;
    typedef ::myra::multifrontal::detail::XContributor<Number> XContributor;
    
    // Constructor - requires const reference to LuContainer, mutable reference to XContainer
    JobGraphBase(const LUContainer* in_lucontainer, XContainer* in_xcontainer, const XContributor* in_xcontributor)
      : Base(&in_xcontainer->tree(), in_xcontainer, in_xcontributor), lucontainer(in_lucontainer), xcontainer(in_xcontainer), xcontributor(in_xcontributor) { }

    // Virtual copy constructor.
    virtual ::myra::JobGraphBase* clone() const
      { return new JobGraphBase(*this); }

  protected:

    // Accessors for L Kernel's / Block's
    const Kernel& u(int n, int ij) const
      { return lucontainer->lu(this->s2a(n),ij); }
    CMatrixRange<Number> u(int n, int i, int j) const
      { return lucontainer->lu(this->s2a(n),i,j); }
  
    // Solves Un(k,k)'*Xn(k,j) = Bn(k,j)
    virtual uint64_t backsolve(int n, int k, int j)
      {
      // Assign internode contributions on first k-iteration.
      if (k == 0) this->assign_contributions(n,k,j);
      MatrixRange<Number> Bn_kj = this->b(n,k,j);
      const Kernel& Un_kk = this->u(n,k);
      return Un_kk.solveU(Bn_kj,'L','T');
      }

    // Downdates Bn(i,j) -= Un(k,i)'*Xn(k,j)
    virtual uint64_t downdate(int n, int k, int i, int j) 
      {
      // Useful constants.
      Number one(1);
      Number zero(0);
      // Downdate Bn(i,j) -= Un(k,i)'*Xn(k,j) [gemm]
       MatrixRange<Number> Bn_ij = this->b(n,i,j);
      CMatrixRange<Number> Un_ki = this->u(n,k,i);
      CMatrixRange<Number> Xn_kj = this->x(n,k,j);
      // Carefully choose beta and order of steps to avoid workspace initialiation.
      Number beta = k ? one : zero;
      uint64_t w = gemm_nwork(Bn_ij, Un_ki, 'T', Xn_kj, 'N', -one, beta);
      // Add internode contributions if on first k-iteration.
      if (k == 0) this->add_contributions(n,i,j);
      return w;
      }

  private:
  
    // Reference to LUContainer (const).
    const LUContainer* lucontainer;

    // Reference to XContainer (mutable).
    XContainer* xcontainer;
    const XContributor* xcontributor;

  }; // class JobGraph

} } } } // namespace

#endif