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

// mkl_backend.hpp
//
// MKLバックエンド実装
//
// このファイルでは、Intel MKLライブラリを使用した線形代数演算のバックエンド実装を
// 提供します。MKLが利用できない環境ではコンパイルエラーを回避するための実装も含まれています。
//
// 主な機能:
// - MKL版BLAS/LAPACKへのC++ラッパー
// - ベクトル/行列演算関数
// - MKLが利用不可の場合の代替実装

#ifndef CALX_MKL_BACKEND_HPP
#define CALX_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>

// MKL利用可能性の確認
#if CALX_HAS_MKL
#include <mkl.h>
#endif

namespace calx {
    namespace computation {
        namespace mkl {

            // MKLが利用できない場合の定義を追加
#if !CALX_HAS_MKL
// CBLAS定義が存在しない場合のモックアップ
            typedef enum { CblasRowMajor = 101, CblasColMajor = 102 } CBLAS_LAYOUT;
            typedef enum { CblasNoTrans = 111, CblasTrans = 112, CblasConjTrans = 113 } CBLAS_TRANSPOSE;
            typedef int MKL_INT;
#endif

            // MKLバージョン情報の取得
            std::string get_mkl_version();

            // MKLスレッド数の設定
            void set_mkl_num_threads(int num_threads);

            // 現在のMKLスレッド数の取得
            int get_mkl_num_threads();

            // MKLのスレッドモードの設定
            void set_mkl_threading_mode(bool sequential);

            // MKL環境の初期化
            bool initialize_mkl();

            // MKL環境の終了処理
            void finalize_mkl();

            //-----------------------------------------------------------------------------
            // BLAS Level 1関数 (ベクトル演算)
            //-----------------------------------------------------------------------------

            // 内積計算
            template<typename T>
            T dot(MKL_INT n, const T* x, MKL_INT incx, const T* y, MKL_INT incy) {
#if CALX_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 {
                    // MKLがサポートしない型
                    T result = T{};
                    for (MKL_INT i = 0; i < n; ++i) {
                        result += x[i * incx] * y[i * incy];
                    }
                    return result;
                }
#else
                // MKLなしの標準的な実装
                T result = T{};
                for (MKL_INT i = 0; i < n; ++i) {
                    result += x[i * incx] * y[i * incy];
                }
                return result;
#endif
            }

            // 複素数内積計算
            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 CALX_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 {
                    // 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
                // 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
            }

            // ベクトルのユークリッドノルム（L2ノルム）
            template<typename T>
            T nrm2(MKL_INT n, const T* x, MKL_INT incx) {
#if CALX_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 {
                    // 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
                // 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
            }

            // 複素数ベクトルのユークリッドノルム
            template<typename T>
            T nrm2_complex(MKL_INT n, const std::complex<T>* x, MKL_INT incx) {
#if CALX_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 {
                    // 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
                // 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
            }

            // ベクトルの絶対値の総和（L1ノルム）
            template<typename T>
            T asum(MKL_INT n, const T* x, MKL_INT incx) {
#if CALX_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 {
                    // MKLがサポートしない型
                    T sum_abs = T{};
                    for (MKL_INT i = 0; i < n; ++i) {
                        sum_abs += std::abs(x[i * incx]);
                    }
                    return sum_abs;
                }
#else
                // MKLなしの標準的な実装
                T sum_abs = T{};
                for (MKL_INT i = 0; i < n; ++i) {
                    sum_abs += std::abs(x[i * incx]);
                }
                return sum_abs;
#endif
            }

            // 複素数ベクトルの絶対値の総和
            template<typename T>
            T asum_complex(MKL_INT n, const std::complex<T>* x, MKL_INT incx) {
#if CALX_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 {
                    // 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
                // 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
            }

            // ベクトルの最大値要素のインデックスを取得
            template<typename T>
            MKL_INT iamax(MKL_INT n, const T* x, MKL_INT incx) {
#if CALX_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 {
                    // 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
                // 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
            }

