mex fortran allocate

Discussion:

(too old to reply)

Nathaniel

2011-10-13 18:49:27 UTC

I recently came across this example on how to allocate arrays in Fortran. Up until I've been hard-coding the dimensions, which is works fine except when you want to change them. The example below works on a 64-bit windows machine with 2010a repeatedly, it works once on a 64-bit windows machine with 2011a (MATLAB freezes without an error message on the second execution) and it does not work at all on a 32-bit windows machine with 2011a (MATLAB crashes with generic Windows "program has stopped working" message, not freezes). It is my understanding that allocate should work with mex. I wonder, is this a bug in 2011a?

This is the script I run:
********************************************************************\

mex xf_allocatable_5D.F

x = rand(5,5,5,5,10);
y = xf_allocatable_5D(x);

********************************************************************

This is the Fortran code:
********************************************************************
#include "fintrf.h"
C ALLOCATABLE arrays are Fortran's version of malloc/free or new/delete.

C Copyright 1984-2010 MathWorks, Inc.
C $Revision: 1.1.6.1 $

module work_array
mwSize, dimension(:), allocatable :: m
real*8, dimension(:,:,:,:,:), allocatable :: x
real*8, dimension(:,:,:,:,:), allocatable :: y
end module

subroutine timestwo()
use work_array
implicit none
mwSize i,j,k,l,p
C
do i = 1,m(1)
do j = 1,m(2)
do k = 1,m(3)
do l = 1,m(4)
do p = 1,m(5)
y(i,j,k,l,p) = 2.0 * x(i,j,k,l,p)
end do
end do
end do
end do
end do
return
end

subroutine mexFunction(nlhs, plhs, nrhs, prhs)
use work_array
implicit none
mwPointer plhs(*), prhs(*)
integer nlhs, nrhs

integer*4 mxClassIDFromClassName
mwPointer mxGetPr, mxCreateNumericArray, mxGetDimensions
mwPointer x_pr, y_pr, m_pr

integer mxIsNumeric
mwSize mxGetM, mxGetN, mxGetNumberofDimensions
mwSize size, mdim
integer*4 myclassid
integer status

C Check for proper number of arguments.
if(nrhs .ne. 1) then
call mexErrMsgTxt('One input required.')
elseif(nlhs .ne. 1) then
call mexErrMsgTxt('One output required.')
endif

C Get the size of the input array.
m_pr = mxGetDimensions(prhs(1))
mdim = mxGetNumberofDimensions(prhs(1))
allocate(m(mdim))
call mxCopyPtrToReal8(m_pr,m,5)
size = m(1)*m(2)*m(3)*m(4)*m(5)

C Check to insure the input is a number.
if(mxIsNumeric(prhs(1)) .eq. 0) then
call mexErrMsgTxt('Input must be a number.')
endif

C Create matrix for the return argument.
myclassid = mxClassIDFromClassName('double')
plhs(1) = mxCreateNumericArray(5,m,myclassid,0)
x_pr = mxGetPr(prhs(1))
y_pr = mxGetPr(plhs(1))

C Allocate work arrays
allocate(x(m(1),m(2),m(3),m(4),m(5)), STAT=status)
allocate(y(m(1),m(2),m(3),m(4),m(5)), STAT=status)

C Call the computational subroutine.
call mxCopyPtrToReal8(x_pr,x,size)
call timestwo()
call mxCopyReal8ToPtr(y,y_pr,size)

C Free work arrays
deallocate(x)
deallocate(y)
deallocate(m)

return
end
*******************************************************************

As an experiment to try to get it to work on the 32-bit windows system with 2011a I hard-coded the vector m. This runs once and on the second time MATLAB freezes. I also tested it on the 64-bit 2010a and 64-bit 2011a and get the same result. None of the three systems can run this function a second time.

Here is the Fortran code (used the same script above to mex and call the function):
********************************************************************

#include "fintrf.h"
C ALLOCATABLE arrays are Fortran's version of malloc/free or new/delete.

C Copyright 1984-2010 MathWorks, Inc.
C $Revision: 1.1.6.1 $

