// Copyright (C) 2026 Kiyotsugu Arai
// SPDX-License-Identifier: LGPL-3.0-or-later

// mkl_backend.hpp
//
// MKL backend implementation
//
// This file provides backend implementations of linear-algebra operations using the
// Intel MKL library. It also includes stub implementations so the code compiles where
// MKL is unavailable.
//
// Main features:
// - C++ wrappers for the MKL versions of BLAS/LAPACK
// - Vector/matrix operation routines
// - Alternative implementations for environments without MKL

#ifndef SANGI_MKL_BACKEND_HPP
#define SANGI_MKL_BACKEND_HPP

#include "math/core/common.hpp"
#include "math/core/traits.hpp"

#include <cmath>
#include <complex>
#include <cstddef>
#include <string>
#include <type_traits>
#include <vector>

// Check MKL availability
#if SANGI_HAS_MKL
#include <mkl.h>
#endif

namespace sangi {
    namespace computation {
        namespace mkl {

            // Add definitions for when MKL is unavailable
#if !SANGI_HAS_MKL
// Mock definitions when CBLAS is unavailable
            typedef enum { CblasRowMajor = 101, CblasColMajor = 102 } CBLAS_LAYOUT;
            typedef enum { CblasNoTrans = 111, CblasTrans = 112, CblasConjTrans = 113 } CBLAS_TRANSPOSE;
            typedef int MKL_INT;
#endif

            // Retrieve the MKL version string
            std::string get_mkl_version();

            // Set the number of MKL threads
            void set_mkl_num_threads(int num_threads);

            // Get the current number of MKL threads
            int get_mkl_num_threads();

            // Set the MKL threading mode
            void set_mkl_threading_mode(bool sequential);

            // Initialize the MKL environment
            bool initialize_mkl();

            // Finalize the MKL environment
            void finalize_mkl();

            //-----------------------------------------------------------------------------
            // BLAS Level 1 functions (vector operations)
            //-----------------------------------------------------------------------------

            // Inner product
            template<typename T>
            T dot(MKL_INT n, const T* x, MKL_INT incx, const T* y, MKL_INT incy) {
#if SANGI_HAS_MKL
                if constexpr (std::is_same_v<T, float>) {
                    return cblas_sdot(n, x, incx, y, incy);
                }
                else if constexpr (std::is_same_v<T, double>) {
                    return cblas_ddot(n, x, incx, y, incy);
                }
                else {
                    // Type not supported by MKL
                    T result = T{};
                    for (MKL_INT i = 0; i < n; ++i) {
                        result += x[i * incx] * y[i * incy];
                    }
                    return result;
                }
#else
                // Standard implementation without MKL
                T result = T{};
                for (MKL_INT i = 0; i < n; ++i) {
                    result += x[i * incx] * y[i * incy];
                }
                return result;
#endif
            }

            // Complex inner product
            template<typename T>
            std::complex<T> dot_complex(MKL_INT n, const std::complex<T>* x, MKL_INT incx,
                const std::complex<T>* y, MKL_INT incy) {
#if SANGI_HAS_MKL
                if constexpr (std::is_same_v<T, float>) {
                    std::complex<float> result;
                    cblas_cdotu_sub(n, x, incx, y, incy, &result);
                    return result;
                }
                else if constexpr (std::is_same_v<T, double>) {
                    std::complex<double> result;
                    cblas_zdotu_sub(n, x, incx, y, incy, &result);
                    return result;
                }
                else {
                    // Type not supported by MKL
                    std::complex<T> result(0, 0);
                    for (MKL_INT i = 0; i < n; ++i) {
                        result += x[i * incx] * y[i * incy];
                    }
                    return result;
                }
#else
                // Standard implementation without MKL
                std::complex<T> result(0, 0);
                for (MKL_INT i = 0; i < n; ++i) {
                    result += x[i * incx] * y[i * incy];
                }
                return result;
#endif
            }

