如何在 fortran mex 文件中使用 mxCalloc

Posted 2023-02-16

【中文标题】如何在 fortran mex 文件中使用 mxCalloc

【英文标题】：How to use mxCalloc in fortran mex-file

【发布时间】：2016-10-19 14:41:46

【问题描述】：

我必须在 mex 文件中使用 mxCalloc 函数而不是常规分配，以避免 matlab 在使用 dgesv 时崩溃。 我尝试了很多方法，但都没有奏效。 这是其中一个示例

 #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]

测试(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)，这是我遇到的错误：

错误 #6410：此名称尚未声明为数组或 功能。 [一]

【参考方案1】：

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

这是 @ftiaronsem 和 Dave 之前的 one 答案的一些延续。通过查看这个tutorial，似乎可以在mexfunction() 中使用Fortran 可分配数组。但是，OP 的目的是使用mxCalloc() 而不是可分配的数组进行内存管理。

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

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

（B 也是如此）。接下来，我们要以二维 Fortran 数组的形式访问从 iA 开始的内存。这可以通过使用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

【讨论】：

非常感谢您的帮助，好像没问题

但我还是不能使用 dgesv，你知道吗？

请发布您使用 dgesv 的代码的最小示例。如果不查看实际代码，就不可能知道可能出了什么问题。而且我认为在另一个问题中执行此操作是个好主意，因为您原来的问题“如何在 fortran mex 文件中使用 mxCalloc”似乎已通过此答案解决。

【参考方案2】：

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

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

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

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

【讨论】：

我用过这个语法，建mex文件时没有错误，但是用mex文件matlab会崩溃

你知道matlab在fortran文件的哪一行崩溃了吗？您可以将打印语句插入到 fortran 代码中，或者如果这些语句出于某种原因未在 matlab 中显示，请创建文件以查看它崩溃的位置。你能告诉我们调用这个函数的matlab代码吗？

还尝试在 sizea = mxGetM(prhs(1)) 行之后打印（或写入）sizea 的值。它有正确的价值吗？

这个问题有一些变化，我认为你的 cmets 有参考

如果我使用通常的分配，matlab 崩溃，我确定问题是内存分配。有一个我想在 C 中做的例子：mathworks.com/matlabcentral/fx_files/8006/1/dgesv.c