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_LU_KERNEL_H
#define MYRAMATH_MULTIFRONTAL_LU_KERNEL_H

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

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

#include <myramath/dense/MatrixRange.h>
#include <myramath/dense/getrf.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 MYRAMATH_EXPORT LUKernel
  {
  public:

    explicit LUKernel() 
      { }
  
    explicit LUKernel(MatrixRange<Number>& A) : LU(A) 
      { }

    explicit LUKernel(MatrixRange<Number>& A, InputStream& in) : LU(A) 
      { in  >> swaps; }

    void write(OutputStream& out) const 
      { out << swaps; }

    uint64_t factor()
      {
      swaps = getrf_inplace(LU);
      uint64_t k_work = this->size();
      return k_work*(k_work+1)*(2*k_work+1)/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();
      if (side == 'L')
        {
        // Check size.
        if (B.I != N) throw eprintf("LUKernel::solveL('L'eft), size mismatch B.I != N [%d != %d]", B.I, N);
        if (op == 'N') 
          {
          swap_rows(swaps,B);
          uint64_t w = trsm_nwork('L','L','N',LU,B,'U');
          return w;
          }
        else if (op == 'T') 
          {
          uint64_t w = trsm_nwork('L','L','T',LU,B,'U');
          iswap_rows(swaps,B);
          return w;
          }
        else if (op == 'H') 
          {
          uint64_t w = trsm_nwork('L','L','H',LU,B,'U');
          iswap_rows(swaps,B);
          return w;
          }
        else if (op == 'C') 
          {
          swap_rows(swaps,B);
          uint64_t w = trsm_nwork('L','L','C',LU,B,'U');
          return w;
          }
        else throw eprintf("LUKernel::solveL('L'eft), didn't understand op = %c", op);
        }
      else if (side == 'R')
        {
        // Check size.
        if (B.J != N) throw eprintf("LUKernel::solveL('R'ight), size mismatch B.J != N [%d != %d]", B.J, N);
        if (op == 'N') 
          {
          uint64_t w = trsm_nwork('R','L','N',LU,B,'U');
          swap_columns(swaps,B);
          return w;
          }
        else if (op == 'T') 
          { 
          iswap_columns(swaps,B);
          uint64_t w = trsm_nwork('R','L','T',LU,B,'U');
          return w;
          }
        else if (op == 'H') 
          {
          iswap_columns(swaps,B);
          uint64_t w = trsm_nwork('R','L','H',LU,B,'U');
          return w;
          }
        else if (op == 'C') 
          {
          uint64_t w = trsm_nwork('R','L','C',LU,B,'U');
          swap_columns(swaps,B);
          return w;
          }
        else throw eprintf("LUKernel::solveL('R'ight), didn't understand op = %c", op);
        }
      else throw eprintf("LUKernel::solveL(), didn't understand side = %c", side);
      }

    //    side  = Solve by U from the 'L'eft or from the 'R'ight?
    //    op    = Apply an operation to U? ('T'ranspose, 'H'ermitian, 'C'onjugate or 'N'othing)
    uint64_t solveU(const MatrixRange<Number>& B, char side, char op) const
      {
      int N = this->size();
      if (side == 'L')
        {
        // Check size.
        if (B.I != N) throw eprintf("LUKernel::solveU('L'eft), size mismatch B.I != N [%d != %d]", B.I, N);
        return trsm_nwork(side,'U',op,LU,B,'N'); 
        }
      else if (side == 'R')
        {
        // Check size.
        if (B.J != N) throw eprintf("LUKernel::solveU('R'ight), size mismatch B.J != N [%d != %d]", B.J, N);
        return trsm_nwork(side,'U',op,LU,B,'N');
        }
      else throw eprintf("LUKernel::solveU(), didn't understand side = %c", side);
      }

    int size() const
      { return LU.size().first; }
   
  private:
  
    // Points to underlying data.
    MatrixRange<Number> LU;

    // Pivoting metadata - encodes row swaps.
    std::vector<int> swaps;
  };

template<class Number> class ReflectNumber< LUKernel<Number> >
  { public: typedef Number type; };

} // namespace myra

#endif