Sacado_Fad_BLAS.hpp

Go to the documentation of this file.

// $Id$ 
// $Source$ 
// @HEADER
// ***********************************************************************
// 
//                           Sacado Package
//                 Copyright (2006) Sandia Corporation
// 
// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
// the U.S. Government retains certain rights in this software.
// 
// This library is free software; you can redistribute it and/or modify
// it under the terms of the GNU Lesser General Public License as
// published by the Free Software Foundation; either version 2.1 of the
// License, or (at your option) any later version.
//  
// This library is distributed in the hope that it will be useful, but
// WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Lesser General Public License for more details.
//  
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
// USA
// Questions? Contact David M. Gay (dmgay@sandia.gov) or Eric T. Phipps
// (etphipp@sandia.gov).
// 
// ***********************************************************************
// @HEADER

#ifndef SACADO_FAD_BLAS_HPP
#define SACADO_FAD_BLAS_HPP

#include "Teuchos_BLAS.hpp"
#include "Sacado_No_Kokkos.hpp"
#include "Sacado_CacheFad_DFad.hpp"
#include "Sacado_dummy_arg.hpp"

namespace Sacado {

  namespace Fad {

    template <typename OrdinalType, typename FadType>
    class ArrayTraits {

      typedef typename Sacado::ValueType<FadType>::type ValueType;
      typedef typename Sacado::ScalarType<FadType>::type scalar_type;
      typedef typename Sacado::dummy<ValueType,scalar_type>::type ScalarType;
      
    public:
      
      ArrayTraits(bool use_dynamic = true,
      OrdinalType workspace_size = 0);

      ArrayTraits(const ArrayTraits& a);

      ~ArrayTraits();
      
      void unpack(const FadType& a, OrdinalType& n_dot, ValueType& val, 
      const ValueType*& dot) const;
      
      void unpack(const FadType* a, OrdinalType n, OrdinalType inc,
      OrdinalType& n_dot, OrdinalType& inc_val, 
      OrdinalType& inc_dot,
      const ValueType*& val, const ValueType*& dot) const;
      
      void unpack(const FadType* A, OrdinalType m, OrdinalType n, 
      OrdinalType lda, OrdinalType& n_dot, 
      OrdinalType& lda_val, OrdinalType& lda_dot,
      const ValueType*& val, const ValueType*& dot) const;

      void unpack(const ValueType& a, OrdinalType& n_dot, ValueType& val, 
      const ValueType*& dot) const;
      
      void unpack(const ValueType* a, OrdinalType n, OrdinalType inc,
      OrdinalType& n_dot, OrdinalType& inc_val, 
      OrdinalType& inc_dot,
      const ValueType*& val, const ValueType*& dot) const;
      
      void unpack(const ValueType* A, OrdinalType m, OrdinalType n, 
      OrdinalType lda, OrdinalType& n_dot, 
      OrdinalType& lda_val, OrdinalType& lda_dot,
      const ValueType*& val, const ValueType*& dot) const;

      void unpack(const ScalarType& a, OrdinalType& n_dot, ScalarType& val, 
      const ScalarType*& dot) const;
      
      void unpack(const ScalarType* a, OrdinalType n, OrdinalType inc,
      OrdinalType& n_dot, OrdinalType& inc_val, 
      OrdinalType& inc_dot,
      const ScalarType*& val, const ScalarType*& dot) const;
      
      void unpack(const ScalarType* A, OrdinalType m, OrdinalType n, 
      OrdinalType lda, OrdinalType& n_dot, 
      OrdinalType& lda_val, OrdinalType& lda_dot,
      const ScalarType*& val, const ScalarType*& dot) const;

      void unpack(FadType& a, OrdinalType& n_dot, OrdinalType& final_n_dot, 
      ValueType& val, ValueType*& dot) const;
      
      void unpack(FadType* a, OrdinalType n, OrdinalType inc,
      OrdinalType& n_dot, OrdinalType& final_n_dot, 
      OrdinalType& inc_val, OrdinalType& inc_dot,
      ValueType*& val, ValueType*& dot) const;
      
      void unpack(FadType* A, OrdinalType m, OrdinalType n, OrdinalType lda, 
      OrdinalType& n_dot, OrdinalType& final_n_dot, 
      OrdinalType& lda_val, OrdinalType& lda_dot,
      ValueType*& val, ValueType*& dot) const;

      void pack(FadType& a, OrdinalType n_dot, const ValueType& val, 
    const ValueType* dot) const;
      