module work_array
mwSize :: m(5)
real*8, dimension(:,:,:,:,:), allocatable :: x
real*8, dimension(:,:,:,:,:), allocatable :: y
end module

subroutine timestwo()
use work_array
implicit none
mwSize i,j,k,l,p
C
do i = 1,m(1)
do j = 1,m(2)
do k = 1,m(3)
do l = 1,m(4)
do p = 1,m(5)
y(i,j,k,l,p) = 2.0 * x(i,j,k,l,p)
end do
end do
end do
end do
end do
return
end

subroutine mexFunction(nlhs, plhs, nrhs, prhs)
use work_array
implicit none
mwPointer plhs(*), prhs(*)
integer nlhs, nrhs

integer*4 mxClassIDFromClassName
mwPointer mxGetPr, mxCreateNumericArray, mxGetDimensions
mwPointer x_pr, y_pr, m_pr

integer mxIsNumeric
mwSize mxGetM, mxGetN, mxGetNumberofDimensions
mwSize size, mdim
integer*4 myclassid
integer status

C Check for proper number of arguments.
if(nrhs .ne. 1) then
call mexErrMsgTxt('One input required.')
elseif(nlhs .ne. 1) then
call mexErrMsgTxt('One output required.')
endif

C Get the size of the input array.
m_pr = mxGetDimensions(prhs(1))
call mxCopyPtrToReal8(m_pr,m,5)
size = m(1)*m(2)*m(3)*m(4)*m(5)

C Check to insure the input is a number.
if(mxIsNumeric(prhs(1)) .eq. 0) then
call mexErrMsgTxt('Input must be a number.')
endif

C Create matrix for the return argument.
myclassid = mxClassIDFromClassName('double')
plhs(1) = mxCreateNumericArray(5,m,myclassid,0)
x_pr = mxGetPr(prhs(1))
y_pr = mxGetPr(plhs(1))

C Allocate work arrays
allocate(x(m(1),m(2),m(3),m(4),m(5)), STAT=status)
allocate(y(m(1),m(2),m(3),m(4),m(5)), STAT=status)

C Call the computational subroutine.
call mxCopyPtrToReal8(x_pr,x,size)
call timestwo()
call mxCopyReal8ToPtr(y,y_pr,size)

C Free work arrays
deallocate(x)
deallocate(y)

return
end
********************************************************************

Please let me know if you have any ideas on what could be happening that prevents a successful second call of this function.

Much appreciation,
Nate

TideMan

2011-10-13 22:07:38 UTC

Permalink

