Using elemental for functions taking parameterised derived types

[Andrew_Smith]

I have been attempting to convert some code to use parameterised derived types. My old code used include files to define the core of the functions and I was able to define multiple implementations of each function, each with a different length array or different REAL kind within the derived type.

I still need include files for the different kind types but the lengths are now generalised. I wrote functions like this:

module Bar
implicit none

integer, parameter :: SP = SELECTED_REAL_KIND(6)
integer, parameter :: DP = SELECTED_REAL_KIND(13)
integer, parameter :: QP = SELECTED_REAL_KIND(24)

type VarA(k, n)
   integer, kind :: k
   integer, len :: n
   real(k) :: v
   real(k) :: d(n)
end type

interface foo
   module procedure foo_SP, foo_DP, foo_QP
end interface

contains

   pure elemental function foo_SP(a, b) result(c)
      type(VarA(SP,*)), intent(in) :: a
      type(VarA(SP,a%n)), intent(in) :: b
      type(VarA(SP,a%n)) c
      include 'foo_comomon.f90'
   end function

   pure elemental function foo_DP(a, b) result(c)
      type(VarA(DP,*)), intent(in) :: a
      type(VarA(DP,a%n)), intent(in) :: b
      type(VarA(DP,a%n)) c
      include 'foo_comomon.f90'
   end function

   pure elemental function foo_QP(a, b) result(c)
      type(VarA(QP,*)), intent(in) :: a
      type(VarA(QP,a%n)), intent(in) :: b
      type(VarA(QP,a%n)) c
      include 'foo_comomon.f90'
   end function
   
end module

Compiling with XE 15.0.1.148 gives

error #7618: In an elemental subprogram, a dummy argument, or a subobject thereof, shall not appear in a specification-expr except as the argument to certain inquiry functions.

Is there a way to avoid this limitation within the standard otherwise I will need to write array variants of every method and I have hundreds to do ?

At the moment I don't see any advantage of using parameterised derived types

[Steven_L_Intel1]

The problem here is that elemental functions have additional restrictions so that the types and sizes of the dummy arguments are the same for all calls when expanded from an array reference. Because you say that a%n is a LEN parameter, this potentially makes B and the result C different sizes across the elemental calls. What you did with the KIND parameter is just fine. If these weren't elemental functions, it would be ok.

For a problem like this, I'd think that making d an allocatable with deferred-shape is more appropriate. 

[Andrew_Smith]

But I ensured B and C match the size of A by my declarations of B and C. An array of the types would all have the same size so elemental should be fine. Perhaps there is another way I can declare that the sizes of A B and C match? Unless the compiler knows this how will it handle the array assignments ? It seams something is missing from the standard.

I did try using allocatable but the drop in performance was horrible. It is quite a small object that is created and destroyed rapidly and the overhead of allocating and deallocating in a threaded environment is too much. Typical sizes of the array component "d" are 1 and 6.

I also suspect that vectorising is limited with allocatable.

[Steven_L_Intel1]

You didn't make the sizes match across all the calls to the scalar function if you have an array of these, since each array element might have a different value for N.

You can use a KIND parameter to declare lengths, this has to be a compile-time constant so you can't make it A%N.

[Andrew_Smith]

How would you declare and array of VarA where n differed in elements of the array ?

[JVanB]

It's not a question of what can happen in a given case. It might require a very complicated set of rules to pin down exactly when a sequence of specification expressions could potentially lead to the length-type parameter of the function result acquiring different values for different array elements. As it is, the presence of three kinds of expressions greatly complicates the language, and specification expressions in elemental functions are actually a fourth, minor, class of expressions. We sort of have to put up with this oversimplification of rules in other contexts, such as with generic procedure resolution, and there is another one where dummy arguments can't be invoked elementally, which last rule makes no sense to me whatsoever.