            // Euclidean norm (L2 norm) of a vector
            template<typename T>
            T nrm2(MKL_INT n, const T* x, MKL_INT incx) {
#if SANGI_HAS_MKL
                if constexpr (std::is_same_v<T, float>) {
                    return cblas_snrm2(n, x, incx);
                }
                else if constexpr (std::is_same_v<T, double>) {
                    return cblas_dnrm2(n, x, incx);
                }
                else {
                    // Type not supported by MKL
                    T sum_squares = T{};
                    for (MKL_INT i = 0; i < n; ++i) {
                        sum_squares += x[i * incx] * x[i * incx];
                    }
                    return std::sqrt(sum_squares);
                }
#else
                // Standard implementation without MKL
                T sum_squares = T{};
                for (MKL_INT i = 0; i < n; ++i) {
                    sum_squares += x[i * incx] * x[i * incx];
                }
                return std::sqrt(sum_squares);
#endif
            }

            // Euclidean norm of a complex vector
            template<typename T>
            T nrm2_complex(MKL_INT n, const std::complex<T>* x, MKL_INT incx) {
#if SANGI_HAS_MKL
                if constexpr (std::is_same_v<T, float>) {
                    return cblas_scnrm2(n, x, incx);
                }
                else if constexpr (std::is_same_v<T, double>) {
                    return cblas_dznrm2(n, x, incx);
                }
                else {
                    // Type not supported by MKL
                    T sum_squares = T{};
                    for (MKL_INT i = 0; i < n; ++i) {
                        const auto& val = x[i * incx];
                        sum_squares += val.real() * val.real() + val.imag() * val.imag();
                    }
                    return std::sqrt(sum_squares);
                }
#else
                // Standard implementation without MKL
                T sum_squares = T{};
                for (MKL_INT i = 0; i < n; ++i) {
                    const auto& val = x[i * incx];
                    sum_squares += val.real() * val.real() + val.imag() * val.imag();
                }
                return std::sqrt(sum_squares);
#endif
            }

            // Sum of absolute values of a vector (L1 norm)
            template<typename T>
            T asum(MKL_INT n, const T* x, MKL_INT incx) {
#if SANGI_HAS_MKL
                if constexpr (std::is_same_v<T, float>) {
                    return cblas_sasum(n, x, incx);
                }
                else if constexpr (std::is_same_v<T, double>) {
                    return cblas_dasum(n, x, incx);
                }
                else {
                    // Type not supported by MKL
                    T sum_abs = T{};
                    for (MKL_INT i = 0; i < n; ++i) {
                        sum_abs += std::abs(x[i * incx]);
                    }
                    return sum_abs;
                }
#else
                // Standard implementation without MKL
                T sum_abs = T{};
                for (MKL_INT i = 0; i < n; ++i) {
                    sum_abs += std::abs(x[i * incx]);
                }
                return sum_abs;
#endif
            }

            // Sum of absolute values of a complex vector
            template<typename T>
            T asum_complex(MKL_INT n, const std::complex<T>* x, MKL_INT incx) {
#if SANGI_HAS_MKL
                if constexpr (std::is_same_v<T, float>) {
                    return cblas_scasum(n, x, incx);
                }
                else if constexpr (std::is_same_v<T, double>) {
                    return cblas_dzasum(n, x, incx);
                }
                else {
                    // Type not supported by MKL
                    T sum_abs = T{};
                    for (MKL_INT i = 0; i < n; ++i) {
                        const auto& val = x[i * incx];
                        sum_abs += std::abs(val.real()) + std::abs(val.imag());
                    }
                    return sum_abs;
                }
#else
                // Standard implementation without MKL
                T sum_abs = T{};
                for (MKL_INT i = 0; i < n; ++i) {
                    const auto& val = x[i * incx];
                    sum_abs += std::abs(val.real()) + std::abs(val.imag());
                }
                return sum_abs;
#endif
            }

