Как использовать mxCalloc в mex-файле fortran

Question

Мне нужно использовать функцию mxCalloc вместо регулярного выделения в файле mex, чтобы избежать сбой mlabab при использовании 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]

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), и имеет ошибку, с которой я сталкиваюсь:

ошибка # 6410: Это имя не было объявлено как массив или функция . [A]

Answer 1

(Отказ от ответственности: У меня нет опыта работы с Mex, поэтому, пожалуйста, следующий код как таковой ...)

Это некоторое продолжение ответов по @ftiaronsem и предыдущей one Дэйв. Если посмотреть на это tutorial, кажется, что использовать распределенные массивы Fortran в mexfunction(). Но целью OP является использование mxCalloc() вместо выделенных массивов для управления памятью.

Согласно интерактивному руководству по Matlab, mxCalloc(), кажется, возвращает адрес выделенной памяти, как mwpointer (это просто псевдоним integer*8 на 64-разрядных машинах, в соответствии с файлом заголовка fintrf.h). Таким образом, мы сначала получить адрес в качестве

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

(а так же для B). Затем мы хотим получить доступ к памяти, начиная с iA, в виде массива 2-d 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 как обычный 2-й 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 

Answer 2

A в коде объявляется как одна переменная, но размер (A, 1) пытается получить к нему доступ, как если бы он был массивом. Вам нужно будет прямо указать Fortran, что A имеет форму массива. После того, как быстрый взгляд на остальной части кода я предполагаю, что А и В должны быть 2 мерной

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

Теперь я не знаю слишком много о том, как работает mxCalloc и как она обрабатывает выделение памяти и интерфейс к Fortran, но может быть, вам также нужно объявить A и B как указатели.

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