Post by Nathaniel
I recently came across this example on how to allocate arrays in Fortran. Up until I've been hard-coding the dimensions, which is works fine except when you want to change them. The example below works on a 64-bit windows machine with 2010a repeatedly, it works once on a 64-bit windows machine with 2011a (MATLAB freezes without an error message on the second execution) and it does not work at all on a 32-bit windows machine with 2011a (MATLAB crashes with generic Windows "program has stopped working" message, not freezes). It is my understanding that allocate should work with mex. I wonder, is this a bug in 2011a?
********************************************************************\
mex xf_allocatable_5D.F
x = rand(5,5,5,5,10);
y = xf_allocatable_5D(x);
********************************************************************
********************************************************************
#include "fintrf.h"
C ALLOCATABLE arrays are Fortran's version of malloc/free or new/delete.
C Copyright 1984-2010 MathWorks, Inc.
C $Revision: 1.1.6.1 $
module work_array
mwSize, dimension(:), allocatable :: m
real*8, dimension(:,:,:,:,:), allocatable :: x
real*8, dimension(:,:,:,:,:), allocatable :: y
end module
subroutine timestwo()
use work_array
implicit none
mwSize i,j,k,l,p
C
do i = 1,m(1)
do j = 1,m(2)
do k = 1,m(3)
do l = 1,m(4)
do p = 1,m(5)
y(i,j,k,l,p) = 2.0 * x(i,j,k,l,p)
end do
end do
end do
end do
end do
return
end
subroutine mexFunction(nlhs, plhs, nrhs, prhs)
use work_array
implicit none
mwPointer plhs(*), prhs(*)
integer nlhs, nrhs
integer*4 mxClassIDFromClassName
mwPointer mxGetPr, mxCreateNumericArray, mxGetDimensions
mwPointer x_pr, y_pr, m_pr
integer mxIsNumeric
mwSize mxGetM, mxGetN, mxGetNumberofDimensions
mwSize size, mdim
integer*4 myclassid
integer status
C Check for proper number of arguments.
if(nrhs .ne. 1) then
call mexErrMsgTxt('One input required.')
elseif(nlhs .ne. 1) then
call mexErrMsgTxt('One output required.')
endif
C Get the size of the input array.
m_pr = mxGetDimensions(prhs(1))
mdim = mxGetNumberofDimensions(prhs(1))
allocate(m(mdim))
call mxCopyPtrToReal8(m_pr,m,5)
size = m(1)*m(2)*m(3)*m(4)*m(5)
C Check to insure the input is a number.
if(mxIsNumeric(prhs(1)) .eq. 0) then
call mexErrMsgTxt('Input must be a number.')
endif
C Create matrix for the return argument.
myclassid = mxClassIDFromClassName('double')
plhs(1) = mxCreateNumericArray(5,m,myclassid,0)
x_pr = mxGetPr(prhs(1))
y_pr = mxGetPr(plhs(1))
C Allocate work arrays
allocate(x(m(1),m(2),m(3),m(4),m(5)), STAT=status)
allocate(y(m(1),m(2),m(3),m(4),m(5)), STAT=status)
C Call the computational subroutine.
call mxCopyPtrToReal8(x_pr,x,size)
call timestwo()
call mxCopyReal8ToPtr(y,y_pr,size)
C Free work arrays
deallocate(x)
deallocate(y)
deallocate(m)
return
end
*******************************************************************
As an experiment to try to get it to work on the 32-bit windows system with 2011a I hard-coded the vector m. This runs once and on the second time MATLAB freezes. I also tested it on the 64-bit 2010a and 64-bit 2011a and get the same result. None of the three systems can run this function a second time.
********************************************************************
#include "fintrf.h"
C ALLOCATABLE arrays are Fortran's version of malloc/free or new/delete.
C Copyright 1984-2010 MathWorks, Inc.
C $Revision: 1.1.6.1 $
module work_array
mwSize :: m(5)
real*8, dimension(:,:,:,:,:), allocatable :: x
real*8, dimension(:,:,:,:,:), allocatable :: y
end module
subroutine timestwo()
use work_array
implicit none
mwSize i,j,k,l,p
C
do i = 1,m(1)
do j = 1,m(2)
do k = 1,m(3)
do l = 1,m(4)
do p = 1,m(5)
y(i,j,k,l,p) = 2.0 * x(i,j,k,l,p)
end do
end do
end do
end do
end do
return
end
subroutine mexFunction(nlhs, plhs, nrhs, prhs)
use work_array
implicit none
mwPointer plhs(*), prhs(*)
integer nlhs, nrhs
integer*4 mxClassIDFromClassName
mwPointer mxGetPr, mxCreateNumericArray, mxGetDimensions
mwPointer x_pr, y_pr, m_pr
integer mxIsNumeric
mwSize mxGetM, mxGetN, mxGetNumberofDimensions
mwSize size, mdim
integer*4 myclassid
integer status
C Check for proper number of arguments.
if(nrhs .ne. 1) then
call mexErrMsgTxt('One input required.')
elseif(nlhs .ne. 1) then
call mexErrMsgTxt('One output required.')
endif
C Get the size of the input array.
m_pr = mxGetDimensions(prhs(1))
call mxCopyPtrToReal8(m_pr,m,5)
size = m(1)*m(2)*m(3)*m(4)*m(5)
C Check to insure the input is a number.
if(mxIsNumeric(prhs(1)) .eq. 0) then
call mexErrMsgTxt('Input must be a number.')
endif
C Create matrix for the return argument.
myclassid = mxClassIDFromClassName('double')
plhs(1) = mxCreateNumericArray(5,m,myclassid,0)
x_pr = mxGetPr(prhs(1))
y_pr = mxGetPr(plhs(1))
C Allocate work arrays
allocate(x(m(1),m(2),m(3),m(4),m(5)), STAT=status)
allocate(y(m(1),m(2),m(3),m(4),m(5)), STAT=status)
C Call the computational subroutine.
call mxCopyPtrToReal8(x_pr,x,size)
call timestwo()
call mxCopyReal8ToPtr(y,y_pr,size)
C Free work arrays
deallocate(x)
deallocate(y)
return
end
********************************************************************
Please let me know if you have any ideas on what could be happening that prevents a successful second call of this function.
Much appreciation,
Nate