      void pack(FadType* a, OrdinalType n, OrdinalType inc,
    OrdinalType n_dot, OrdinalType inc_val, OrdinalType inc_dot,
    const ValueType* val, const ValueType* dot) const;
      
      void pack(FadType* A, OrdinalType m, OrdinalType n, 
    OrdinalType lda, OrdinalType n_dot, 
    OrdinalType lda_val, OrdinalType lda_dot,
    const ValueType* val, const ValueType* dot) const;

      void free(const FadType& a, OrdinalType n_dot, 
    const ValueType* dot) const;
      
      void free(const FadType* a, OrdinalType n, OrdinalType n_dot,
    OrdinalType inc_val, OrdinalType inc_dot,
    const ValueType* val, const ValueType* dot) const;

      void free(const FadType* A, OrdinalType m, OrdinalType n, 
    OrdinalType n_dot, OrdinalType lda_val, OrdinalType lda_dot,
    const ValueType* val, const ValueType* dot) const;

      void free(const ValueType& a, OrdinalType n_dot, 
    const ValueType* dot) const {}
      
      void free(const ValueType* a, OrdinalType n, OrdinalType n_dot,
    OrdinalType inc_val, OrdinalType inc_dot,
    const ValueType* val, const ValueType* dot) const {}

      void free(const ValueType* A, OrdinalType m, OrdinalType n, 
    OrdinalType n_dot, OrdinalType lda_val, OrdinalType lda_dot,
    const ValueType* val, const ValueType* dot) const {}

      void free(const ScalarType& a, OrdinalType n_dot, 
    const ScalarType* dot) const {}
      
      void free(const ScalarType* a, OrdinalType n, OrdinalType n_dot,
    OrdinalType inc_val, OrdinalType inc_dot,
    const ScalarType* val, const ScalarType* dot) const {}

      void free(const ScalarType* A, OrdinalType m, OrdinalType n, 
    OrdinalType n_dot, OrdinalType lda_val, OrdinalType lda_dot,
    const ScalarType* val, const ScalarType* dot) const {}

      ValueType* allocate_array(OrdinalType size) const;

      void free_array(const ValueType* ptr, OrdinalType size) const;

      bool is_array_contiguous(const FadType* a, OrdinalType n, 
             OrdinalType n_dot) const;

    protected:

      bool use_dynamic;

      OrdinalType workspace_size;

      mutable ValueType *workspace;

      mutable ValueType *workspace_pointer;
    
    };

    template <typename T> struct ArrayValueType { typedef T type; };

    template <typename OrdinalType, typename FadType>
    class BLAS : public Teuchos::DefaultBLASImpl<OrdinalType,FadType> {    
      
      typedef typename Teuchos::ScalarTraits<FadType>::magnitudeType MagnitudeType;
      typedef typename Sacado::ValueType<FadType>::type ValueType;
      typedef typename Sacado::ScalarType<FadType>::type scalar_type;
      typedef typename Sacado::dummy<ValueType,scalar_type>::type ScalarType;
      typedef Teuchos::DefaultBLASImpl<OrdinalType,FadType> BLASType;
    
    public:

    
      BLAS(bool use_default_impl = true,
     bool use_dynamic = true, OrdinalType static_workspace_size = 0);


      BLAS(const BLAS& x);

      virtual ~BLAS();

      

      
      void ROTG(FadType* da, FadType* db, MagnitudeType* c, FadType* s) const { 
  BLASType::ROTG(da,db,c,s); 
      }
      
      void ROT(const OrdinalType n, FadType* dx, const OrdinalType incx, 
         FadType* dy, const OrdinalType incy, MagnitudeType* c, 
         FadType* s) const { 
  BLASType::ROT(n,dx,incx,dy,incy,c,s); 
      }
      
      void SCAL(const OrdinalType n, const FadType& alpha, FadType* x, 
    const OrdinalType incx) const;

      void COPY(const OrdinalType n, const FadType* x, 
    const OrdinalType incx, FadType* y, 
    const OrdinalType incy) const;

      template <typename alpha_type, typename x_type>
      void AXPY(const OrdinalType n, const alpha_type& alpha, 
    const x_type* x, const OrdinalType incx, FadType* y, 
    const OrdinalType incy) const;

      typename Teuchos::ScalarTraits<FadType>::magnitudeType 
      ASUM(const OrdinalType n, const FadType* x, 
     const OrdinalType incx) const {
  return BLASType::ASUM(n,x,incx);
      }