            // Get the index of the maximum-magnitude element of a vector
            template<typename T>
            MKL_INT iamax(MKL_INT n, const T* x, MKL_INT incx) {
#if SANGI_HAS_MKL
                if constexpr (std::is_same_v<T, float>) {
                    return cblas_isamax(n, x, incx);
                }
                else if constexpr (std::is_same_v<T, double>) {
                    return cblas_idamax(n, x, incx);
                }
                else {
                    // Type not supported by MKL
                    if (n <= 0 || incx <= 0) {
                        return -1;
                    }

                    MKL_INT max_index = 0;
                    T max_abs = std::abs(x[0]);

                    for (MKL_INT i = 1; i < n; ++i) {
                        T abs_val = std::abs(x[i * incx]);
                        if (abs_val > max_abs) {
                            max_abs = abs_val;
                            max_index = i;
                        }
                    }

                    return max_index;
                }
#else
                // Standard implementation without MKL
                if (n <= 0 || incx <= 0) {
                    return -1;
                }

                MKL_INT max_index = 0;
                T max_abs = std::abs(x[0]);

                for (MKL_INT i = 1; i < n; ++i) {
                    T abs_val = std::abs(x[i * incx]);
                    if (abs_val > max_abs) {
                        max_abs = abs_val;
                        max_index = i;
                    }
                }

                return max_index;
#endif
            }

            // Get the index of the maximum-magnitude element of a complex vector
            template<typename T>
            MKL_INT iamax_complex(MKL_INT n, const std::complex<T>* x, MKL_INT incx) {
#if SANGI_HAS_MKL
                if constexpr (std::is_same_v<T, float>) {
                    return cblas_icamax(n, x, incx);
                }
                else if constexpr (std::is_same_v<T, double>) {
                    return cblas_izamax(n, x, incx);
                }
                else {
                    // Type not supported by MKL
                    if (n <= 0 || incx <= 0) {
                        return -1;
                    }

                    MKL_INT max_index = 0;
                    T max_abs = std::abs(x[0].real()) + std::abs(x[0].imag());

                    for (MKL_INT i = 1; i < n; ++i) {
                        const auto& val = x[i * incx];
                        T abs_val = std::abs(val.real()) + std::abs(val.imag());
                        if (abs_val > max_abs) {
                            max_abs = abs_val;
                            max_index = i;
                        }
                    }

                    return max_index;
                }
#else
                // Standard implementation without MKL
                if (n <= 0 || incx <= 0) {
                    return -1;
                }

                MKL_INT max_index = 0;
                T max_abs = std::abs(x[0].real()) + std::abs(x[0].imag());

                for (MKL_INT i = 1; i < n; ++i) {
                    const auto& val = x[i * incx];
                    T abs_val = std::abs(val.real()) + std::abs(val.imag());
                    if (abs_val > max_abs) {
                        max_abs = abs_val;
                        max_index = i;
                    }
                }

                return max_index;
#endif
            }

            // Scalar multiplication (x = alpha*x)
            template<typename T, typename Alpha>
            void scal(MKL_INT n, Alpha alpha, T* x, MKL_INT incx) {
#if SANGI_HAS_MKL
                if constexpr (std::is_same_v<T, float> && std::is_convertible_v<Alpha, float>) {
                    cblas_sscal(n, static_cast<float>(alpha), x, incx);
                }
                else if constexpr (std::is_same_v<T, double> && std::is_convertible_v<Alpha, double>) {
                    cblas_dscal(n, static_cast<double>(alpha), x, incx);
                }
                else {
                    // Type not supported by MKL
                    for (MKL_INT i = 0; i < n; ++i) {
                        x[i * incx] *= static_cast<T>(alpha);
                    }
                }
#else
                // Standard implementation without MKL
                for (MKL_INT i = 0; i < n; ++i) {
                    x[i * incx] *= static_cast<T>(alpha);
                }
#endif
            }

            // Scalar multiplication of a complex vector
            template<typename T, typename Alpha>
            void scal_complex(MKL_INT n, Alpha alpha, std::complex<T>* x, MKL_INT incx) {
#if SANGI_HAS_MKL
                if constexpr (std::is_same_v<T, float> && std::is_convertible_v<Alpha, float>) {
                    cblas_csscal(n, static_cast<float>(alpha), x, incx);
                }
                else if constexpr (std::is_same_v<T, double> && std::is_convertible_v<Alpha, double>) {
                    cblas_zdscal(n, static_cast<double>(alpha), x, incx);
                }
                else {
                    // Type not supported by MKL
                    for (MKL_INT i = 0; i < n; ++i) {
                        x[i * incx] *= static_cast<T>(alpha);
                    }
                }
#else
                // Standard implementation without MKL
                for (MKL_INT i = 0; i < n; ++i) {
                    x[i * incx] *= static_cast<T>(alpha);
                }
#endif
            }

