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_ZLDLT_KERNEL_H
#define MYRAMATH_MULTIFRONTAL_ZLDLT_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/io/detail/complex.h>

#include <myramath/dense/MatrixRange.h>
#include <myramath/dense/LowerMatrixRange.h>
#include <myramath/dense/hetrf.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 MYRAMATH_EXPORT ZLDLTKernel
  {
  public:
    typedef std::complex<Precision> Number;

    explicit ZLDLTKernel() 
      { }
  
    explicit ZLDLTKernel(LowerMatrixRange<Number>& A) : L(A) 
      { }

    explicit ZLDLTKernel(LowerMatrixRange<Number>& 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<Number> result = sytrf_inplace(L);
      LDLSwaps<Number> 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 the 'R'ight?
    //    op    = Apply an operation to L? ('T'ranspose, 'H'ermitian, 'C'onjugate or 'N'othing)
    uint64_t solveL(const MatrixRange<Number>& 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.
      if (side == 'L')
        {
        // Check size.
        if (B.I != N) throw eprintf("ZLDLTKernel::solveL('L'eft), size mismatch B.I != N [%d != %d]", B.I, N);
        // Solve L*X = B?
        if (op == 'N')
          { 
          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;
          }
        // Solve transpose(L)*X = B?
        else if (op == 'T')
          { 
          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;
          }
        // Solve hermitian(L)*X = B?
        else if (op == 'H')
          { 
          iswap_rows(Q_swaps,B);
          R.solve(B,'L','H');
          uint64_t w = trsm_nwork('L','H',L,B); 
          iswap_rows(P_swaps,B);
          return w;
          }
        // Solve conjugate(L)*X = B?
        else if (op == 'C')
          { 
          swap_rows(P_swaps,B);
          uint64_t w = trsm_nwork('L','C',L,B);
          R.solve(B,'L','C'); 
          swap_rows(Q_swaps,B);
          return w;
          }
        else throw eprintf("ZLDLTKernel::solveL('L'eft), didn't understand op = %c", op);
        }
      else if (side == 'R')
        {
        // Check size.
        if (B.J != N) throw eprintf("ZLDLTKernel::solveL('R'ight), size mismatch B.J != N [%d != %d]", B.J, N);
        // Solve X*L = B?
        if (op == 'N')
          {
          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;
          }
        // Solve X*transpose(L) = B?
        else if (op == 'T')
          {
          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;
          }
        // Solve X*hermitian(L) = B?
        else if (op == 'H')
          {
          iswap_columns(P_swaps,B);
          uint64_t w = trsm_nwork('R','H',L,B); 
          R.solve(B,'R','H');
          iswap_columns(Q_swaps,B);
          return w;
          }
        // Solve X*conjugate(L) = B?
        else if (op == 'C')
          {
          swap_columns(Q_swaps,B);
          R.solve(B,'R','C'); 
          uint64_t w = trsm_nwork('R','C',L,B); 
          swap_columns(P_swaps,B); 
          return w;
          }
        else throw eprintf("ZLDLTKernel::solveL('R'ight), didn't understand op = %c", op);
        }
      else throw eprintf("ZLDLTKernel::solveL(), didn't understand side = %c", side);
      }

    void solveI(const MatrixRange<Number>& B, char side) const
      {
      if (side == 'L')
        B.bottom(n_minus) *= -Number(1);
      else if (side == 'R')
        B.right(n_minus) *= -Number(1);
      else throw eprintf("ZLDLTKernel::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<Number> 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<Number> R;

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

template<class Precision> class ReflectNumber< ZLDLTKernel<Precision> >
  { public: typedef std::complex<Precision> type; };

} // namespace myra

#endif