Whenever I've had this problem, it is because there is something that
has not been deallocated, but I cannot see where this is a problem in
your code.
Have you tried doing a clear fun before calling the 2nd time?

BTW, you don't need all those do loops
y=2*x
will do it.

Nathaniel

2011-10-13 23:35:26 UTC

Permalink

Post by TideMan
Whenever I've had this problem, it is because there is something that
has not been deallocated, but I cannot see where this is a problem in
your code.
Have you tried doing a clear fun before calling the 2nd time?
BTW, you don't need all those do loops
y=2*x
will do it.

In fact, I have tried clear functions after the first call and this also crashes MATLAB. I believe everything that has been allocated is then deallocated and is consistent with other examples out there.

The loops are a place holder for a more complicated routine. This is just as simple example where I am trying to isolate the underlying cause of the freeze/crash.

Please give it a try on your system and see if you can replicate the bug. Thanks for taking a look!

James Tursa

2011-10-14 07:30:31 UTC

Permalink

Alas ... yet another reminder that I am extremely deficient in getting the 64-bit changes to my Fortran 95 interface package uploaded to the FEX. I will point out this package to you with the caveat that the examples provided in the package, and possibly some of the functions, will not work in 64-bit. Feel free to contact me on the side if you want a "beta" version that hasn't been fully tested yet:

http://www.mathworks.com/matlabcentral/fileexchange/25934-fortran-95-interface-to-matlab-api-with-extras

The timestwo example can be accomplished in just a few lines if you use my package. E.g., the following functions are provided:

fpGetPr : Gets a Fortran pointer that points directly at the data area of the input variable (no copy is made). The result can be used just like a regular Fortran array.

fpGetPrCopy : Gets a Fortran pointer to a *copy* of the data of the input variable. The result can be used just like a regular Fortran array. Uses the MATLAB API allocating function mxMalloc behind the scenes, so the memory is completely known to the MATLAB memory manager.

fpDeallocate : Deallocates memory allocated with fpGetPrCopy. But even if your mex function exits with an error before you call fpDeallocate the memory will not leak since the MATLAB memory manager knows about it and will garbage collect it for you.

And now on to specific comments about the code you posted ...

Post by Nathaniel
********************************************************************\
mex xf_allocatable_5D.F
x = rand(5,5,5,5,10);
y = xf_allocatable_5D(x);
********************************************************************
********************************************************************
#include "fintrf.h"
C ALLOCATABLE arrays are Fortran's version of malloc/free or new/delete.
C Copyright 1984-2010 MathWorks, Inc.
C $Revision: 1.1.6.1 $
module work_array
mwSize, dimension(:), allocatable :: m
real*8, dimension(:,:,:,:,:), allocatable :: x
real*8, dimension(:,:,:,:,:), allocatable :: y
end module
subroutine timestwo()
use work_array
implicit none
mwSize i,j,k,l,p
C
do i = 1,m(1)
do j = 1,m(2)
do k = 1,m(3)
do l = 1,m(4)
do p = 1,m(5)
y(i,j,k,l,p) = 2.0 * x(i,j,k,l,p)
end do
end do
end do
end do
end do
return
end
subroutine mexFunction(nlhs, plhs, nrhs, prhs)
use work_array
implicit none
mwPointer plhs(*), prhs(*)
integer nlhs, nrhs
integer*4 mxClassIDFromClassName
mwPointer mxGetPr, mxCreateNumericArray, mxGetDimensions
mwPointer x_pr, y_pr, m_pr
integer mxIsNumeric
mwSize mxGetM, mxGetN, mxGetNumberofDimensions
mwSize size, mdim
integer*4 myclassid
integer status
C Check for proper number of arguments.
if(nrhs .ne. 1) then
call mexErrMsgTxt('One input required.')
elseif(nlhs .ne. 1) then
call mexErrMsgTxt('One output required.')
endif
C Get the size of the input array.
m_pr = mxGetDimensions(prhs(1))
mdim = mxGetNumberofDimensions(prhs(1))
allocate(m(mdim))
call mxCopyPtrToReal8(m_pr,m,5)