            // Vector copy (y = x)
            template<typename T>
            void copy(MKL_INT n, const T* x, MKL_INT incx, T* y, MKL_INT incy) {
#if SANGI_HAS_MKL
                if constexpr (std::is_same_v<T, float>) {
                    cblas_scopy(n, x, incx, y, incy);
                }
                else if constexpr (std::is_same_v<T, double>) {
                    cblas_dcopy(n, x, incx, y, incy);
                }
                else {
                    // Type not supported by MKL
                    for (MKL_INT i = 0; i < n; ++i) {
                        y[i * incy] = x[i * incx];
                    }
                }
#else
                // Standard implementation without MKL
                for (MKL_INT i = 0; i < n; ++i) {
                    y[i * incy] = x[i * incx];
                }
#endif
            }

            // Complex vector copy
            template<typename T>
            void copy_complex(MKL_INT n, const std::complex<T>* x, MKL_INT incx,
                std::complex<T>* y, MKL_INT incy) {
#if SANGI_HAS_MKL
                if constexpr (std::is_same_v<T, float>) {
                    cblas_ccopy(n, x, incx, y, incy);
                }
                else if constexpr (std::is_same_v<T, double>) {
                    cblas_zcopy(n, x, incx, y, incy);
                }
                else {
                    // Type not supported by MKL
                    for (MKL_INT i = 0; i < n; ++i) {
                        y[i * incy] = x[i * incx];
                    }
                }
#else
                // Standard implementation without MKL
                for (MKL_INT i = 0; i < n; ++i) {
                    y[i * incy] = x[i * incx];
                }
#endif
            }

            // Vector addition (y = alpha*x + y)
            template<typename T, typename Alpha>
            void axpy(MKL_INT n, Alpha alpha, const T* x, MKL_INT incx, T* y, MKL_INT incy) {
#if SANGI_HAS_MKL
                if constexpr (std::is_same_v<T, float> && std::is_convertible_v<Alpha, float>) {
                    cblas_saxpy(n, static_cast<float>(alpha), x, incx, y, incy);
                }
                else if constexpr (std::is_same_v<T, double> && std::is_convertible_v<Alpha, double>) {
                    cblas_daxpy(n, static_cast<double>(alpha), x, incx, y, incy);
                }
                else {
                    // Type not supported by MKL
                    for (MKL_INT i = 0; i < n; ++i) {
                        y[i * incy] += static_cast<T>(alpha) * x[i * incx];
                    }
                }
#else
                // Standard implementation without MKL
                for (MKL_INT i = 0; i < n; ++i) {
                    y[i * incy] += static_cast<T>(alpha) * x[i * incx];
                }
#endif
            }

            // Complex vector addition
            template<typename T>
            void axpy_complex(MKL_INT n, const std::complex<T>& alpha,
                const std::complex<T>* x, MKL_INT incx,
                std::complex<T>* y, MKL_INT incy) {
#if SANGI_HAS_MKL
                if constexpr (std::is_same_v<T, float>) {
                    cblas_caxpy(n, &alpha, x, incx, y, incy);
                }
                else if constexpr (std::is_same_v<T, double>) {
                    cblas_zaxpy(n, &alpha, x, incx, y, incy);
                }
                else {
                    // Type not supported by MKL
                    for (MKL_INT i = 0; i < n; ++i) {
                        y[i * incy] += alpha * x[i * incx];
                    }
                }
#else
                // Standard implementation without MKL
                for (MKL_INT i = 0; i < n; ++i) {
                    y[i * incy] += alpha * x[i * incx];
                }
#endif
            }

