// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.

#include <mkl_dfti.h>

#include "./InternalHeaderCheck.h"

#include <complex>

namespace Eigen {
namespace internal {
namespace imklfft {

#define RUN_OR_ASSERT(EXPR, ERROR_MSG)                   \
  {                                                      \
    MKL_LONG status = (EXPR);                            \
    eigen_assert(status == DFTI_NO_ERROR && (ERROR_MSG)); \
  };

inline MKL_Complex16* complex_cast(const std::complex<double>* p) {
  return const_cast<MKL_Complex16*>(reinterpret_cast<const MKL_Complex16*>(p));
}

inline MKL_Complex8* complex_cast(const std::complex<float>* p) {
  return const_cast<MKL_Complex8*>(reinterpret_cast<const MKL_Complex8*>(p));
}

/*
 * Parameters:
 * precision: enum, Precision of the transform: DFTI_SINGLE or DFTI_DOUBLE.
 * forward_domain: enum, Forward domain of the transform: DFTI_COMPLEX or
 * DFTI_REAL. dimension: MKL_LONG Dimension of the transform. sizes: MKL_LONG if
 * dimension = 1.Length of the transform for a one-dimensional transform. sizes:
 * Array of type MKL_LONG otherwise. Lengths of each dimension for a
 * multi-dimensional transform.
 */
inline void configure_descriptor(DFTI_DESCRIPTOR_HANDLE* handl,
                                 enum DFTI_CONFIG_VALUE precision,
                                 enum DFTI_CONFIG_VALUE forward_domain,
                                 MKL_LONG dimension, MKL_LONG* sizes) {
  eigen_assert(dimension == 1 ||
               dimension == 2 &&
                   "Transformation dimension must be less than 3.");

  if (dimension == 1) {
    RUN_OR_ASSERT(DftiCreateDescriptor(handl, precision, forward_domain,
                                       dimension, *sizes),
                  "DftiCreateDescriptor failed.")
    if (forward_domain == DFTI_REAL) {
      // Set CCE storage
      RUN_OR_ASSERT(DftiSetValue(*handl, DFTI_CONJUGATE_EVEN_STORAGE,
                                 DFTI_COMPLEX_COMPLEX),
                    "DftiSetValue failed.")
    }
  } else {
    RUN_OR_ASSERT(
        DftiCreateDescriptor(handl, precision, DFTI_COMPLEX, dimension, sizes),
        "DftiCreateDescriptor failed.")
  }

  RUN_OR_ASSERT(DftiSetValue(*handl, DFTI_PLACEMENT, DFTI_NOT_INPLACE),
                "DftiSetValue failed.")
  RUN_OR_ASSERT(DftiCommitDescriptor(*handl), "DftiCommitDescriptor failed.")
}

template <typename T>
struct plan {};

template <>
struct plan<float> {
  typedef float scalar_type;
  typedef MKL_Complex8 complex_type;

  DFTI_DESCRIPTOR_HANDLE m_plan;

  plan() : m_plan(0) {}
  ~plan() {
    if (m_plan) DftiFreeDescriptor(&m_plan);
  };

  enum DFTI_CONFIG_VALUE precision = DFTI_SINGLE;

  inline void forward(complex_type* dst, complex_type* src, MKL_LONG nfft) {
    if (m_plan == 0) {
      configure_descriptor(&m_plan, precision, DFTI_COMPLEX, 1, &nfft);
    }
    RUN_OR_ASSERT(DftiComputeForward(m_plan, src, dst),
                  "DftiComputeForward failed.")
  }

  inline void inverse(complex_type* dst, complex_type* src, MKL_LONG nfft) {
    if (m_plan == 0) {
      configure_descriptor(&m_plan, precision, DFTI_COMPLEX, 1, &nfft);
    }
    RUN_OR_ASSERT(DftiComputeBackward(m_plan, src, dst),
                  "DftiComputeBackward failed.")
  }

  inline void forward(complex_type* dst, scalar_type* src, MKL_LONG nfft) {
    if (m_plan == 0) {
      configure_descriptor(&m_plan, precision, DFTI_REAL, 1, &nfft);
    }
    RUN_OR_ASSERT(DftiComputeForward(m_plan, src, dst),
                  "DftiComputeForward failed.")
  }

  inline void inverse(scalar_type* dst, complex_type* src, MKL_LONG nfft) {
    if (m_plan == 0) {
      configure_descriptor(&m_plan, precision, DFTI_REAL, 1, &nfft);
    }
    RUN_OR_ASSERT(DftiComputeBackward(m_plan, src, dst),
                  "DftiComputeBackward failed.")
  }

  inline void forward2(complex_type* dst, complex_type* src, int n0, int n1) {
    if (m_plan == 0) {
      MKL_LONG sizes[2] = {n0, n1};
      configure_descriptor(&m_plan, precision, DFTI_COMPLEX, 2, sizes);
    }
    RUN_OR_ASSERT(DftiComputeForward(m_plan, src, dst),
                  "DftiComputeForward failed.")
  }

  inline void inverse2(complex_type* dst, complex_type* src, int n0, int n1) {
    if (m_plan == 0) {
      MKL_LONG sizes[2] = {n0, n1};
      configure_descriptor(&m_plan, precision, DFTI_COMPLEX, 2, sizes);
    }
    RUN_OR_ASSERT(DftiComputeBackward(m_plan, src, dst),
                  "DftiComputeBackward failed.")
  }
};

template <>
struct plan<double> {
  typedef double scalar_type;
  typedef MKL_Complex16 complex_type;