Fortran does not automatically type-promote arguments like C does. So don't hard code a literal constant into the API calls like you do with the 5 above ... it may not match mwSize in bit length. Better to use a variable or parameter typed to mwSize, set it to 5, and pass that in the argument list. Also, I am not sure why the mxCopyPtrToReal8 routine is used here ... this copies 8 bytes per element which may not match the number of bytes of an mwSize type. I know it won't on a 32-bit system (will copy into invalid memory and cause a crash) and may do the same on a 64-bit system as well. So you will need additional checks here to only copy the number of bytes that m actually has available based on the number of bytes of an mwSize type. E.g., mxCopyPtrToInteger4.

Post by Nathaniel
size = m(1)*m(2)*m(3)*m(4)*m(5)

Poor choice of variable name, size, since this matches a Fortran intrinsic function.

Post by Nathaniel
C Check to insure the input is a number.
if(mxIsNumeric(prhs(1)) .eq. 0) then
call mexErrMsgTxt('Input must be a number.')
endif
C Create matrix for the return argument.
myclassid = mxClassIDFromClassName('double')
plhs(1) = mxCreateNumericArray(5,m,myclassid,0)

Same comment as above re the hard-coded 5.

Post by Nathaniel
x_pr = mxGetPr(prhs(1))
y_pr = mxGetPr(plhs(1))
C Allocate work arrays
allocate(x(m(1),m(2),m(3),m(4),m(5)), STAT=status)
allocate(y(m(1),m(2),m(3),m(4),m(5)), STAT=status)
C Call the computational subroutine.
call mxCopyPtrToReal8(x_pr,x,size)
call timestwo()
call mxCopyReal8ToPtr(y,y_pr,size)
C Free work arrays
deallocate(x)
deallocate(y)
deallocate(m)
return
end

(snip)

Post by Nathaniel
Please let me know if you have any ideas on what could be happening that prevents a successful second call of this function.
Much appreciation,
Nate

James Tursa

Nathaniel

2011-10-14 14:47:15 UTC

Permalink

Thanks James. As usual, you are very helpful.

I made the changes 1) using mxCopyPtrToInteger4 instead of mxCopyPtrToReal8 (for m_pr to m) and 2) first assigning dimensions to a variable mdim and then passing that into mxCopyPtrToInteger4 and mxCreateNumericArray (MATLAB spoils me when it comes to such subtleties). Now the example works repeatedly on both the 2010a 64-bit and 2011a 32-bit. They even work inside parfor loops as well, which is important for some of my applications.

However, I tried making these changes (for 3 dimensions, not 5) in the actual function I am working with and I get different errors. The function below runs once without an error, but on the second execution it gives the following errors:

32-bit
severe (179): Cannot allocate array - overflow on array size calculation

64-bit
forrtl: severe (41): insufficient virtual memory
Stack trace terminated abnormally

Which are Fortran errors, not MATLAB (not a seg error) or windows errors. I tried to clear the function between executions to no avail. When I rerun the script (with a clear/clc at the beginning) and it only calls the function once it does not produce the error. I should be very close. Structurally, the primary difference between this code and the example is that it doesn't initialize the allocatable arrays in work_array. Does this make a difference?

The function is below:
************************************************************************
#include "fintrf.h"
subroutine mexFunction(nlhs, plhs, nrhs, prhs)

! Declarations
implicit none

