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_RLDLT_FACTOR_H
#define MYRAMATH_MULTIFRONTAL_RLDLT_FACTOR_H

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

#include <myramath/multifrontal/rldlt/Kernel.h>

#include <myramath/dense/Matrix.h>
#include <myramath/dense/MatrixRange.h>
#include <myramath/dense/LowerMatrix.h>
#include <myramath/dense/LowerMatrixRange.h>
#include <myramath/dense/trsm.h>
#include <myramath/dense/syrk.h>
#include <myramath/dense/gemm.h>
#include <myramath/dense/detail/nwork.h>

namespace myra {
namespace multifrontal {
namespace rldlt {
namespace factor {

template<class Precision> class JobGraphBase : public ::myra::multifrontal::detail::llt::factor::JobGraphBase2<RLDLTKernel<Precision> >
  {
  public:

    // Structural typedefs for base class.
    typedef Precision Number;
    typedef RLDLTKernel<Precision> Kernel;
    typedef ::myra::multifrontal::detail::llt::factor::JobGraphBase2<Kernel> Base;

    // Typedefs related to numeric storage on LContainer
    typedef ::myra::multifrontal::detail::llt::LContainer<Kernel> LContainer;
    typedef LowerMatrix<Number> Diagonal;

    // Constructor - requires reference to LContainer under construction.
    JobGraphBase(LContainer* in_lcontainer)
       : Base(in_lcontainer), lcontainer(in_lcontainer) { }

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

  protected:
   
    // Factors An(k,k) = Ln(k,k)*Ln(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); 
      return lcontainer->factor(n,k);
      }

    // Solves Ln(i,k)*Ln(k,k)' = An(i,k)
    virtual uint64_t trsm(int n, int k, int i)
      {
      // Assign internode contributions if on first k-iteration.
      if (k == 0) this->assign_contributions(n,i,k); 
      MatrixRange<Number> An_ik = this->a(n,i,k);
      const Kernel& Ln_kk = this->l(n,k);
      uint64_t w = Ln_kk.solveL(An_ik,'R','T');
      Ln_kk.solveI(An_ik,'R');
      return w;
      }

    // Downdates An(ij,ij) -= Ln(ij,k)*Ln(ij,k)'
    virtual uint64_t rankk(int n, int k, int ij)
      {
      // Useful constants.
      Number one(1);
      Number zero(0);
      // Downdate An(ij,ij) -= Ln(ij,k)*Ln(ij,k)' [syrk]
      LowerMatrixRange<Number> An_ij = this->a(n,ij);
      CMatrixRange<Number> Ln_ijk = this->l(n,ij,k);
      // Carefully choose beta and order of steps to avoid workspace initialization.
      Number beta = k ? one : zero;
      // Two syrk's required, with alpha = +/- one.
      std::pair<int,int> inertia = this->l(n,k).inertia();
      uint64_t w = 0;
      w += syrk_nwork(An_ij, Ln_ijk. left(inertia.first) , 'N', -one, beta);
      w += syrk_nwork(An_ij, Ln_ijk.right(inertia.second), 'N',  one, one );
      // Add internode contributions if on first k-iteration.
      if (k == 0) this->add_contributions(n,ij);
      return w;
      }
    
    // Downdates An(i,j) -= Ln(i,k)*Ln(j,k)'
    virtual uint64_t gemm(int n, int k, int i, int j)
      {
      // Useful constants.
      Number one(1);
      Number zero(0);
      // Downdate An(i,j) -= Ln(i,k)*Ln(j,k)' [gemm]
      MatrixRange<Number>  An_ij = this->a(n,i,j);
      CMatrixRange<Number> Ln_ik = this->l(n,i,k);
      CMatrixRange<Number> Ln_jk = this->l(n,j,k);
      // Carefully choose beta and order of steps to avoid workspace initialization.
      Number beta = k ? one : zero;
      // Two gemm's required, with alpha = +/- one.
      std::pair<int,int> inertia = this->l(n,k).inertia();
      uint64_t w = 0;
      w += gemm_nwork(An_ij, Ln_ik. left(inertia.first ), 'N', Ln_jk. left(inertia.first ), 'T', -one, beta);
      w += gemm_nwork(An_ij, Ln_ik.right(inertia.second), 'N', Ln_jk.right(inertia.second), 'T',  one, one );
      // Add internode contributions if on first k-iteration.
      if (k == 0) this->add_contributions(n,i,j);
      return w;
      }

  private:
  
    // Reference to LContainer under construction.
    LContainer* lcontainer;

  }; // class JobGraph

} } } } // namespace

#endif