            // 複素数ベクトルの最大値要素のインデックスを取得
            template<typename T>
            MKL_INT iamax_complex(MKL_INT n, const std::complex<T>* x, MKL_INT incx) {
#if CALX_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 {
                    // 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
                // 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
            }

            // スカラー倍（x = alpha*x）
            template<typename T, typename Alpha>
            void scal(MKL_INT n, Alpha alpha, T* x, MKL_INT incx) {
#if CALX_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 {
                    // MKLがサポートしない型
                    for (MKL_INT i = 0; i < n; ++i) {
                        x[i * incx] *= static_cast<T>(alpha);
                    }
                }
#else
                // MKLなしの標準的な実装
                for (MKL_INT i = 0; i < n; ++i) {
                    x[i * incx] *= static_cast<T>(alpha);
                }
#endif
            }

            // 複素数ベクトルのスカラー倍
            template<typename T, typename Alpha>
            void scal_complex(MKL_INT n, Alpha alpha, std::complex<T>* x, MKL_INT incx) {
#if CALX_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 {
                    // MKLがサポートしない型
                    for (MKL_INT i = 0; i < n; ++i) {
                        x[i * incx] *= static_cast<T>(alpha);
                    }
                }
#else
                // MKLなしの標準的な実装
                for (MKL_INT i = 0; i < n; ++i) {
                    x[i * incx] *= static_cast<T>(alpha);
                }
#endif
            }

            // ベクトルのコピー（y = x）
            template<typename T>
            void copy(MKL_INT n, const T* x, MKL_INT incx, T* y, MKL_INT incy) {
#if CALX_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 {
                    // MKLがサポートしない型
                    for (MKL_INT i = 0; i < n; ++i) {
                        y[i * incy] = x[i * incx];
                    }
                }
#else
                // MKLなしの標準的な実装
                for (MKL_INT i = 0; i < n; ++i) {
                    y[i * incy] = x[i * incx];
                }
#endif
            }

            // 複素数ベクトルのコピー
            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 CALX_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 {
                    // MKLがサポートしない型
                    for (MKL_INT i = 0; i < n; ++i) {
                        y[i * incy] = x[i * incx];
                    }
                }
#else
                // MKLなしの標準的な実装
                for (MKL_INT i = 0; i < n; ++i) {
                    y[i * incy] = x[i * incx];
                }
#endif
            }

            // ベクトルの加算（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 CALX_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 {
                    // MKLがサポートしない型
                    for (MKL_INT i = 0; i < n; ++i) {
                        y[i * incy] += static_cast<T>(alpha) * x[i * incx];
                    }
                }
#else
                // MKLなしの標準的な実装
                for (MKL_INT i = 0; i < n; ++i) {
                    y[i * incy] += static_cast<T>(alpha) * x[i * incx];
                }
#endif
            }

            // 複素数ベクトルの加算
            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 CALX_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 {
                    // MKLがサポートしない型
                    for (MKL_INT i = 0; i < n; ++i) {
                        y[i * incy] += alpha * x[i * incx];
                    }
                }
#else
                // MKLなしの標準的な実装
                for (MKL_INT i = 0; i < n; ++i) {
                    y[i * incy] += alpha * x[i * incx];
                }