! mexFunction argument
mwPointer plhs(*), prhs(*)
integer*4 nlhs, nrhs

! Function declarations:
mwSize mxGetNumberofDimensions
mwpointer mxGetPr, mxCreateNumericArray, mxGetDimensions
double precision mxGetScalar
integer*4 mxClassIDFromClassName

! Pointers to input/output mxArrays
mwpointer z1i_pr,z2i_pr
mwpointer z1_bnd_pr,z2_bnd_pr
mwpointer T_pr,P_pr
mwpointer A1_pr,A2_pr
mwpointer f1_pr,f2_pr
mwpointer m_pr,n_pr

! Array information
mwSize :: m(3),n(2)
mwSize mdim,mmax,ndim
mwSize coefs,nodes
integer*4 myclassid
double precision, allocatable, dimension(:,:,:) :: z1i,z2i
double precision, dimension(2) :: z1_bnd,z2_bnd
double precision, allocatable, dimension(:,:) :: T,P
double precision, allocatable, dimension(:,:) :: A1,A2
double precision, allocatable, dimension(:,:,:) :: f1,f2

! Load Inputs
! Points to evaluate
mdim = 3
m_pr = mxGetDimensions(prhs(1))
call mxCopyPtrToInteger4(m_pr,m,mdim)
mmax = max(m(1),m(2))
nodes = m(1)*m(2)*m(3)
allocate(z1i(m(1),m(2),m(3)))
allocate(z2i(m(1),m(2),m(3)))
z1i_pr = mxGetPr(prhs(1))
z2i_pr = mxGetPr(prhs(2))
call mxCopyPtrToReal8(z1i_pr,z1i,nodes)
call mxCopyPtrToReal8(z2i_pr,z2i,nodes)

! Bounds
z1_bnd_pr = mxGetPr(prhs(3))
z2_bnd_pr = mxGetPr(prhs(4))
call mxCopyPtrToReal8(z1_bnd_pr,z1_bnd,2)
call mxCopyPtrToReal8(z2_bnd_pr,z2_bnd,2)

! Chebyshev Polynomial Parameters
allocate(T(mmax,mmax))
allocate(P(mmax,mmax))
T_pr = mxGetPr(prhs(5))
P_pr = mxGetPr(prhs(6))
call mxCopyPtrToReal8(T_pr,T,mmax*mmax)
call mxCopyPtrToReal8(P_pr,P,mmax*mmax)

! Least square weights
ndim = 2
n_pr = mxGetDimensions(prhs(7))
call mxCopyPtrToInteger4(n_pr,n,ndim)
coefs = n(1)*n(2)
allocate(A1(n(1),n(2)))
allocate(A2(n(1),n(2)))
A1_pr = mxGetPr(prhs(7))
A2_pr = mxGetPr(prhs(8))
call mxCopyPtrToReal8(A1_pr,A1,coefs)
call mxCopyPtrToReal8(A2_pr,A2,coefs)

!Create array for return argument
myclassid = mxClassIDFromClassName('double')
allocate(f1(m(1),m(2),m(3)))
allocate(f2(m(1),m(2),m(3)))
plhs(1) = mxCreateNumericArray(mdim,m,myclassid,0)
plhs(2) = mxCreateNumericArray(mdim,m,myclassid,0)
f1_pr = mxGetPr(plhs(1))
f2_pr = mxGetPr(plhs(2))

! Call subroutine for assignment
call allcheb(f1,f2,nodes,z1i,z2i,z1_bnd,z2_bnd,mmax,n(1),n(2),m(1),m(2),m(3),T,P,A1,A2)

!Load fortran array to pointer (output to MATLAB)
call mxCopyReal8ToPtr(f1,f1_pr,nodes)
call mxCopyReal8ToPtr(f2,f2_pr,nodes)

!Deallocate arrays
deallocate(z1i)
deallocate(z2i)
deallocate(T)
deallocate(P)
deallocate(A1)
deallocate(A2)
deallocate(f1)
deallocate(f2)

end subroutine mexFunction