            // Vector swap (x <-> y)
            template<typename T>
            void swap(MKL_INT n, T* x, MKL_INT incx, T* y, MKL_INT incy) {
#if SANGI_HAS_MKL
                if constexpr (std::is_same_v<T, float>) {
                    cblas_sswap(n, x, incx, y, incy);
                }
                else if constexpr (std::is_same_v<T, double>) {
                    cblas_dswap(n, x, incx, y, incy);
                }
                else {
                    // Type not supported by MKL
                    for (MKL_INT i = 0; i < n; ++i) {
                        T temp = x[i * incx];
                        x[i * incx] = y[i * incy];
                        y[i * incy] = temp;
                    }
                }
#else
                // Standard implementation without MKL
                for (MKL_INT i = 0; i < n; ++i) {
                    T temp = x[i * incx];
                    x[i * incx] = y[i * incy];
                    y[i * incy] = temp;
                }
#endif
            }

            // Complex vector swap
            template<typename T>
            void swap_complex(MKL_INT n, std::complex<T>* x, MKL_INT incx,
                std::complex<T>* y, MKL_INT incy) {
#if SANGI_HAS_MKL
                if constexpr (std::is_same_v<T, float>) {
                    cblas_cswap(n, x, incx, y, incy);
                }
                else if constexpr (std::is_same_v<T, double>) {
                    cblas_zswap(n, x, incx, y, incy);
                }
                else {
                    // Type not supported by MKL
                    for (MKL_INT i = 0; i < n; ++i) {
                        std::complex<T> temp = x[i * incx];
                        x[i * incx] = y[i * incy];
                        y[i * incy] = temp;
                    }
                }
#else
                // Standard implementation without MKL
                for (MKL_INT i = 0; i < n; ++i) {
                    std::complex<T> temp = x[i * incx];
                    x[i * incx] = y[i * incy];
                    y[i * incy] = temp;
                }
#endif
            }

            //-----------------------------------------------------------------------------
            // BLAS Level 2 functions (matrix-vector operations)
            //-----------------------------------------------------------------------------

            // Matrix-vector product (y = alpha*A*x + beta*y)
            template<typename T>
            inline void gemv(bool transa, std::size_t m, std::size_t n, T alpha,
                const T* a, std::size_t lda, const T* x, std::size_t incx,
                T beta, T* y, std::size_t incy) {
#if SANGI_HAS_MKL
                CBLAS_LAYOUT layout = CblasColMajor;
                CBLAS_TRANSPOSE cblas_transa = transa ? CblasTrans : CblasNoTrans;

                if constexpr (std::is_same_v<T, float>) {
                    cblas_sgemv(layout, cblas_transa, m, n, alpha, a, lda, x, incx, beta, y, incy);
                }
                else if constexpr (std::is_same_v<T, double>) {
                    cblas_dgemv(layout, cblas_transa, m, n, alpha, a, lda, x, incx, beta, y, incy);
                }
                else {
                    // Fall back to the standard implementation
                    for (std::size_t i = 0; i < m; ++i) {
                        T sum = 0;
                        for (std::size_t j = 0; j < n; ++j) {
                            sum += a[i + j * lda] * x[j * incx];
                        }
                        y[i * incy] = alpha * sum + beta * y[i * incy];
                    }
                }
#else
                // Standard implementation
                for (std::size_t i = 0; i < m; ++i) {
                    T sum = 0;
                    for (std::size_t j = 0; j < n; ++j) {
                        if (transa) {
                            sum += a[j + i * lda] * x[j * incx];
                        }
                        else {
                            sum += a[i + j * lda] * x[j * incx];
                        }
                    }
                    y[i * incy] = alpha * sum + beta * y[i * incy];
                }
#endif
            }

