Pure subroutine error with do concurrent for external routine

I am getting this error:

for my main program

program test
  implicit none 

  integer ( kind = 4 ), parameter :: Nx = 64, Ny = 64
  real ( kind = 8 )   , dimension ( Nx, Ny ) :: phi, dfdphi
  integer ( kind = 4 )            :: i, j

  call random_number (phi)    

  do concurrent ( i = 1 : Nx , j = 1 : Ny )
     call Get_dphi ( phi, dfdphi, i, j )
  end do
  
end program test

The external routine is

  pure subroutine Get_dphi (  phi_, dfdphi_, i_, j_ )
    implicit none

    integer ( kind = 4 ), parameter                         :: Nx_ = 64,Ny_ = 64
    real ( kind = 8 ), dimension ( Nx_, Ny_ ), intent ( in )  :: phi_
    real ( kind = 8 ), dimension ( Nx_, Ny_ ), intent ( out ) :: dfdphi_
    integer ( kind = 4 ), intent ( in )                     :: i_, j_

    dfdphi_(i_,j_) =  2.0*phi_(i_,j_)*( 1.0 - phi_(i_,j_) )**2 &
         *( 1.0 - 2*phi_(i_,j_) ) 

  end subroutine Get_dphi

I do not see any error if it is used as an internal routine. So what makes it not pure here?

I think gfortran cannot see that the subroutine is pure because it is not in a module and an INTERFACE was not provided. Gfortran 13.0.0 20221218 compiles the following code, in which the subroutine is put in a module:

module m
contains
  pure subroutine Get_dphi (  phi_, dfdphi_, i_, j_ )
    implicit none

    integer ( kind = 4 ), parameter                         :: Nx_ = 64,Ny_ = 64
    real ( kind = 8 ), dimension ( Nx_, Ny_ ), intent ( in )  :: phi_
    real ( kind = 8 ), dimension ( Nx_, Ny_ ), intent ( out ) :: dfdphi_
    integer ( kind = 4 ), intent ( in )                     :: i_, j_

    dfdphi_(i_,j_) =  2.0*phi_(i_,j_)*( 1.0 - phi_(i_,j_) )**2 &
         *( 1.0 - 2*phi_(i_,j_) ) 

  end subroutine Get_dphi
end module m

program test
  use m
  implicit none 

  integer ( kind = 4 ), parameter :: Nx = 64, Ny = 64
  real ( kind = 8 )   , dimension ( Nx, Ny ) :: phi, dfdphi
  integer ( kind = 4 )            :: i, j

  call random_number (phi)    

  do concurrent ( i = 1 : Nx , j = 1 : Ny )
     call Get_dphi ( phi, dfdphi, i, j )
  end do
  
end program test
3 Likes

@Shahid ,

As illustrated by @Beliavsky , a reference to a PURE subprogram (e.g., in your main program) requires an explicit interface.

Try to avoid using implicit interfaces as much as you can. Consider using an explicit interface as imperative with all features introduced in the language starting Fortran 90; PURE procedures a good example of this.

By the way, you may want to consider moving away from “magic numbers” such as 8 for the KINDs of intrinsic types, instead try using named constants with language intrinsics such as SELECTED_REAL_KIND.

1 Like

Interface block usually serves that purpose, as to my understanding.

But using it still produces the same error!

program test
  implicit none 

  interface
     subroutine Get_dphi (  phi_, dfdphi_, i_, j_ )
       implicit none

       integer ( kind = 4 ), parameter                         :: Nx = 64,Ny = 64
       real ( kind = 8 ), dimension ( Nx, Ny ), intent ( in )  :: phi_
       real ( kind = 8 ), dimension ( Nx, Ny ), intent ( out ) :: dfdphi_
       integer ( kind = 4 ), intent ( in )                     :: i_, j_
     end subroutine Get_dphi
  end interface
  
  integer ( kind = 4 ), parameter :: Nx = 64, Ny = 64
  real ( kind = 8 )   , dimension ( Nx, Ny ) :: phi, dfdphi
  integer ( kind = 4 )            :: i, j

  call random_number (phi)    

  do concurrent ( i = 1 : Nx , j = 1 : Ny )
     call Get_dphi ( phi, dfdphi, i, j )
  end do

end program

The declared interface is not marked as PURE

1 Like