  DFTI_DESCRIPTOR_HANDLE m_plan;

  plan() : m_plan(0) {}
  ~plan() {
    if (m_plan) DftiFreeDescriptor(&m_plan);
  };

  enum DFTI_CONFIG_VALUE precision = DFTI_DOUBLE;

  inline void forward(complex_type* dst, complex_type* src, MKL_LONG nfft) {
    if (m_plan == 0) {
      configure_descriptor(&m_plan, precision, DFTI_COMPLEX, 1, &nfft);
    }
    RUN_OR_ASSERT(DftiComputeForward(m_plan, src, dst),
                  "DftiComputeForward failed.")
  }

  inline void inverse(complex_type* dst, complex_type* src, MKL_LONG nfft) {
    if (m_plan == 0) {
      configure_descriptor(&m_plan, precision, DFTI_COMPLEX, 1, &nfft);
    }
    RUN_OR_ASSERT(DftiComputeBackward(m_plan, src, dst),
                  "DftiComputeBackward failed.")
  }

  inline void forward(complex_type* dst, scalar_type* src, MKL_LONG nfft) {
    if (m_plan == 0) {
      configure_descriptor(&m_plan, precision, DFTI_REAL, 1, &nfft);
    }
    RUN_OR_ASSERT(DftiComputeForward(m_plan, src, dst),
                  "DftiComputeForward failed.")
  }

  inline void inverse(scalar_type* dst, complex_type* src, MKL_LONG nfft) {
    if (m_plan == 0) {
      configure_descriptor(&m_plan, precision, DFTI_REAL, 1, &nfft);
    }
    RUN_OR_ASSERT(DftiComputeBackward(m_plan, src, dst),
                  "DftiComputeBackward failed.")
  }

  inline void forward2(complex_type* dst, complex_type* src, int n0, int n1) {
    if (m_plan == 0) {
      MKL_LONG sizes[2] = {n0, n1};
      configure_descriptor(&m_plan, precision, DFTI_COMPLEX, 2, sizes);
    }
    RUN_OR_ASSERT(DftiComputeForward(m_plan, src, dst),
                  "DftiComputeForward failed.")
  }

  inline void inverse2(complex_type* dst, complex_type* src, int n0, int n1) {
    if (m_plan == 0) {
      MKL_LONG sizes[2] = {n0, n1};
      configure_descriptor(&m_plan, precision, DFTI_COMPLEX, 2, sizes);
    }
    RUN_OR_ASSERT(DftiComputeBackward(m_plan, src, dst),
                  "DftiComputeBackward failed.")
  }
};

template <typename Scalar_>
struct imklfft_impl {
  typedef Scalar_ Scalar;
  typedef std::complex<Scalar> Complex;

  inline void clear() { m_plans.clear(); }

  // complex-to-complex forward FFT
  inline void fwd(Complex* dst, const Complex* src, int nfft) {
    MKL_LONG size = nfft;
    get_plan(nfft, dst, src)
        .forward(complex_cast(dst), complex_cast(src), size);
  }

  // real-to-complex forward FFT
  inline void fwd(Complex* dst, const Scalar* src, int nfft) {
    MKL_LONG size = nfft;
    get_plan(nfft, dst, src)
        .forward(complex_cast(dst), const_cast<Scalar*>(src), nfft);
  }

  // 2-d complex-to-complex
  inline void fwd2(Complex* dst, const Complex* src, int n0, int n1) {
    get_plan(n0, n1, dst, src)
        .forward2(complex_cast(dst), complex_cast(src), n0, n1);
  }

  // inverse complex-to-complex
  inline void inv(Complex* dst, const Complex* src, int nfft) {
    MKL_LONG size = nfft;
    get_plan(nfft, dst, src)
        .inverse(complex_cast(dst), complex_cast(src), nfft);
  }

  // half-complex to scalar
  inline void inv(Scalar* dst, const Complex* src, int nfft) {
    MKL_LONG size = nfft;
    get_plan(nfft, dst, src)
        .inverse(const_cast<Scalar*>(dst), complex_cast(src), nfft);
  }

  // 2-d complex-to-complex
  inline void inv2(Complex* dst, const Complex* src, int n0, int n1) {
    get_plan(n0, n1, dst, src)
        .inverse2(complex_cast(dst), complex_cast(src), n0, n1);
  }

 private:
  std::map<int64_t, plan<Scalar>> m_plans;

  inline plan<Scalar>& get_plan(int nfft, void* dst,
                                const void* src) {
    int inplace = dst == src ? 1 : 0;
    int aligned = ((reinterpret_cast<size_t>(src) & 15) |
                   (reinterpret_cast<size_t>(dst) & 15)) == 0
                      ? 1
                      : 0;
    int64_t key = ((nfft << 2) | (inplace << 1) | aligned)
                  << 1;

    // Create element if key does not exist.
    return m_plans[key];
  }

  inline plan<Scalar>& get_plan(int n0, int n1, void* dst,
                                const void* src) {
    int inplace = (dst == src) ? 1 : 0;
    int aligned = ((reinterpret_cast<size_t>(src) & 15) |
                   (reinterpret_cast<size_t>(dst) & 15)) == 0
                      ? 1
                      : 0;
    int64_t key = (((((int64_t)n0) << 31) | (n1 << 2) |
                    (inplace << 1) | aligned)
                   << 1) +
                  1;

    // Create element if key does not exist.
    return m_plans[key];
  }
};

#undef RUN_OR_ASSERT

}  // namespace imklfft
}  // namespace internal
}  // namespace Eigen