            // Complex matrix-vector product
            template<typename T>
            void gemv_complex(bool transa, std::size_t m, std::size_t n,
                const std::complex<T>& alpha, const std::complex<T>* a, std::size_t lda,
                const std::complex<T>* x, std::size_t incx,
                const std::complex<T>& beta, std::complex<T>* y, std::size_t incy) {
#if SANGI_HAS_MKL
                CBLAS_LAYOUT layout = CblasColMajor;
                CBLAS_TRANSPOSE cblas_transa = transa ? CblasTrans : CblasNoTrans;

                if constexpr (std::is_same_v<T, float>) {
                    cblas_cgemv(layout, cblas_transa, m, n, &alpha, a, lda, x, incx, &beta, y, incy);
                }
                else if constexpr (std::is_same_v<T, double>) {
                    cblas_zgemv(layout, cblas_transa, m, n, &alpha, a, lda, x, incx, &beta, y, incy);
                }
                else {
                    // Type not supported by MKL
                    for (std::size_t i = 0; i < m; ++i) {
                        std::complex<T> sum(0, 0);
                        for (std::size_t j = 0; j < n; ++j) {
                            if (transa) {
                                sum += a[j + i * lda] * x[j * incx];
                            }
                            else {
                                sum += a[i + j * lda] * x[j * incx];
                            }
                        }
                        y[i * incy] = alpha * sum + beta * y[i * incy];
                    }
                }
#else
                // Standard implementation without MKL
                for (std::size_t i = 0; i < m; ++i) {
                    std::complex<T> sum(0, 0);
                    for (std::size_t j = 0; j < n; ++j) {
                        if (transa) {
                            sum += a[j + i * lda] * x[j * incx];
                        }
                        else {
                            sum += a[i + j * lda] * x[j * incx];
                        }
                    }
                    y[i * incy] = alpha * sum + beta * y[i * incy];
                }
#endif
            }

            //-----------------------------------------------------------------------------
            // BLAS Level 3 functions (matrix-matrix operations)
            //-----------------------------------------------------------------------------

            // Matrix-matrix product (C = alpha*A*B + beta*C)
            template<typename T>
            inline void gemm(bool transa, bool transb, std::size_t m, std::size_t n, std::size_t k,
                T alpha, const T* a, std::size_t lda, const T* b, std::size_t ldb,
                T beta, T* c, std::size_t ldc) {
#if SANGI_HAS_MKL
                CBLAS_LAYOUT layout = CblasColMajor;
                CBLAS_TRANSPOSE cblas_transa = transa ? CblasTrans : CblasNoTrans;
                CBLAS_TRANSPOSE cblas_transb = transb ? CblasTrans : CblasNoTrans;

                if constexpr (std::is_same_v<T, float>) {
                    cblas_sgemm(layout, cblas_transa, cblas_transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc);
                }
                else if constexpr (std::is_same_v<T, double>) {
                    cblas_dgemm(layout, cblas_transa, cblas_transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc);
                }
                else {
                    // Fall back to the standard implementation
                    for (std::size_t i = 0; i < m; ++i) {
                        for (std::size_t j = 0; j < n; ++j) {
                            T sum = 0;
                            for (std::size_t l = 0; l < k; ++l) {
                                sum += a[i + l * lda] * b[l + j * ldb];
                            }
                            c[i + j * ldc] = alpha * sum + beta * c[i + j * ldc];
                        }
                    }
                }
#else
                // Standard implementation
                for (std::size_t i = 0; i < m; ++i) {
                    for (std::size_t j = 0; j < n; ++j) {
                        T sum = 0;
                        for (std::size_t l = 0; l < k; ++l) {
                            if (transa && transb) {
                                sum += a[l + i * lda] * b[j + l * ldb];
                            }
                            else if (transa) {
                                sum += a[l + i * lda] * b[l + j * ldb];
                            }
                            else if (transb) {
                                sum += a[i + l * lda] * b[j + l * ldb];
                            }
                            else {
                                sum += a[i + l * lda] * b[l + j * ldb];
                            }
                        }
                        c[i + j * ldc] = alpha * sum + beta * c[i + j * ldc];
                    }
                }
#endif
            }

