Hi,
I have a Fortran mex file that has been running fine on OSX 64 bit or 32 bit, and also on 32 bit Windows. However, I've tried compiling it on Linux (using gfortran 4.3 or 4.4), and it crashes with a segmentation error. I'm at a total loss at the moment, because surely if there was a problem with 32/64 bit variable types etc, then it would have crashed on the 64 bit mac too??
Can anyone help, and see any obvious mistakes in the code?? I'm really a Linux newbie and I'm really lost here...
#include "fintrf.h"
C=======================================================================
C Abeles Calculation for Matlab
C
C out = abeles(xdata,sld,nbair,nbsub,ssub,repeats,s_between_repeats,resol)
C
C=======================================================================
subroutine mexfunction(nlhs, plhs, nrhs, prhs)
C
C
mwpointer plhs(*), prhs(*)
mwpointer mxGetPr, mxCreateDoubleMatrix
mwpointer x_pr, y_pr
C
C
integer pr_in, pr_in2, pr_in3, pr_in4, pr_in5, pr_in6, pr_in7
integer pr_out, pr_in8
mwsize mxGetM, mxGetN
mwsize m_in_q, n_in_q, points, layers
mwsize m_in_sld, n_in_sld
if(nrhs .ne. 8) then
call mexErrMsgTxt('8 inputs required to Abeles.')
endif
if(nlhs .gt. 1) then
call mexErrMsgTxt('Less than one output required to Abeles.')
endif
! work out points
m_in_q = mxGetM(prhs(1))
n_in_q = mxGetN(prhs(1))
if (n_in_q .gt. 1) then
call mexErrMsgTxt('Qz needs to be a column vector')
endif
points = m_in_q
! do the same for layers
m_in_sld = mxGetM(prhs(2))
n_in_sld = mxGetN(prhs(2))
if (n_in_sld .ne. 3) then
call mexErrMsgTxt('SLD needs three columns: d r s')
endif
layers = m_in_sld
pr_in = mxGetPr(prhs(1))
pr_in2 = mxGetPr(prhs(2))
pr_in3 = mxGetPr(prhs(3))
pr_in4 = mxGetPr(prhs(4))
pr_in5 = mxGetPr(prhs(5))
pr_in6 = mxGetPr(prhs(6))
pr_in7 = mxGetPr(prhs(7))
pr_in8 = mxGetPr(prhs(8))
plhs(1) = mxCreateDoubleMatrix(m_in_q, n_in_q, 0)
pr_out = mxGetPr(plhs(1))
C Call the computational routine.
call comp(%val(pr_out),%val(pr_in),%val(pr_in2),%val(pr_in3),
+ %val(pr_in4),%val(pr_in5),%val(pr_in6),
+ %val(pr_in7),%val(pr_in8),layers,points)
return
end
C=========================================================================
c
C computational routine.....
C
C=========================================================================
subroutine comp(out,x,sld,nbair,nbsub,ssub,nrepeats,resType,
+ resol,layers,points)
integer size, i
integer layers,points,j,reploop, nl, loop, nrepeats
double complex c0 , c1 , ci , beta , rij , N , M , MI , rimaj
double complex blast , num , den , quo
double complex pimag , pair , psub , plast
real*8 pi , twopi , lam , ilow, ihi
real*8 thick , q , theta, rho, rough, resol
real*8 preal , snair , snsub , snlay, rfac, resType
real*8 out(points,1),x(points,1),sld(layers,3),nbair,nbsub,ssub
real*8 dummydata(points,1), g(points,1)
dimension N(2,2) , M(2,2) , MI(2,2)
c0 = Dcmplx(0 , 0)
c1 = Dcmplx(1 , 0)
ci = Dcmplx(0 , 1)
lam = 4.5
pi = 3.141592653589
rfac = ((4*pi)*(4*pi))/2
twopi = 2*pi
snair = (1.0  dble(nbair*((lam*lam) / twopi)))
snsub = (1.0  dble(nbsub*((lam*lam) / twopi)))
do loop = 1 , points
q = x(loop,1)
theta = asin(q*lam / (4*pi))
preal = ((snsub)*(snsub))  ((snair)*(snair))*(cos(theta)**2)
psub = cdsqrt(preal*c1)
pair = snair*sin(theta)*c1
plast = pair
blast = 0.0
rlast = sld(1,3)
MI(1 , 1) = c1
MI(2 , 2) = c1
MI(1 , 2) = c0
MI(2 , 1) = c0
do reploop = 1, nrepeats
do nl = 1 , 3
thick = 1000
rho = 20e6
rough = 1
snlay = (1  dble(rho*((lam*lam) / twopi)))
preal = (snlay*snlay)  (snair*snair) *cos(theta)**2
pimag = cdsqrt(preal*c1)
beta = (twopi / lam)*thick*pimag
rij = dcmplx(plast  pimag) / dcmplx(plast + pimag)
rij = rij * cdexp(rfac*plast*pimag*(rough*rough)/(lam*lam))
N(1 , 1) = cdexp(blast*ci)
N(2 , 1) = rij * cdexp(  blast*ci)
N(2 , 2) = cdexp(  blast*ci)
N(1 , 2) = rij * cdexp(blast*ci)
plast = pimag
blast = beta
rlast = rough
call mcopy(MI , M)
call mmult(M , N , MI)
end do
end do
rij = (dcmplx(plast  psub)) / (dcmplx(plast + psub))
rij = rij * cdexp(rfac*plast*psub*(ssub*ssub)/(lam*lam))
N(1 , 1) = cdexp(blast*ci)
N(2 , 1) = rij*cdexp(  blast*ci)
N(2 , 2) = cdexp(  blast*ci)
N(1 , 2) = rij*cdexp(blast*ci)
call mcopy(MI , M)
call mmult(M , N , MI)
num = MI(2 , 1)*dconjg(MI(2 , 1))
den = MI(1 , 1)*dconjg(MI(1 , 1))
call cdiv(num , den , quo)
out(loop,1) = quo
end
return
end subroutine
!  
SUBROUTINE CDIV(NUM , DEN , QUO)
! SUBROUTINE TO DIVIDE COMPLEX NUMBERS
! AND TRAP DIVIDE CHECKS
double COMPLEX NUM , DEN , QUO , C1 , C0
C1 = DCMPLX(1.0 , 0.0)
C0 = DCMPLX(0.0 , 0.0)
IF(DEN.EQ.NUM)THEN
QUO = C1
RETURN
ENDIF
IF(DEN.EQ.C0 )THEN
QUO = C0
RETURN
ENDIF
QUO = NUM / DEN
RETURN
END SUBROUTINE
!  
SUBROUTINE MMULT(A1 , A2 , A3)
! SUBROUTINE FOR MATRIX MUTIPLICATION
DOUBLE COMPLEX A1(2 , 2) , A2(2 , 2) , A3(2 , 2)
A3(1 , 1) = A1(1 , 1)*A2(1 , 1) + A1(1 , 2)*A2(2 , 1)
A3(1 , 2) = A1(1 , 1)*A2(1 , 2) + A1(1 , 2)*A2(2 , 2)
A3(2 , 1) = A1(2 , 1)*A2(1 , 1) + A1(2 , 2)*A2(2 , 1)
A3(2 , 2) = A1(2 , 1)*A2(1 , 2) + A1(2 , 2)*A2(2 , 2)
RETURN
END
SUBROUTINE MCOPY(A1 , A2)
! SUBROUTINE TO COPY A MATRIX
DOUBLE COMPLEX A1(2 , 2) , A2(2 , 2)
A2(1 , 1) = A1(1 , 1)
A2(1 , 2) = A1(1 , 2)
A2(2 , 1) = A1(2 , 1)
A2(2 , 2) = A1(2 , 2)
RETURN
END SUBROUTINE