#endif
            }

            // ベクトルの交換（x <-> y）
            template<typename T>
            void swap(MKL_INT n, T* x, MKL_INT incx, T* y, MKL_INT incy) {
#if CALX_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 {
                    // 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
                // 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
            }

            // 複素数ベクトルの交換
            template<typename T>
            void swap_complex(MKL_INT n, std::complex<T>* x, MKL_INT incx,
                std::complex<T>* y, MKL_INT incy) {
#if CALX_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 {
                    // 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
                // 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関数 (行列-ベクトル演算)
            //-----------------------------------------------------------------------------

            // 行列-ベクトル積（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 CALX_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 {
                    // 標準実装にフォールバック
                    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
                // 標準実装
                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
            }

            // 複素数行列-ベクトル積
            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 CALX_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 {
                    // 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
                // 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関数 (行列-行列演算)
            //-----------------------------------------------------------------------------

            // 行列-行列積（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 CALX_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 {
                    // 標準実装にフォールバック
                    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
                // 標準実装
                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
            }

            // 複素数行列-行列積
            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 CALX_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 {
                    // 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
                // 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
            }

            //-----------------------------------------------------------------------------
            // ヘルパー関数（計算ポリシーから使用）
            //-----------------------------------------------------------------------------

            // SIMD対応の内積計算
            template<typename T>
            T dot_product_simd(const T* x, const T* y, std::size_t size) {
#if CALX_HAS_MKL
                return dot(static_cast<MKL_INT>(size), x, 1, y, 1);
#else
                // SIMDバックエンドまたは標準実装にフォールバック
                T result = T{};
                for (std::size_t i = 0; i < size; ++i) {
                    result += x[i] * y[i];
                }
                return result;
#endif
            }

            // SIMD対応のベクトルノルム計算
            template<typename T>
            T norm2_simd(const T* x, std::size_t size) {
#if CALX_HAS_MKL
                return nrm2(static_cast<MKL_INT>(size), x, 1);
#else
                // SIMDバックエンドまたは標準実装にフォールバック
                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対応のaxpy演算
            template<typename T, typename Alpha>
            void axpy_simd(T* y, Alpha alpha, const T* x, std::size_t size) {
#if CALX_HAS_MKL
                axpy(static_cast<MKL_INT>(size), alpha, x, 1, y, 1);
#else
                // SIMDバックエンドまたは標準実装にフォールバック
                for (std::size_t i = 0; i < size; ++i) {
                    y[i] += static_cast<T>(alpha) * x[i];
                }
#endif
            }

            // SIMD対応のスケーリング演算
            template<typename T, typename Alpha>
            void scale_simd(T* x, Alpha alpha, std::size_t size) {
#if CALX_HAS_MKL
                scal(static_cast<MKL_INT>(size), alpha, x, 1);
#else
                // SIMDバックエンドまたは標準実装にフォールバック
                for (std::size_t i = 0; i < size; ++i) {
                    x[i] *= static_cast<T>(alpha);
                }
#endif
            }

            // SIMD対応のベクトル加算
            template<typename T>
            void add_simd(const T* x, const T* y, T* z, std::size_t size) {
#if CALX_HAS_MKL
                // MKLではaxpyを使用
                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
                // SIMDバックエンドまたは標準実装にフォールバック
                for (std::size_t i = 0; i < size; ++i) {
                    z[i] = x[i] + y[i];
                }
#endif
            }

            // SIMD対応のベクトル減算
            template<typename T>
            void subtract_simd(const T* x, const T* y, T* z, std::size_t size) {
#if CALX_HAS_MKL
                // MKLではaxpyを使用
                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
                // SIMDバックエンドまたは標準実装にフォールバック
                for (std::size_t i = 0; i < size; ++i) {
                    z[i] = x[i] - y[i];
                }
#endif
            }

            // SIMD対応の要素ごとの乗算
            template<typename T>
            void multiply_elements_simd(const T* x, const T* y, T* z, std::size_t size) {
                // MKLにはないので、SIMDまたは標準実装を使用
                for (std::size_t i = 0; i < size; ++i) {
                    z[i] = x[i] * y[i];
                }
            }

            // SIMD対応の要素ごとの除算
            template<typename T>
            void divide_elements_simd(const T* x, const T* y, T* z, std::size_t size) {
                // MKLにはないので、SIMDまたは標準実装を使用
                for (std::size_t i = 0; i < size; ++i) {
                    z[i] = x[i] / y[i];
                }
            }

        } // namespace mkl
    } // namespace computation
} // namespace calx

#endif // CALX_MKL_BACKEND_HPP