            // Matrix-matrix product (row-major layout): C = alpha*op(A)*op(B) + beta*C
            template<typename T>
            inline void gemm_row_major(bool transa, bool transb,
                std::size_t m, std::size_t n, std::size_t k,
                T alpha, const T* a, std::size_t lda, const T* b, std::size_t ldb,
                T beta, T* c, std::size_t ldc) {
#if SANGI_HAS_MKL
                CBLAS_TRANSPOSE cblas_transa = transa ? CblasTrans : CblasNoTrans;
                CBLAS_TRANSPOSE cblas_transb = transb ? CblasTrans : CblasNoTrans;

                if constexpr (std::is_same_v<T, float>) {
                    cblas_sgemm(CblasRowMajor, cblas_transa, cblas_transb,
                                m, n, k, alpha, a, lda, b, ldb, beta, c, ldc);
                }
                else if constexpr (std::is_same_v<T, double>) {
                    cblas_dgemm(CblasRowMajor, cblas_transa, cblas_transb,
                                m, n, k, alpha, a, lda, b, ldb, beta, c, ldc);
                }
                else {
                    // Fall back to the standard implementation (row-major)
                    for (std::size_t i = 0; i < m; ++i) {
                        for (std::size_t j = 0; j < n; ++j) {
                            T sum = 0;
                            for (std::size_t l = 0; l < k; ++l) {
                                T a_val = transa ? a[l * lda + i] : a[i * lda + l];
                                T b_val = transb ? b[j * ldb + l] : b[l * ldb + j];
                                sum += a_val * b_val;
                            }
                            c[i * ldc + j] = alpha * sum + beta * c[i * ldc + j];
                        }
                    }
                }
#else
                // Standard implementation (row-major)
                for (std::size_t i = 0; i < m; ++i) {
                    for (std::size_t j = 0; j < n; ++j) {
                        T sum = 0;
                        for (std::size_t l = 0; l < k; ++l) {
                            T a_val = transa ? a[l * lda + i] : a[i * lda + l];
                            T b_val = transb ? b[j * ldb + l] : b[l * ldb + j];
                            sum += a_val * b_val;
                        }
                        c[i * ldc + j] = alpha * sum + beta * c[i * ldc + j];
                    }
                }
#endif
            }

            // Complex matrix-matrix product
            template<typename T>
            void gemm_complex(bool transa, bool transb, std::size_t m, std::size_t n, std::size_t k,
                const std::complex<T>& alpha, const std::complex<T>* a, std::size_t lda,
                const std::complex<T>* b, std::size_t ldb,
                const std::complex<T>& beta, std::complex<T>* c, std::size_t ldc) {
#if SANGI_HAS_MKL
                CBLAS_LAYOUT layout = CblasColMajor;
                CBLAS_TRANSPOSE cblas_transa = transa ? CblasTrans : CblasNoTrans;
                CBLAS_TRANSPOSE cblas_transb = transb ? CblasTrans : CblasNoTrans;

                if constexpr (std::is_same_v<T, float>) {
                    cblas_cgemm(layout, cblas_transa, cblas_transb, m, n, k, &alpha, a, lda, b, ldb, &beta, c, ldc);
                }
                else if constexpr (std::is_same_v<T, double>) {
                    cblas_zgemm(layout, cblas_transa, cblas_transb, m, n, k, &alpha, a, lda, b, ldb, &beta, c, ldc);
                }
                else {
                    // Type not supported by MKL
                    for (std::size_t i = 0; i < m; ++i) {
                        for (std::size_t j = 0; j < n; ++j) {
                            std::complex<T> sum(0, 0);
                            for (std::size_t l = 0; l < k; ++l) {
                                if (transa && transb) {
                                    sum += a[l + i * lda] * b[j + l * ldb];
                                }
                                else if (transa) {
                                    sum += a[l + i * lda] * b[l + j * ldb];
                                }
                                else if (transb) {
                                    sum += a[i + l * lda] * b[j + l * ldb];
                                }
                                else {
                                    sum += a[i + l * lda] * b[l + j * ldb];
                                }
                            }
                            c[i + j * ldc] = alpha * sum + beta * c[i + j * ldc];
                        }
                    }
                }
#else
                // Standard implementation without MKL
                for (std::size_t i = 0; i < m; ++i) {
                    for (std::size_t j = 0; j < n; ++j) {
                        std::complex<T> sum(0, 0);
                        for (std::size_t l = 0; l < k; ++l) {
                            if (transa && transb) {
                                sum += a[l + i * lda] * b[j + l * ldb];
                            }
                            else if (transa) {
                                sum += a[l + i * lda] * b[l + j * ldb];
                            }
                            else if (transb) {
                                sum += a[i + l * lda] * b[j + l * ldb];
                            }
                            else {
                                sum += a[i + l * lda] * b[l + j * ldb];
                            }
                        }
                        c[i + j * ldc] = alpha * sum + beta * c[i + j * ldc];
                    }
                }
#endif
            }

