如何在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 和 Visual Studio 2012Microsoft。

Answer 1

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

彻底检查DGEMM的文档，并检查一个又一个的调用参数，以完全遵循界面。我非常确信传递 1 的值而不是 K 的 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