如何在 fortran mex 文件中使用 mxCalloc

Question

我必须使用 mxCalloc 函数而不是 mex 文件中的常规分配，以避免 matlab 在使用 dgesv 时崩溃。我尝试了很多方法，但其中 none 行得通。这是示例之一

 #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
      mwPointer mxGetM

C     Pointers to input/output mxArrays:
      mwPointer pr_A, pr_B

C     Array information:

      mwPointer sizea,mxCalloc
      real*8 :: A,B
      character*120 :: line

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

      A=mxCalloc(sizea*sizea,8)
      B=mxCalloc(sizea*sizea,8)

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

      call mxCopyPtrToReal8(pr_A,A,sizea**2)

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

      write(line,*), sizea
      call mexPrintf(line) 
      B=A

      call mxCopyReal8ToPtr(B,pr_B,sizea*sizea)

      return
      end

当我运行这段代码时，我得到以下结果

A = [0.9575 , 0.1576 ; 0.9649 , 0.9706]

test(A) = [0.9575 , 0 ; 0.9649 , 0]

但如果我换行 call mxCopyReal8ToPtr(B,pr_B,sizea*sizea) 至

call mxCopyReal8ToPtr(A,pr_B,sizea*sizea), 结果正确

变量 sizea 等于 2，这是正确的，但我无法访问 A 的任何成员，比如 A(1,1)，这是我遇到的错误：

error #6410: This name has not been declared as an array or a function. [A]

Answer 1

代码中的 A 被声明为单个变量，但 size(A,1) 试图像访问数组一样访问它。您必须明确地告诉 Fortran A 具有数组的形状。快速浏览完其余代码后，我假设 A 和 B 应该是二维的

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

现在我不太了解 mxCalloc 的工作原理以及它如何处理内存分配和与 Fortran 的接口，但可能您还需要将 A 和 B 声明为指针。

real*8, pointer :: A(:,:),B(:,:)

Answer 2

（免责声明：我没有使用 Mex 的经验，所以请按以下代码操作...）

这是@ftiaronsem 和之前 one by Dave. By looking at this tutorial 回答的延续，似乎可以在 mexfunction() 中使用 Fortran 可分配数组。但是，OP 的目的是使用 mxCalloc() 而不是用于内存管理的可分配数组。

根据Matlab的在线手册，mxCalloc()似乎return分配内存的地址为mwpointer（这只是integer*8的别名64 位机器，根据头文件fintrf.h）。所以，我们首先得到地址为

mwpointer iA
iA = mxCalloc( sizea * sizea, 8 )

（B 也类似）。接下来，我们要访问从 iA 开始的内存作为二维 Fortran 数组。这可以通过使用 c_f_pointer()（在 F2003 中）来完成。但是，由于此函数的第一个参数接收 type(c_ptr)，我们按如下方式进行：

use iso_c_binding
type(c_ptr) cA
real*8, pointer :: A(:,:)

cA = transfer( iA, cA )                      !! cast integer*8 to c_ptr
call c_f_pointer( cA, A, [ sizea, sizea ] )  !! init an array pointer with c_ptr

在这条语句之后，我们可以将 A 用作普通的二维 Fortran 数组（例如，我们可以将其传递给 LAPACK 例程）。将这些修改包含在 OP 的代码中会给出以下内容。那你能不能试试看能不能用...?

更新代码（ver1）：

#include "fintrf.h"

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

      use iso_c_binding    !<---

C     Declarations
      implicit none

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

C     Function declarations:
      mwPointer mxGetPr
      mwPointer mxCreateDoubleMatrix
      mwsize    mxGetM, sizea          !<--- changed to mwsize (may not be necessary, though)

C     Pointers to input/output mxArrays:
      mwPointer pr_A, pr_B

C     Array information:

      mwPointer mxCalloc

      mwPointer       :: iA, iB           !<---
      type(c_ptr)     :: cA, cB           !<---
      real*8, pointer :: A(:,:), B(:,:)   !<---

      character*120 :: line

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

      iA = mxCalloc( sizea**2, 8 )   !<---
      iB = mxCalloc( sizea**2, 8 )   !<---

      cA = transfer( iA, cA )        !<---
      cB = transfer( iB, cB )        !<---

      call c_f_pointer( cA, A, [ nsizea, nsizea ] )   !<---
      call c_f_pointer( cB, B, [ nsizea, nsizea ] )   !<---

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

      call mxCopyPtrToReal8( pr_A, A, sizea**2 )

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

      write(line,*), sizea
      call mexPrintf(line) 
      B = A

      call mxCopyReal8ToPtr( B, pr_B, sizea**2 )

      !! may need to call mxFree( iA ) and mxFree( iB ) manually?
      !! (or Matlab may automatically do it upon exit)

      return
      end

如何在 fortran mex 文件中使用 mxCalloc

How to use mxCalloc in fortran mex-file

matlab

fortran

mex