      template <typename x_type, typename y_type>
      FadType DOT(const OrdinalType n, const x_type* x, 
      const OrdinalType incx, const y_type* y, 
      const OrdinalType incy) const;

      MagnitudeType NRM2(const OrdinalType n, const FadType* x, 
       const OrdinalType incx) const;

      OrdinalType IAMAX(const OrdinalType n, const FadType* x, 
      const OrdinalType incx) const {
  return BLASType::IAMAX(n,x,incx); 
      }

      

      
      template <typename alpha_type, typename A_type, typename x_type,
    typename beta_type>
      void GEMV(Teuchos::ETransp trans, const OrdinalType m, 
    const OrdinalType n, 
    const alpha_type& alpha, const A_type* A, 
    const OrdinalType lda, const x_type* x, 
    const OrdinalType incx, const beta_type& beta, 
    FadType* y, const OrdinalType incy) const;

      template <typename A_type>
      void TRMV(Teuchos::EUplo uplo, Teuchos::ETransp trans, 
    Teuchos::EDiag diag, const OrdinalType n, 
    const A_type* A, const OrdinalType lda, FadType* x, 
    const OrdinalType incx) const;

      template <typename alpha_type, typename x_type, typename y_type>
      void GER(const OrdinalType m, const OrdinalType n, 
         const alpha_type& alpha, 
         const x_type* x, const OrdinalType incx, 
         const y_type* y, const OrdinalType incy, 
         FadType* A, const OrdinalType lda) const;
      
      

      
      template <typename alpha_type, typename A_type, typename B_type,
    typename beta_type>
      void GEMM(Teuchos::ETransp transa, Teuchos::ETransp transb, 
    const OrdinalType m, const OrdinalType n, const OrdinalType k, 
    const alpha_type& alpha, const A_type* A, const OrdinalType lda,
    const B_type* B, const OrdinalType ldb, const beta_type& beta, 
    FadType* C, const OrdinalType ldc) const;
      
      template <typename alpha_type, typename A_type, typename B_type,
    typename beta_type>
      void SYMM(Teuchos::ESide side, Teuchos::EUplo uplo, const OrdinalType m, 
    const OrdinalType n, 
    const alpha_type& alpha, const A_type* A, 
    const OrdinalType lda, const B_type* B, 
    const OrdinalType ldb,
    const beta_type& beta, FadType* C, 
    const OrdinalType ldc) const;
      
      template <typename alpha_type, typename A_type>
      void TRMM(Teuchos::ESide side, Teuchos::EUplo uplo, 
    Teuchos::ETransp transa, Teuchos::EDiag diag, 
    const OrdinalType m, const OrdinalType n, 
    const alpha_type& alpha, 
    const A_type* A, const OrdinalType lda, 
    FadType* B, const OrdinalType ldb) const;

      template <typename alpha_type, typename A_type>
      void TRSM(Teuchos::ESide side, Teuchos::EUplo uplo, 
    Teuchos::ETransp transa, Teuchos::EDiag diag, 
    const OrdinalType m, const OrdinalType n, 
    const alpha_type& alpha, 
    const A_type* A, const OrdinalType lda, 
    FadType* B, const OrdinalType ldb) const;


    protected:

      ArrayTraits<OrdinalType,FadType> arrayTraits;

      Teuchos::BLAS<OrdinalType, ValueType> blas;

      bool use_default_impl;
      
      mutable std::vector<ValueType> gemv_Ax;

      mutable std::vector<ValueType> gemm_AB;

    protected:

      template <typename x_type, typename y_type>
      void Fad_DOT(const OrdinalType n,
       const x_type* x, 
       const OrdinalType incx, 
       const OrdinalType n_x_dot,  
       const x_type* x_dot, 
       const OrdinalType incx_dot,
       const y_type* y, 
       const OrdinalType incy, 
       const OrdinalType n_y_dot, 
       const y_type* y_dot,
       const OrdinalType incy_dot,
       ValueType& z,
       const OrdinalType n_z_dot,
       ValueType* zdot) const;