When everything that seems sensible is allowed, it can make life difficult for vendors. If you look at the f03/f08 scorecard (I can't remember the link), one of the last features not fully implemented in ifort is transformational intrinsics in initialization expressions, and in fact a recent test I posted showed some still missing, some available but not documented as such, and even an ICE and a runtime failure. Allowing transformational intrinsics in initialization expressions meant that the standards committee didn't have to answer questions about why they weren't allowed where they made sense, but forced vendors to implement them even when they weren't necessarily all that useful.

[IanH]

Restricting any specification expressions for the type parameters of a function result of an elemental procedure to constant expressions (which is the remnant constraint in F2008, the compiler is complaining about the more restrictive constraint that exists in F2003) also allows the type and type parameters of the function result to be evaluated, and hence the underlying storage units required for the function result to be determined, independently of an invocation of the elemental procedure, and hence independently of the number of elements that need to be evaluated for any elemental reference.

[JVanB]

Thanks for pointing this quirk out, Ian. Actually, what is allowed in f2008 is not a subset of what f2003 allows.

module M
   implicit none
   integer N
   contains
      elemental function F(x)
         character(LEN=N) F
         integer, intent(in) :: x
         F = REPEAT(ACHAR(MODULO(x,128)),N)
      end function F
end module M

program P
   use M
   implicit none
   integer x
   x = 64
   N = 5
   write(*,*) F(x)
end program P

This is allowed by f2003 rules, but not by f2008. gfortran with -std=f2008 and ifort with /stand:f08 both accept it, however.

[IanH]

That's F2003 as published, but that was "fixed" with corrigendum five.

(So that then breaks F95 compatibility?)

[FortranFan]

Repeat Offender wrote:

 ..  Actually, what is allowed in f2008 is not a subset of what f2003 allows.

 ...

Not sure if I understand what this means and even whether that is true; moreover it is unclear if the code snippet in Quote #8 is explicitly allowed by Fortran 2003.  Generally for ELEMENTAL procedures, Fortran 2003 is more restrictive and Fortran 2008 has removed these restrictions.

If one compares the Fortran standard documents, Fortran 2003 includes

2 C1279 In the scoping unit of an elemental subprogram, an object designator with a dummy argument
3 as the base object shall not appear in a specification-expr except as the argument to one of the
4 intrinsic functions BIT SIZE, KIND, LEN, or the numeric inquiry functions (13.5.6).

which is not present in Fortran 2008.

So it appears Andrew has to either wait until Intel Fortran becomes fully compliant with Fortran 2008 for his PDT example in the original post to work or he has to take an alternate approach.

Andrew,

Have you considered a derived type with type-bound procedures including a finalizer, generics, and some "state" management control for the allocation aspect (to prevent what you indicate a performance hit)?

[IanH]

Andrew's example includes use of a specification expression (a%n) for a length type parameter in a function result.  Fortran 2008 doesn't permit that - type parameters have to be specified by constant expressions. 

RO's point was that prior to F2003 corrigendum five, a type parameter for a function result was allowed to be a specification expression, so long as that expression didn't refer to a dummy argument (bar a few exceptions which included LEN - which is where I think things got problematic).  The requirement to be a constant expression is more restrictive than that.  Just removing LEN from the allowed exceptions would have been less severe.

I think something accidentally got thrown out with the bath water.

It is a shame that you cannot use a "direct" length type parameter of a dummy argument (versus a length type parameter of a subobject of a dummy argument - that wouldn't work - hence the change in F2003 corrigendum five) in the specifications for of the type parameters of a function result of an elemental function, but I can understand why the simplicity of specification that RO refers to might be chosen instead.

[Andrew_Smith]

I am open to using type bound procedures if it would help. An interesting compiler error occures then:

module Bar
   implicit none
   type VarA(k, n)
      integer, kind :: k
      integer, len :: n
      real(k) :: v
      real(k) :: d(n)
   contains
      procedure foo
   end type
   integer, parameter :: DP = selected_real_kind(13)
contains
   pure subroutine foo(b, a)
      type(VarA(DP,*)), intent(in) :: a
      type(VarA(DP,a%n)), intent(inout) :: b
      b%v = a%v
      b%d = a%d
   end subroutine
end module

 error #8739: All length type parameters of the passed object dummy argument must be assumed.  
 

But when I comply with that in the more useful example in my other post I get fortcom: Fatal: There has been an internal compiler error (C0000005).fortcom: Fatal: There has been an internal compiler error (C0000005).

It is a serious problem if we can not define multiple arguments have the same length. How does the compler handle assignment right now ? Try running this test of assignment. With or without elemental it crashes with array bound error !

module Bar
   implicit none
   integer, parameter :: DP = selected_real_kind(13)

   type VarA(k, n)
      integer, kind :: k
      integer, len :: n
      real(k) :: v
      real(k) :: d(n)
   end type
   
   interface assignment(=)
      module procedure ass_pair
   end interface
   
contains
   pure elemental subroutine ass_pair(b, a)
      type(VarA(DP,*)), intent(in) :: a
      type(VarA(DP,*)), intent(inout) :: b
      b%v = a%v
      b%d = a%d
   end subroutine
end module
   
program test
   use Bar
   type (VarA(DP,1)) :: v1, v2
   v1%v = 1.0_DP
   v1%d = 2.0_DP
   v2 = v1
end program