            //-----------------------------------------------------------------------------
            // Helper functions (used by the computation policy)
            //-----------------------------------------------------------------------------

            // SIMD-friendly inner product
            template<typename T>
            T dot_product_simd(const T* x, const T* y, std::size_t size) {
#if SANGI_HAS_MKL
                return dot(static_cast<MKL_INT>(size), x, 1, y, 1);
#else
                // Fall back to the SIMD backend or the standard implementation
                T result = T{};
                for (std::size_t i = 0; i < size; ++i) {
                    result += x[i] * y[i];
                }
                return result;
#endif
            }

            // SIMD-friendly vector-norm computation
            template<typename T>
            T norm2_simd(const T* x, std::size_t size) {
#if SANGI_HAS_MKL
                return nrm2(static_cast<MKL_INT>(size), x, 1);
#else
                // Fall back to the SIMD backend or the standard implementation
                T sum_squares = T{};
                for (std::size_t i = 0; i < size; ++i) {
                    sum_squares += x[i] * x[i];
                }
                return std::sqrt(sum_squares);
#endif
            }

            // SIMD-friendly axpy
            template<typename T, typename Alpha>
            void axpy_simd(T* y, Alpha alpha, const T* x, std::size_t size) {
#if SANGI_HAS_MKL
                axpy(static_cast<MKL_INT>(size), alpha, x, 1, y, 1);
#else
                // Fall back to the SIMD backend or the standard implementation
                for (std::size_t i = 0; i < size; ++i) {
                    y[i] += static_cast<T>(alpha) * x[i];
                }
#endif
            }

            // SIMD-friendly scaling
            template<typename T, typename Alpha>
            void scale_simd(T* x, Alpha alpha, std::size_t size) {
#if SANGI_HAS_MKL
                scal(static_cast<MKL_INT>(size), alpha, x, 1);
#else
                // Fall back to the SIMD backend or the standard implementation
                for (std::size_t i = 0; i < size; ++i) {
                    x[i] *= static_cast<T>(alpha);
                }
#endif
            }

            // SIMD-friendly vector addition
            template<typename T>
            void add_simd(const T* x, const T* y, T* z, std::size_t size) {
#if SANGI_HAS_MKL
                // Use axpy under MKL
                copy(static_cast<MKL_INT>(size), x, 1, z, 1);
                axpy(static_cast<MKL_INT>(size), static_cast<T>(1), y, 1, z, 1);
#else
                // Fall back to the SIMD backend or the standard implementation
                for (std::size_t i = 0; i < size; ++i) {
                    z[i] = x[i] + y[i];
                }
#endif
            }

            // SIMD-friendly vector subtraction
            template<typename T>
            void subtract_simd(const T* x, const T* y, T* z, std::size_t size) {
#if SANGI_HAS_MKL
                // Use axpy under MKL
                copy(static_cast<MKL_INT>(size), x, 1, z, 1);
                axpy(static_cast<MKL_INT>(size), static_cast<T>(-1), y, 1, z, 1);
#else
                // Fall back to the SIMD backend or the standard implementation
                for (std::size_t i = 0; i < size; ++i) {
                    z[i] = x[i] - y[i];
                }
#endif
            }

            // SIMD-friendly element-wise multiplication
            template<typename T>
            void multiply_elements_simd(const T* x, const T* y, T* z, std::size_t size) {
                // Not in MKL; use the SIMD or standard implementation
                for (std::size_t i = 0; i < size; ++i) {
                    z[i] = x[i] * y[i];
                }
            }

            // SIMD-friendly element-wise division
            template<typename T>
            void divide_elements_simd(const T* x, const T* y, T* z, std::size_t size) {
                // Not in MKL; use the SIMD or standard implementation
                for (std::size_t i = 0; i < size; ++i) {
                    z[i] = x[i] / y[i];
                }
            }

        } // namespace mkl
    } // namespace computation
} // namespace sangi

#endif // SANGI_MKL_BACKEND_HPP