      template <typename alpha_type, typename A_type, typename x_type,
    typename beta_type>
      void Fad_GEMV(Teuchos::ETransp trans, 
        const OrdinalType m, 
        const OrdinalType n, 
        const alpha_type& alpha, 
        const OrdinalType n_alpha_dot, 
        const alpha_type* alpha_dot,
        const A_type* A, 
        const OrdinalType lda, 
        const OrdinalType n_A_dot, 
        const A_type* A_dot,
        const OrdinalType lda_dot,
        const x_type* x, 
        const OrdinalType incx, 
        const OrdinalType n_x_dot, 
        const x_type* x_dot, 
        const OrdinalType incx_dot, 
        const beta_type& beta, 
        const OrdinalType n_beta_dot, 
        const beta_type* beta_dot,
        ValueType* y, 
        const OrdinalType incy, 
        const OrdinalType n_y_dot, 
        ValueType* y_dot,
        const OrdinalType incy_dot,
        const OrdinalType n_dot) const;

      template <typename alpha_type, typename x_type, typename y_type>
      void Fad_GER(const OrdinalType m, 
       const OrdinalType n, 
       const alpha_type& alpha, 
       const OrdinalType n_alpha_dot, 
       const alpha_type* alpha_dot,
       const x_type* x, 
       const OrdinalType incx, 
       const OrdinalType n_x_dot, 
       const x_type* x_dot, 
       const OrdinalType incx_dot, 
       const y_type* y, 
       const OrdinalType incy, 
       const OrdinalType n_y_dot, 
       const y_type* y_dot,
       const OrdinalType incy_dot,
       ValueType* A, 
       const OrdinalType lda, 
       const OrdinalType n_A_dot, 
       ValueType* A_dot,
       const OrdinalType lda_dot,
       const OrdinalType n_dot) const;

      template <typename alpha_type, typename A_type, typename B_type,
    typename beta_type>
      void Fad_GEMM(Teuchos::ETransp transa,
        Teuchos::ETransp transb,
        const OrdinalType m, 
        const OrdinalType n, 
        const OrdinalType k,
        const alpha_type& alpha, 
        const OrdinalType n_alpha_dot, 
        const alpha_type* alpha_dot,
        const A_type* A, 
        const OrdinalType lda, 
        const OrdinalType n_A_dot, 
        const A_type* A_dot,
        const OrdinalType lda_dot,
        const B_type* B, 
        const OrdinalType ldb, 
        const OrdinalType n_B_dot, 
        const B_type* B_dot, 
        const OrdinalType ldb_dot, 
        const beta_type& beta, 
        const OrdinalType n_beta_dot, 
        const beta_type* beta_dot,
        ValueType* C, 
        const OrdinalType ldc, 
        const OrdinalType n_C_dot, 
        ValueType* C_dot,
        const OrdinalType ldc_dot,
        const OrdinalType n_dot) const;

      template <typename alpha_type, typename A_type, typename B_type,
    typename beta_type>
      void Fad_SYMM(Teuchos::ESide side, 
        Teuchos::EUplo uplo,
        const OrdinalType m, 
        const OrdinalType n, 
        const alpha_type& alpha, 
        const OrdinalType n_alpha_dot, 
        const alpha_type* alpha_dot,
        const A_type* A, 
        const OrdinalType lda, 
        const OrdinalType n_A_dot, 
        const A_type* A_dot,
        const OrdinalType lda_dot,
        const B_type* B, 
        const OrdinalType ldb, 
        const OrdinalType n_B_dot, 
        const B_type* B_dot, 
        const OrdinalType ldb_dot, 
        const beta_type& beta, 
        const OrdinalType n_beta_dot, 
        const beta_type* beta_dot,
        ValueType* C, 
        const OrdinalType ldc, 
        const OrdinalType n_C_dot, 
        ValueType* C_dot,
        const OrdinalType ldc_dot,
        const OrdinalType n_dot) const;

      template <typename alpha_type, typename A_type>
      void Fad_TRMM(Teuchos::ESide side, 
        Teuchos::EUplo uplo,
        Teuchos::ETransp transa, 
        Teuchos::EDiag diag, 
        const OrdinalType m, 
        const OrdinalType n, 
        const alpha_type& alpha, 
        const OrdinalType n_alpha_dot, 
        const alpha_type* alpha_dot,
        const A_type* A, 
        const OrdinalType lda, 
        const OrdinalType n_A_dot, 
        const A_type* A_dot,
        const OrdinalType lda_dot,
        ValueType* B, 
        const OrdinalType ldb, 
        const OrdinalType n_B_dot, 
        ValueType* B_dot, 
        const OrdinalType ldb_dot, 
        const OrdinalType n_dot) const;

