I apologise for not being clear, I’m afraid I myself am not sure what the problem is either. I’m afraid the smallest example that contains all parts that are causing errors I can post is a tad messy and not particularly small.
LNP is an array that is on the main kernal code.
The subroutines and functions that are involved are listed here.
attributes(device) SUBROUTINE integrate(func,a,b,ss)
INTEGER KMAX,KMAXP,J,JM
real*8 :: a,b,func,ss,tol
EXTERNAL :: func
PARAMETER (KMAX=14, KMAXP=KMAX+1, J=5, JM=J-1, tol=1.d-1) !#CONTROL Should work even for EPS=1e-1
INTEGER k
real*8 :: dss,h(KMAXP),s(KMAXP)
h(1)=1.d0
do k=1,KMAX
call halfway(func,a,b,s(k),k)
if (k.ge.J) then
call pol_interpol(h(k-JM),s(k-JM),J,0.d0,ss,dss)
if (dabs(dss).le.tol*dabs(ss)) return
endif
s(k+1)=s(k)
h(k+1)=h(k)/9.d0
end do
!write(*,*) 'ERROR in subroutine Integrate'
stop
END SUBROUTINE integrate
attributes(device) SUBROUTINE pol_interpol(xa,ya,n,x,y,dy)
INTEGER n,NMAX
real*8 :: dy,x,y,xa(n),ya(n)
PARAMETER (NMAX=10)
INTEGER i,m,ns
real*8 :: den,dif,dift,ho,hp,w,c(NMAX),d(NMAX)
ns=1
dif=dabs(x-xa(1))
do i=1,n
dift=abs(x-xa(i))
if (dift.lt.dif) then
ns=i
dif=dift
endif
c(i)=ya(i)
d(i)=ya(i)
end do
y=ya(ns)
ns=ns-1
do m=1,n-1
do i=1,n-m
ho=xa(i)-x
hp=xa(i+m)-x
w=c(i+1)-d(i)
den=ho-hp
if(den.eq.0.d0)then
!write(*,*) 'ERROR in subroutine interpol'
stop
end if
den=w/den
d(i)=hp*den
c(i)=ho*den
end do
if (2*ns.lt.n-m)then
dy=c(ns+1)
else
dy=d(ns)
ns=ns-1
endif
y=y+dy
end do
return
END SUBROUTINE pol_interpol
attributes(device) SUBROUTINE halfway(func,a,b,s,n)
INTEGER n
real*8 :: a,b,s,func
EXTERNAL :: func
INTEGER :: it,j
real*8 :: ddel,del,sum,tnm,x
if (n.eq.1) then
s=(b-a)*func(0.5d0*(a+b))
else
it=3**(n-2)
tnm=it
del=(b-a)/(3.d0*tnm)
ddel=del+del
x=a+0.5d0*del
sum=0.d0
do j=1,it
sum=sum+func(x)
x=x+ddel
sum=sum+func(x)
x=x+del
end do
s=(s+(b-a)*sum/tnm)/3.d0
endif
return
END SUBROUTINE halfway
attributes(device) function funcxi(z)
implicit none
integer N
parameter(N=50)
real*8 :: funcxi,z,xvec(N),yvec(N),divdif
COMMON/VEC/xvec,yvec
funcxi=dexp(divdif(yvec,xvec,N,z,3))*dexp(4.d0*z)*dexp(z)/dsqrt(1.d0+dexp(2.d0*z))
return
end function funcxi
The main device code that calls this is as such (this is just one of many such calls to funcxi or integrate within this main code.)
attributes(device) subroutine doCalculation
!all variables used are defined way up top including LNP(N)
...
do i=1,N
do j=1,N
xvec(j)=lnp(j)
yvec(j)=dlog(f0old(j))+dlog(Kxold(i,j))
enddo
Y1=Ymin
do j=1,N
if(funcxi(lnp(N+1-j)).ge.1.d-100)then
Y1=lnp(N+1-j)
endif
enddo
call integrate(funcxi,Y1,Ymax,csi_c(i))
enddo
...
I was under the impression device code could call other device code.
All of the functions and subroutines are within the same module file.
I hope this is helpful.
As for what I’m trying to do, I believe (I’ve not written the original version of this code, so there are parts of it I’m not familiar with such as that line) I’m doing both passing the results and trying to pass it as an argument to the subroutines. I believe what this is meant to do for passing it to the subroutine is for instance passing a function to integrate is supposed to integrate that function, instead of having an integrate subroutine for each function we’ll need to do.
I can provide more information if this is still insufficient. Also, thank you for your assistance.