与 LAPACK、BLAS、dgemm、intiger 参数类型的矩阵乘法答案

【问题标题】：Matrix multiplication with LAPACK, BLAS, dgemm, intiger parameter type与 LAPACK、BLAS、dgemm、intiger 参数类型的矩阵乘法
【发布时间】：2013-12-21 00:12:28
【问题描述】：

我正在使用不同的矩阵参数类型测试矩阵乘法的不同选项。其中之一是 BLAS 中的 dgemm 例程。当我想创建一个定义为整数（kind=1）、大小为 1000x1000（nxp）的矩阵时，它与 dgemm 一起崩溃，但 matmul 做得很好。当我将矩阵的大小减小到 500x500 时，两者都运行良好。此外，我将所有矩阵定义为 real(8)，两者都计算了矩阵乘积，但结果不同。我正在使用的代码是：

    program test
    implicit none

    real(8), allocatable :: x(:,:),xi(:,:),xt(:,:)
    integer(kind=1), allocatable :: z(:,:)
    integer :: i,j,n,p
    real(8):: u,myone= 1.d0

    n=1000
    p=1000
    allocate(x(n,n),z(n,p),xi(n,n),xt(n,n))
    do i=1,n
      do j=1,p
        call random_number(u)
        z(i,j)=real(floor(u*3),8)
      enddo
    enddo

    print*,"matmul"
    x=matmul(z,transpose(z))

    do i=1,min(10,n)
       write(*,'(10(g10.3,x))') x(i,1:min(10,n))
    enddo 

    print*,"dgemm"
    call dgemm('n'   ,'t'   ,n,n,p,myone  ,Z,n  ,Z,n  ,myone ,X,n)

    do i=1,min(10,n)
       write(*,'(10(g10.3,x))') x(i,1:min(10,n))
    enddo 

    end program test

我用运行以下代码的 make 语句编译代码（我将其命名为：Makefile）：

    f90=ifort
    optf90=-O0 -heap-arrays
    optdir=-I
    mkl=-L/opt/intel/mkl/lib -lmkl_intel_lp64 -lmkl_intel_thread -lmkl_core -openmp -lpthread
    prog=test
    dir= .
    a.out:  $(prog).o
            $(f90)  $(optf90) \
                    $(prog).o  \
                    $(mkl)  $(libf77)

    $(prog).o:      $(prog).f90
            $(f90)  $(optdir)$(dir)  -c $(optf90) $(prog).f90

有谁知道 dgemm 例程对于定义为整数的大矩阵有什么问题，以及 matmul/dgemm 产生不同结果的可能原因是什么？

【问题讨论】：

标签： fortran fortran90 matrix-multiplication lapack blas

【解决方案1】：

DGEMM 适用于double precision 实数数字，而不适用于整数（任何类型）。

如果您在使用integer 数字和DGEMM 时得到任何（正确）结果，我会感到惊讶。

另一方面，MATMUL 接受 integers 作为输入。

【讨论】：