Passing function to another function/subroutine in FORTRAN

Accelerated Computing HPC Compilers Legacy PGI Compilers
#1

I’ve got some code I’m trying to get running on a GPU using CUDAFORTRAN and I’m having some difficulty figuring out why something isn’t working quite right. I’ve got this system that goes 3 layers deep that looks something like this. I can’t tell if it’s one problem or two related problems either.

attributes(device) SUBROUTINE halfway(func,a,b,s,n)
EXTERNAL func
integer :: a,b,s,n
func(a,b,s,n)
…
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.d0z)dexp(z)/dsqrt(1.d0+dexp(2.d0z))
return
end function funcxi

and I’m trying to do funcxi(lnp(N+1-j)) but it doesn’t seem to like that as it puts out the error
Calls from device code to a host function are allowed only in emulation mode (line with funcxi(lnp(…))

And also the same error on the lines that take func as an argument. I’m not sure what the issue is or how to solve it.

I can provide more information should that be useful.

#2

Hi Cattaneo,

Apologies, but I’m not clear on the question. The error:

Calls from device code to a host function are allowed only in emulation mode (line with funcxi(lnp(…))

typically means that you’re calling a function that doesn’t have a device attribute. Though, in the example you posted there’s no “lnp” in the assignment to funccxi, so I’m not sure what “lnp” is.

Also, I’m not clear on the topic " Passing function to another function/subroutine in FORTRAN". Are you meaning procedure pointers or using the result of a function as an argument? Procedure pointers are not currently supported within device code, but given the error your seeing I suspect that you’re trying to pass the result of function as an argument.

Can you post a small reproducing example which illustrates the issue you’re having?

Thanks,
Mat

#3

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.

#4

Device routines do need interfaces so I put the above code into a module and see the following errors:

% nvfortran -c test.cuf
NVFORTRAN-S-0155-Calls from device code to a host function are allowed only in emulation mode - func (test.cuf: 80)
NVFORTRAN-S-0155-Calls from device code to a host function are allowed only in emulation mode - func (test.cuf: 89)
NVFORTRAN-S-0155-Calls from device code to a host function are allowed only in emulation mode - func (test.cuf: 91)
  0 inform,   0 warnings,   3 severes, 0 fatal for halfway
NVFORTRAN-S-0155-Calls from device code to a host function are allowed only in emulation mode - divdif (test.cuf: 106)
  0 inform,   0 warnings,   1 severes, 0 fatal for funcxi

“func” are procedure pointers, which aren’t supported in device code. We’re hoping to add this support in future, but it’s not available at this time. You’ll need to change these to be direct calls to a specific function.

The “divdif” function does not have device routine. I’m not sure if it’s missing from the just the example, or your original version as well.

Hope this helps,
Mat

#5

It is missing from the example, I completely missed it, my apologise. I had read that procedure-pointers were not in yet, but the posts I had seen were old and I was hoping that maybe it had changed, likewise all of the ones I saw were in C so I was hoping (foolishly I’ll admit) that somehow FORTRAN did have them. Thank you very much for your assistance.