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

Question

我尝试在我的 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 2013 与 Microsoft Visual Studio 2012 一起使用。

Answer 1

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

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

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

Answer 2

解决方案如下：

#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