[线性代数] OpenBLAS在c++、Fortran中效率差距这么大？

module inv
contains
    function inverse(A) result(Ainv)
        real(kind=8), dimension(:,:), intent(in) :: A
        real(kind=8), dimension(size(A,1),size(A,2)) :: Ainv
        real(kind=8), dimension(size(A,1)) :: work
        integer, dimension(size(A,1)) :: ipiv 
        integer :: n, info

        print *, "Step 1:: 矩阵求逆中 ... ..."
        !external DGETRF
        !external DGETRI
        Ainv = A
        n = size(A,1)
        call DGETRF(n, n, Ainv, n, ipiv, info)
        if (info /= 0) then
           stop '矩阵是 奇异阵 !'
        end if
        call DGETRI(n, Ainv, n, ipiv, work, n, info)
        if (info /= 0) then
           stop '矩阵求逆 失败 !'
        end if
        print *, "Step 2:: 矩阵相乘中... ..."
    end function inverse   
end module


program MatrixSpeed

    use inv
    implicit none
    real(kind=8), allocatable :: z(:,:), re(:,:)
    integer             :: i, j, n, p
    real(kind=8)        :: u
    real(kind=4)        :: t1, t2

    n = 1500
    p = 1500
    allocate(z(n,p),re(n,n))
    do i=1,n
        do j=1,p
            call random_number(u)
            z(i,j) = u * 3.0
        enddo
    enddo

    ! 随机矩阵求逆，再乘以原矩阵，得到单位阵。
    print *, "N= ", n
    print *, "正在计算... ..."
    call cpu_time(t1)
    re = matmul(z, inverse(z))
    call cpu_time(t2)
    print *, "Time:  ", t2-t1, " 秒"
    do i=1,min(10,n)
        write (*,'(10(f10.4, x))') re(i, 1:min(10,n))
    enddo

end program MatrixSpeed

#include <iostream>
#include <vector>
#include <random>
#include <ctime>
#include <iomanip>
#include <cblas.h>  // CBLAS 头文件
#include <lapacke.h>  // LAPACKE 头文件

using namespace std;

// 生成 N 维随机方阵
void generate_random_matrix(std::vector<double>& A, int N) {
    std::default_random_engine generator(time(0));
    std::uniform_real_distribution<double> distribution(-1.0, 1.0);

    A.resize(N * N);
    for (int i = 0; i < N * N; ++i) {
        A[i] = distribution(generator);
    }
}


// 求方阵的逆阵
bool inverse_matrix(std::vector<double>& A, int N) {
    std::vector<int> ipiv(N);  // 用于记录主元的位置
    int info = LAPACKE_dgetrf(LAPACK_ROW_MAJOR, N, N, A.data(), N, ipiv.data());  // LU 分解
    if (info != 0) {
        std::cerr << "LU decomposition failed.\n";
        return false;
    }

    info = LAPACKE_dgetri(LAPACK_ROW_MAJOR, N, A.data(), N, ipiv.data());  // 求逆阵
    if (info != 0) {
        std::cerr << "Failed to compute inverse matrix.\n";
        return false;
    }

    return true;
}


// 输出矩阵  前 K 行 K 列
void print_matrix(const std::vector<double>& A, int N, int K) {
    for (int i = 0; i < K; ++i) {
        for (int j = 0; j < K; ++j) {
            std::cout << std::setw(10) << std::fixed << std::setprecision(4) << A[i * N + j];
        }
        std::cout << "\n";
    }
}

int main() 
{
    int N = 1500;  // 矩阵的维度
    std::cout << "N =   " << N << endl;

    double start;
    std::vector<double> A(N * N);  // 随机生成的 N 维方阵
    std::vector<double> B(N * N);  // 拷贝 A 的矩阵 B
    std::vector<double> E(N * N);  // 单位阵

    std::cout << "正在计算 ... ..." << endl;
    start = clock(); //初始时间
    
    // 生成随机方阵 A
    generate_random_matrix(A, N);

    // 拷贝 A 到 B
    B = A;

    // 求 B 的逆阵并返回 逆阵 B。 使用 LAPACK
    std::cout << "Step 1:: 逆阵求解中 ... ...\n";
    inverse_matrix(B, N);

    // 计算 A * B，得到单位阵 E，使用 CBLAS
    std::cout << "Step 2:: 矩阵相乘中 ... ...\n";
    cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
                N, N, N, 1.0, A.data(), N, B.data(), N, 0.0, E.data(), N);

    //截止时间
    std::cout << "\nTime = " << clock() - start << " MilliSeconds\n" << endl;    

    // 打印单位阵 E 的前 10 行 10 列
    cout << "Matrix A_Inv * A = E 单位阵:" << endl;
    print_matrix(E, N, 10);

    return 0;
}[/i][i]