factor.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_FACTOR_H
#define MYRAMATH_MULTIFRONTAL_LU_FACTOR_H

#include <myramath/multifrontal/detail/lu/factor.h>

namespace myra {
namespace multifrontal {
namespace lu {
namespace factor {

template<class Number> class JobGraphBase : public ::myra::multifrontal::detail::lu::factor::JobGraphBase2<LUKernel<Number> >
  {
  public:

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

    // Typedefs related to numeric storage.
    typedef LUKernel<Number> Kernel;    
    typedef ::myra::multifrontal::detail::lu::LUContainer<Kernel> LUContainer;

    // Constructor - requires reference to LUContainer under construction.
    JobGraphBase(LUContainer* in_lucontainer)
       : Base(in_lucontainer), lucontainer(in_lucontainer) { }

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

  protected:

    // Factors An(k,k) = Ln(k,k)*Un(k,k)
    virtual uint64_t factor(int n, int k)
      {
      // Assign internode contributions if on first k-iteration.
      if (k == 0) this->assign_contributions(n,k,k);
      return lucontainer->factor(n,k); 
      }

    // Solves Ln(i,k)*Un(k,k) = An(i,k)
    virtual uint64_t backsolveL(int n, int k, int i)
      {
      // Assign internode contributions if on first k-iteration.
      if (k == 0) this->assign_contributions(n,i,k);
      const Kernel& Un_kk = this->lu(n,k);
      MatrixRange<Number> An_ik = this->a(n,i,k);
      return Un_kk.solveU(An_ik,'R','N');
      }

    // Solves Ln(k,k)*Un(k,j) = An(k,j)
    virtual uint64_t backsolveU(int n, int k, int j)
      {
      // Assign internode contributions if on first k-iteration.
      if (k == 0) this->assign_contributions(n,k,j);
      const Kernel& Ln_kk = this->lu(n,k);
      MatrixRange<Number> An_kj = this->a(n,k,j);
      return Ln_kk.solveL(An_kj,'L','N');
      }

    // Downdates An(i,j) -= Ln(i,k)*Un(k,j)
    virtual uint64_t downdate(int n, int k, int i, int j)
      {
      // Useful constants.
      Number one(1);
      Number zero(0);
      // Downdate An(i,j) -= Ln(i,k)*Un(k,j) [gemm]
      MatrixRange<Number> An_ij = this->a(n,i,j);
      CMatrixRange<Number> Ln_ik = this->lu(n,i,k);
      CMatrixRange<Number> Un_kj = this->lu(n,k,j);
      // Carefully choose beta and order of steps to avoid workspace initialization.
      Number beta = k ? one : zero;
      uint64_t w = gemm_nwork(An_ij, Ln_ik, 'N', Un_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 under construction.
    LUContainer* lucontainer;

  }; // class JobGraph

} } } } // namespace

#endif