Kernel.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_KERNEL_H
#define MYRAMATH_MULTIFRONTAL_RLDLT_KERNEL_H

#include <myramath/utility/Number.h>
#include <myramath/utility/eprintf.h>

#include <myramath/io/Streams.h>
#include <myramath/io/ReaderWriter.h>
#include <myramath/io/detail/vector.h>

#include <myramath/dense/LowerMatrixRange.h>
#include <myramath/dense/sytrf.h>
#include <myramath/dense/trsm.h>
#include <myramath/dense/swaps.h>
#include <myramath/dense/detail/nwork.h>

#include <stdint.h>

namespace myra {

// Forward declarations.
class InputStream;
class OutputStream;
template<class Number> class MatrixRange;
template<class Number> class LowerMatrix;

template<class Precision> class LIBPUBLIC RLDLTKernel
  {
  public:

    explicit RLDLTKernel() 
      { }
  
    explicit RLDLTKernel(LowerMatrixRange<Precision>& A) : L(A) 
      { }

    explicit RLDLTKernel(LowerMatrixRange<Precision>& A, InputStream& in) : L(A) 
      { in  >> P_swaps >> Q_swaps >> R >> n_plus >> n_minus; }

    void write(OutputStream& out) const 
      { out << P_swaps << Q_swaps << R << n_plus << n_minus; }

    uint64_t factor()
      {
      // LDLSwaps<Precision> result = sytrf_inplace(L);
      LDLSwaps<Precision> result = sytrf_outplace(L);
      P_swaps = result.P_swaps;
      Q_swaps = result.Q_swaps;
      R = result.R;
      n_plus  = result.n_plus;  
      n_minus = result.n_minus;
      // Return flop count.
      uint64_t n_work = L.size();
      return n_work*(n_work+1)*(n_work+2)/6;
      }

    //    side  = Solve by L from the 'L'eft or from      std::cout << "exit solveL" << std::endl;       the 'R'ight?
    //    op    = Apply an operation to L? ('T'ranspose, 'H'ermitian, 'C'onjugate or 'N'othing)
    uint64_t solveL(const MatrixRange<Precision>& B, char side, char op) const
      {
      int N = this->size();  
      // Internally L is decomposed into P'*L*Q', so solving by it always takes a few steps.
      // Solve by L from the 'L'eft:
      if (side == 'L')
        {
        // Check size.
        if (B.I != N) throw eprintf("RLDLTKernel::solveL('L'eft), size mismatch B.I != N [%d != %d]", B.I, N);
        // Since L is real, L^N == L^C
        if (op == 'N' || op == 'C')
          {
          swap_rows(P_swaps,B); 
          uint64_t w = trsm_nwork('L','N',L,B); 
          R.solve(B,'L','N');
          swap_rows(Q_swaps,B);
          return w;
          }
        // Since L is real, L^T == L^H
        else if (op == 'T' || op == 'H')
          { 
          iswap_rows(Q_swaps,B);
          R.solve(B,'L','T');
          uint64_t w = trsm_nwork('L','T',L,B); 
          iswap_rows(P_swaps,B);
          return w;
          }
        else throw eprintf("RLDLTKernel::solveL('L'eft), didn't understand op = %c", op);
        }
      // Solve by L from the 'R'ight:
      else if (side == 'R')
        {
        // Check size.
        if (B.J != N) throw eprintf("RLDLTKernel::solveL('R'ight), size mismatch B.J != N [%d != %d]", B.J, N);
        // Since L is real, L^N == L^C
        if (op == 'N' || op == 'C')
          { 
          swap_columns(Q_swaps,B); 
          R.solve(B,'R','N');
          uint64_t w = trsm_nwork('R','N',L,B); 
          swap_columns(P_swaps,B);
          return w;
          }
        // Since L is real, L^T == L^H
        else if (op == 'T' || op == 'H')
          { 
          iswap_columns(P_swaps,B);
          uint64_t w = trsm_nwork('R','T',L,B);
          R.solve(B,'R','T');
          iswap_columns(Q_swaps,B);
          return w;
          }
        else throw eprintf("RLDLTKernel::solveL('R'ight), didn't understand op = %c", op);
        }
      else throw eprintf("RLDLTKernel::solveL(), didn't understand side = %c", side);
      }

    void solveI(const MatrixRange<Precision>& B, char side) const
      {
      if (side == 'L')
        B.bottom(n_minus) *= -Precision(1);
      else if (side == 'R')
        B.right(n_minus) *= -Precision(1);
      else throw eprintf("RLDLTKernel::solveI(), didn't understand side = %c", side);
      }

    std::pair<int,int> inertia() const
      { return std::pair<int,int>(n_plus, n_minus); }
      
    int size() const
      { return L.size(); }
  
  private:

    // Points to underlying data.
    LowerMatrixRange<Precision> L;

    // For applying permutation P.
    std::vector<int> P_swaps;

    // For applying permutation Q.
    std::vector<int> Q_swaps;

    // For applying pivot rotations R.
    PivotMatrix<Precision> R;

    // Encodes inertia.
    int n_plus;
    int n_minus;
    
  };

template<class Precision> class ReflectNumber< RLDLTKernel<Precision> >
  { public: typedef Precision type; };

} // namespace myra

#endif