subroutine allcheb(f1,f2,nodes,z1i,z2i,z1_bnd,z2_bnd,zmax,n1,n2,m1,m2,m3,T,P,A1,A2)

implicit none
mwSize :: nodes,zmax,n1,n2,m1,m2,m3
double precision, dimension(m1,m2,m3) :: z1i,z2i,f1,f2
double precision, dimension(2) :: z1_bnd,z2_bnd
double precision, dimension(zmax,zmax) :: T,P
double precision, dimension(n1,n1) :: T1,P1
double precision, dimension(n2,n2) :: T2,P2
double precision, dimension(n1) :: vec1
double precision, dimension(n2) :: vec2
double precision, dimension(n1,n2) :: A1,A2
double precision :: temp1,temp2,x1i,x2i
mwSize :: i1,i2,i3,j1,j2

T1 = T(1:n1,1:n1)
T2 = T(1:n2,1:n2)
P1 = P(1:n1,1:n1)
P2 = P(1:n2,1:n2)
do i1 = 1,m1
do i2 = 1,m2
do i3 = 1,m3
! Transform Z in [a,b] to X in [-1,1]
x1i = 2*(z1i(i1,i2,i3)-z1_bnd(1))/(z1_bnd(2) - z1_bnd(1)) - 1;
x2i = 2*(z2i(i1,i2,i3)-z2_bnd(1))/(z2_bnd(2) - z2_bnd(1)) - 1;

! Evaluate Chebyshev Polynomials
vec1 = sum(T1*x1i**P1,2);
vec2 = sum(T2*x2i**P2,2);

! Evaluate function
temp1 = 0
temp2 = 0
do j1 = 1,n1
do j2 = 1,n2
temp1 = temp1 + A1(j1,j2)*vec1(j1)*vec2(j2)
temp2 = temp2 + A2(j1,j2)*vec1(j1)*vec2(j2)
end do
end do
f1(i1,i2,i3) = temp1
f2(i1,i2,i3) = temp2
end do
end do
end do

end subroutine allcheb

James Tursa

2011-10-14 17:03:10 UTC

Permalink

Post by Nathaniel
Thanks James. As usual, you are very helpful.
I made the changes 1) using mxCopyPtrToInteger4 instead of mxCopyPtrToReal8 (for m_pr to m) and 2) first assigning dimensions to a variable mdim and then passing that into mxCopyPtrToInteger4 and mxCreateNumericArray (MATLAB spoils me when it comes to such subtleties). Now the example works repeatedly on both the 2010a 64-bit and 2011a 32-bit. They even work inside parfor loops as well, which is important for some of my applications.
32-bit
severe (179): Cannot allocate array - overflow on array size calculation
64-bit
forrtl: severe (41): insufficient virtual memory
Stack trace terminated abnormally
Which are Fortran errors, not MATLAB (not a seg error) or windows errors. I tried to clear the function between executions to no avail. When I rerun the script (with a clear/clc at the beginning) and it only calls the function once it does not produce the error. I should be very close. Structurally, the primary difference between this code and the example is that it doesn't initialize the allocatable arrays in work_array. Does this make a difference?
call mxCopyPtrToInteger4(m_pr,m,mdim)

e.g., put the following in your code and see what you get:

integer*4, external :: mexPrintf
integer*4 k
character*80 line
:
(etc)
:
call mxCopyPtrToInteger4(m_pr,m,mdim)
write(line,*) 'm = ',m
k = mexPrintf(line//achar(10))
return

And you might do the same for other size info that you use the mxCopy___ routines with.

James Tursa

James Tursa

2011-10-14 17:04:30 UTC

Permalink

Post by Nathaniel
I made the changes 1) using mxCopyPtrToInteger4 instead of mxCopyPtrToReal8 (for m_pr to m) and 2) first assigning dimensions to a variable mdim and then passing that into mxCopyPtrToInteger4 and mxCreateNumericArray (MATLAB spoils me when it comes to such subtleties).

P.S. Fortran can do this too if you set up a generic interface that accepts various size integers ... that is an upgrage to my Fortran 95 submission that is (alas) not uploaded to the FEX yet.