如何在 Fortran mex 文件中使用 BLAS 函数

How to use BLAS functions in Fortran mex files

我曾尝试在我的 Fortran mex 文件中使用 BLAS 函数,但没有成功。这是我使用 DGEMM 的代码之一的示例:

#include "fintrf.h"

C     Gateway subroutine
  subroutine mexfunction(nlhs, plhs, nrhs, prhs)

C     Declarations
  implicit none

C     mexFunction arguments:
  mwPointer plhs(*), prhs(*)
  integer nlhs, nrhs

C     Function declarations:
  mwPointer mxGetPr
  mwPointer mxCreateDoubleMatrix
  mwsize    mxGetM,mxGetN

C     Pointers to input/output mxArrays:
  mwPointer pr_A, pr_B, pr_C
  mwSize :: sizea,sizeb

C     Array information:

  real*8, allocatable :: A(:,:),B(:,:),C(:,:)

C     Get the size of the input array.
  sizea = mxGetM(prhs(1))
  sizeb = mxGetN(prhs(2))

  allocate(A(sizea,sizea),B(sizea,sizeb),C(sizea,sizeb))

C     Create Fortran array from the input argument.
  pr_A = mxGetPr(prhs(1))
  pr_B = mxGetPr(prhs(2))

  call mxCopyPtrToReal8( pr_A, A, sizea**2 )
  call mxCopyPtrToReal8( pr_B, B, sizea*sizeb )

  call MUL(A,B,C,sizea,sizeb)

C     Create matrix for the return argument.
  plhs(1) = mxCreateDoubleMatrix(sizea, sizeb, 0)

  pr_C = mxGetPr(plhs(1))

  call mxCopyReal8ToPtr(C , pr_C, sizea*sizeb )

  return
  end

  subroutine MUL(A,B,C,sizea,sizeb)

  implicit none

  mwSize :: sizea,sizeb
  real*8 :: A(sizea,sizea),B(sizea,sizeb),C(sizea,sizeb),alpha,beta
  integer*4 :: M,N,K
  character ch1, ch2

  ch1 = 'N'
  ch2 = 'N'
  M=size(A,1)
  N=size(B,2)
  K=size(A,2)
  Alpha=1.
  Beta=0.

  CALL DGEMM(ch1,ch2,M,N,K,ALPHA,A,M,B,K,BETA,C,M)


  return
  end subroutine MUL

我使用以下行创建一个 mex 文件:

mex -lmwblas Test.F

mex 文件构建没有任何错误,但是当我尝试使用该函数时,matlab 关闭且没有任何错误消息。我正在使用 Matlab R2016a 和 Intel Visual Fortran Composer XE 2013Microsoft Visual Studio 2012.

您没有正确调用 DGEMM。你所有的论点都是 real*8 但其中许多应该是 integers.

彻底检查 DGEMM 的文档并检查你的调用参数是否完全遵循接口。我非常确信为 K 传递值 1 而不是 sizea 也是错误的。但实际上,你必须检查每一个参数。

我建议在纯 Fortran 中测试您的子例程,只有在它们正确工作后才能从 Matlab 中使用它们。

解决方法如下:

#include "fintrf.h"

C     Gateway subroutine
  subroutine mexfunction(nlhs, plhs, nrhs, prhs)

C     Declarations
  implicit none

C     mexFunction arguments:
  mwPointer plhs(*), prhs(*)
  integer nlhs, nrhs

C     Function declarations:
  mwPointer mxGetPr
  mwPointer mxCreateDoubleMatrix
  mwsize    mxGetM,mxGetN

C     Pointers to input/output mxArrays:
  mwPointer pr_A, pr_B, pr_C
  mwSignedIndex :: sizea,sizeb

C     Array information:

  real*8, allocatable :: A(:,:),B(:,:),C(:,:)

C     Get the size of the input array.
  sizea = mxGetM(prhs(1))
  sizeb = mxGetN(prhs(2))

  allocate(A(sizea,sizea),B(sizea,sizeb),C(sizea,sizeb))

C     Create Fortran array from the input argument.
  pr_A = mxGetPr(prhs(1))
  pr_B = mxGetPr(prhs(2))

  call mxCopyPtrToReal8( pr_A, A, sizea**2 )
  call mxCopyPtrToReal8( pr_B, B, sizea*sizeb )

  call MUL(A,B,C,sizea,sizeb)

C     Create matrix for the return argument.
  plhs(1) = mxCreateDoubleMatrix(sizea, sizeb, 0)

  pr_C = mxGetPr(plhs(1))

  call mxCopyReal8ToPtr(C , pr_C, sizea*sizeb )

  return
  end

  subroutine MUL(A,B,C,sizea,sizeb)

  implicit none

  mwSignedIndex :: sizea,sizeb
  real*8 :: A(sizea,sizea),B(sizea,sizeb),C(sizea,sizeb),alpha,beta
  mwSignedIndex :: M,N,K
  character ch1, ch2

  ch1 = 'N'
  ch2 = 'N'
  M=size(A,1)
  N=size(B,2)
  K=size(A,2)
  Alpha=1.0
  Beta=0.0
  CALL DGEMM(ch1,ch2,M,N,K,ALPHA,A,M,B,K,BETA,C,M)


  return
  end subroutine MUL

并构建 mex 文件使用

mex -largeArrayDims -lmwblas Test.F