I’m trying to port some solver code onto GPUs, but unfortunately I’m running into an issue that I’m not really sure how to handle. Specifically I’m running into a problem on the integrate part in that it keeps failing on what seems to be the second integration call. Specifically the error I’m getting is

```
FORTRAN STOP: 0: Block (2,1,1), Thread (1,1,1)
stop08.cu:28: pgf90_stop08: block: [1,0,0], thread: [0,0,0] Assertion `FORTRAN_STOP_STATEMENT` failed.
```

This error occurs in this call, it always seems to fail on the a block that has entered but not left this section (N=50).

```
fi = (U-V)
call grad(fi,X,DU)
do i=1,N
DU(i)=max(DU(i),1.d-100)
enddo
do j=1,N
do i=1,N
xvec1(i)=X(i)
yvec1(i)=DU(i)*Kxold(i,j)
enddo
call integrate(2,X(1),X(N),tmp(j),1)
Upold(j)=U(1)-V(1)+tmp(j)
enddo
```

and the subroutines and functions that are causing it are below (all of them are in the module file)

```
attributes(device) function Diff(mom)
implicit none
real*8 :: Diff,mom
Diff=(mom*mom/dsqrt(1.d0+mom*mom))
return
end function Diff
attributes(device) function funcexp(z)
implicit none
real*8 :: funcexp,z,divdif
funcexp=dexp(divdif(yvec,xvec,N,z,3))*dexp(z)
return
end function funcexp
attributes(device) function funcp1(z)
implicit none
real*8 :: funcp1,divdif,z
funcp1=divdif(yvec1,xvec1,N,z,3)
return
end function funcp1
attributes(device) function funcp0(z)
implicit none
real*8 :: funcp0,divdif,z
funcp0=dexp(divdif(yvec0,xvec0,N,z,3))
return
end function funcp0
attributes(device) function funcxi(z)
implicit none
real*8 :: funcxi,z,divdif
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
attributes(device) subroutine grad(fi,xi,dfi)
implicit none
integer i
real*8 :: fi(N),xi(N),dfi(N)
do i=2,N-1
dfi(i)=(fi(i+1)-fi(i-1))/(xi(i+1)-xi(i-1))
enddo
dfi(1)=(fi(2)-fi(1))/(xi(2)-xi(1))
dfi(N)=(fi(N)-fi(N-1))/(xi(N)-xi(N-1))
return
end subroutine grad
!!!! Interpolating subroutine
attributes(device) function DIVDIF(F,A,NN,X,MM)
implicit none
integer N,NN,MM,M,MPLUS,IX,IY,MID,NPTS,IP,L,ISUB,J,I,MMAX, errorNumber
real*8 :: SUM,DIVDIF
real*8 :: A(NN),F(NN),T(20),D(20),X
LOGICAL EXTRA
DATA MMAX/10/
!
! TABULAR INTERPOLATION USING SYMMETRICALLY PLACED ARGUMENT POINTS.
!
! START. FIND SUBSCRIPT IX OF X IN ARRAY A.
IF( (NN.LT.2) .OR. (MM.LT.1) ) GO TO 20
N=NN
M=MM
MPLUS=M+1
IX=0
IY=N+1
! (SEARCH INCREASING ARGUMENTS.)
1 MID=(IX+IY)/2
IF(X.GE.A(MID)) GO TO 2
IY=MID
GO TO 3
! (IF TRUE.)
2 IX=MID
3 IF(IY-IX.GT.1) GO TO 1
GO TO 7
! COPY REORDERED INTERPOLATION POINTS INTO (T(I),D(I)), SETTING
! *EXTRA* TO TRUE IF M+2 POINTS TO BE USED.
7 NPTS=M+2-MOD(M,2)
IP=0
L=0
GO TO 9
8 L=-L
IF(L.GE.0) L=L+1
9 ISUB=IX+L
IF((1.LE.ISUB).AND.(ISUB.LE.N)) GO TO 10
! (SKIP POINT.)
NPTS=MPLUS
GO TO 11
! (INSERT POINT.)
10 IP=IP+1
T(IP)=A(ISUB)
D(IP)=F(ISUB)
11 IF(IP.LT.NPTS) GO TO 8
EXTRA=NPTS.NE.MPLUS
!
! REPLACE D BY THE LEADING DIAGONAL OF A DIVIDED-DIFFERENCE TABLE, SUP-
! PLEMENTED BY AN EXTRA LINE IF *EXTRA* IS TRUE.
DO 14 L=1,M
IF(.NOT.EXTRA) GO TO 12
ISUB=MPLUS-L
D(M+2)=(D(M+2)-D(M))/(T(M+2)-T(ISUB))
12 I=MPLUS
DO 13 J=L,M
ISUB=I-L
D(I)=(D(I)-D(I-1))/(T(I)-T(ISUB))
I=I-1
13 CONTINUE
14 CONTINUE
!
! EVALUATE THE NEWTON INTERPOLATION FORMULA AT X, AVERAGING TWO VALUES
! OF LAST DIFFERENCE IF *EXTRA* IS TRUE.
SUM=D(MPLUS)
IF(EXTRA) SUM=0.5*(SUM+D(M+2))
J=M
DO 15 L=1,M
SUM=D(J)+(X-T(J))*SUM
J=J-1
15 CONTINUE
DIVDIF=SUM
RETURN
20 errorNumber = 1
errorNumber = 2
DIVDIF=999999999999999999.
101 errorNumber = 3
102 errorNumber = 4
END
subroutine get_time(time_sec)
implicit integer(i-n), double precision(a-h, o-z)
character (len=12) tdum1, tdum2, tdum3
dimension i_date(8)
call date_and_time(tdum1, tdum2, tdum3, i_date)
time_sec = i_date(3)*24.0*3600.0 + i_date(5)*3600.0
time_sec = time_sec + i_date(6)*60.0 + i_date(7)+ i_date(8)/1.d3
return
END subroutine get_time
!!!!------------------------------------------------------------------------------
!!!! Subroutine for integrals
attributes(device) SUBROUTINE integrate(func,a,b,ss,v)
INTEGER KMAX,KMAXP,J,JM
real*8 :: a,b,ss,tol
integer :: func,v
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,v)
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'
print *, "ERROR in subroutine integrate"
print *, "Error on thread", threadIdx%x
print *, "Error on block", blockIdx%x
!Keep running on other ones
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'
print *, "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,v)
INTEGER n
real*8 :: a,b,s
integer :: func,v
INTEGER :: it,j
real*8 :: ddel,del,sum,tnm,x
if (v == 1) then
! print *, "block", blockIdx%x, "has entered halfway"
! print *, "a:", a
! print *, "b:", b
v = 1
endif
Select Case (func)
case (1)
if (n.eq.1) then
s=(b-a)*funcexp(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+funcexp(x)
x=x+ddel
sum=sum+funcexp(x)
x=x+del
end do
s=(s+(b-a)*sum/tnm)/3.d0
endif
case (2)
if (n.eq.1) then
s=(b-a)*funcp1(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+funcp1(x)
x=x+ddel
sum=sum+funcp1(x)
x=x+del
end do
s=(s+(b-a)*sum/tnm)/3.d0
endif
case (3)
if (n.eq.1) then
s=(b-a)*funcp0(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+funcp0(x)
x=x+ddel
sum=sum+funcp0(x)
x=x+del
end do
s=(s+(b-a)*sum/tnm)/3.d0
endif
case (4)
if (n.eq.1) then
s=(b-a)*funcxi(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+funcxi(x)
x=x+ddel
sum=sum+funcxi(x)
x=x+del
end do
s=(s+(b-a)*sum/tnm)/3.d0
endif
case default
return
end select
return
END SUBROUTINE halfway
```

Unfortunately I don’t know how to better pinpoint the cause of this error. I know it’s not an issue with the init and evolution file I’m using since they work fine on the non-GPU version. I’ve included all the files needed to compile and run the code should this be insufficient.

I’m compiling it with this command

`pgfortran -Mcuda=cuda8.0 -o CRAFTest CRAFT_CLEAN.f90 -#`

Apologies for the massive amount of text.

CRAFT.zip (12.2 KB)