      template <typename alpha_type, typename A_type>
      void Fad_TRSM(Teuchos::ESide side, 
        Teuchos::EUplo uplo,
        Teuchos::ETransp transa, 
        Teuchos::EDiag diag, 
        const OrdinalType m, 
        const OrdinalType n, 
        const alpha_type& alpha, 
        const OrdinalType n_alpha_dot, 
        const alpha_type* alpha_dot,
        const A_type* A, 
        const OrdinalType lda, 
        const OrdinalType n_A_dot, 
        const A_type* A_dot,
        const OrdinalType lda_dot,
        ValueType* B, 
        const OrdinalType ldb, 
        const OrdinalType n_B_dot, 
        ValueType* B_dot, 
        const OrdinalType ldb_dot, 
        const OrdinalType n_dot) const;

    }; // class FadBLAS

  }  // namespace Fad

  // template <typename FadType> ArrayValueType<FadType> { typedef ValueType type; };
  // template <> ArrayValueType<ValueType> { typedef ValueType type; };
  // template <> ArrayValueType<ScalarType> { typedef ScalarType type; };

} // namespace Sacado

// Here we provide partial specializations for Teuchos::BLAS for each Fad type
#define TEUCHOS_BLAS_FAD_SPEC(FADTYPE)          \
namespace Teuchos {             \
  template <typename OrdinalType, typename ValueT>      \
  class BLAS< OrdinalType, FADTYPE<ValueT> > :        \
    public Sacado::Fad::BLAS< OrdinalType, FADTYPE<ValueT> > {    \
  public:               \
    BLAS(bool use_default_impl = true,  bool use_dynamic = true,  \
   OrdinalType static_workspace_size = 0) :     \
      Sacado::Fad::BLAS< OrdinalType, FADTYPE<ValueT> >(    \
  use_default_impl, use_dynamic,static_workspace_size) {}   \
    BLAS(const BLAS& x) :           \
      Sacado::Fad::BLAS< OrdinalType, FADTYPE<ValueT> >(x) {}   \
    virtual ~BLAS() {}              \
  };                  \
}                 \
namespace Sacado {              \
  namespace Fad {             \
    template <typename ValueT>            \
    struct ArrayValueType< FADTYPE<ValueT> > {        \
      typedef ValueT type;            \
    };                  \
  }                 \
}
#define TEUCHOS_BLAS_SFAD_SPEC(FADTYPE)         \
namespace Teuchos {             \
  template <typename OrdinalType, typename ValueT, int Num>   \
  class BLAS< OrdinalType, FADTYPE<ValueT,Num> > :      \
    public Sacado::Fad::BLAS< OrdinalType, FADTYPE<ValueT,Num> > {  \
  public:               \
    BLAS(bool use_default_impl = true,  bool use_dynamic = true,  \
   OrdinalType static_workspace_size = 0) :     \
      Sacado::Fad::BLAS< OrdinalType, FADTYPE<ValueT,Num> >(    \
  use_default_impl, use_dynamic, static_workspace_size) {}  \
    BLAS(const BLAS& x) :           \
      Sacado::Fad::BLAS< OrdinalType, FADTYPE<ValueT,Num> >(x) {} \
    virtual ~BLAS() {}              \
  };                  \
}                 \
namespace Sacado {              \
  namespace Fad {             \
    template <typename ValueT, int Num>         \
    struct ArrayValueType< FADTYPE<ValueT,Num> > {      \
      typedef ValueT type;            \
    };                  \
  }                 \
}
TEUCHOS_BLAS_FAD_SPEC(Sacado::Fad::DFad)
TEUCHOS_BLAS_SFAD_SPEC(Sacado::Fad::SFad)
TEUCHOS_BLAS_SFAD_SPEC(Sacado::Fad::SLFad)
TEUCHOS_BLAS_FAD_SPEC(Sacado::Fad::DMFad)
TEUCHOS_BLAS_FAD_SPEC(Sacado::Fad::DVFad)
TEUCHOS_BLAS_FAD_SPEC(Sacado::ELRFad::DFad)
TEUCHOS_BLAS_SFAD_SPEC(Sacado::ELRFad::SFad)
TEUCHOS_BLAS_SFAD_SPEC(Sacado::ELRFad::SLFad)
TEUCHOS_BLAS_FAD_SPEC(Sacado::CacheFad::DFad)

#undef TEUCHOS_BLAS_FAD_SPEC
#undef TEUCHOS_BLAS_SFAD_SPEC

#include "Sacado_Fad_BLASImp.hpp"

#endif // SACADO_FAD_BLAS_HPP