Manifesto of the Ninety-Three

franck_reinquin · ‎09-22-2023

Hello again,

Several days ago I posted a question about generic reallocation of arrays of polymorphic elements, which was answered by Steve Lionel. But when switching to the compilation options which are actually used in our software, I found another problem.

At first glance, this is not directly related to polymorphism. However, I am unable to create a smaller program that shows the same problems so that may have something to do with it.

Here is the same program that I submitted last time (with Steve's change), with one addition : I set the first three elements of the measurement array before resizing the array, and print the three elements after the resizing is done.

Depending on the compilation options, the result is OK or not.

module measurements_management

  ! Measurements types
  !-------------------
  type :: generic_measurement
  end type generic_measurement

  type, extends(generic_measurement) :: measurement_as_reals
    real     :: real_measurement
  end type measurement_as_reals

  type, extends(generic_measurement) :: measurement_as_integers
    integer  :: int_measurement
  end type measurement_as_integers


  ! Generic management of measurement arrays
  ! ----------------------------------------
  type measurements_buffer
    integer    :: nb_used = 0
    class(generic_measurement), allocatable, dimension(:) :: tab_measurements

  contains
    procedure   :: init    => buffer_init
    procedure   :: realloc => buffer_realloc
  end type measurements_buffer

contains

    !   create a buffer containing an array with a selected measurement type
    !   the measurement array has an initial dimension of 50 elements
    subroutine buffer_init(this, measurement)

      class(measurements_buffer), intent(inout) :: this
      class(generic_measurement) :: measurement

      allocate(this%tab_measurements(50), mold=measurement)
    end subroutine buffer_init

    !   for an existing buffer created for a specific measurement type,
    !   reallocate the measurements array, adding 100 empty elements and
    !   copying the existing elements
    subroutine buffer_realloc(this)

      class(measurements_buffer), intent(inout), target :: this

      !.. Local variables ..
      integer  :: current_size
      class(generic_measurement), allocatable, dimension(:), target :: tab_tmp

      current_size = size(this%tab_measurements)
      write(*,*) 'current size=', current_size
      tab_tmp = reshape(this%tab_measurements,[current_size+100],this%tab_measurements)
      write(*,*) 'tab_tmp size after copy = ',size(tab_tmp)

      call move_alloc(tab_tmp, this%tab_measurements)
	  
    end subroutine buffer_realloc

end module measurements_management



! example with a buffer of integer measurements
program polymorphic_realloc
  use measurements_management
  implicit none
  type(measurements_buffer), target :: buffer_of_integer_measurements
  type(measurement_as_integers) :: measurement_int
  class(generic_measurement), pointer :: ptr_measurement
  integer  :: i
  
  ! init buffer
  call buffer_of_integer_measurements%init(measurement_int)
  
  ! fill buffer
  do i=1,3
    ptr_measurement => buffer_of_integer_measurements%tab_measurements(i)
    select type (ptr_measurement)
       class is (measurement_as_integers)
         ptr_measurement%int_measurement = 10*i
    end select
  end do

  ! extend buffer
  call buffer_of_integer_measurements%realloc()

  do i=1,3
    ptr_measurement => buffer_of_integer_measurements%tab_measurements(i)
    select type (ptr_measurement)
       class is (measurement_as_integers)
         write (*,*) i, ' =>', ptr_measurement%int_measurement
    end select
  end do
end program polymorphic_realloc

which should display :

 current size=                    50
 tab_tmp size after copy =                    150
                     1  =>          10
                     2  =>          20
                     3  =>          30

This is ok with the following options :

ifort -i8

ifort -qmkl

ifort -i4 -qmkl

but not with :

ifort -i8 -qmkl

as shown below :

 current size=                    50
 tab_tmp size after copy =                    150
                     1  =>       139883383696650
                     2  =>                    20
                     3  =>       139883383696670

Unfortunately, this is the option combination that we are using. Contrary to the example above, we are using a lot of BLAS95 ans LAPACK95 routines in our real software.

I have also tried :

ifort -i8 -qmkl-ilp64

to force the ILP64 back-end of mkl, to no avail.

Strangely (or not) ifx -i8 -mkl gives the correct result.

Also, when I change the code to :

  type, extends(generic_measurement) :: measurement_as_integers
    integer*4  :: int_measurement
  end type measurement_as_integers

Then all combinations are ok. But this means changing a lot of things in our code.

Thanks for any help,

Franck

franck_reinquin · ‎09-26-2023

This smaller program shows the same behaviour.

module measurements_management

  ! Measurements types
  !-------------------
  type :: generic_measurement
  end type generic_measurement

  type, extends(generic_measurement) :: measurement_as_integers
    integer  :: int_measurement
  end type measurement_as_integers

contains

    subroutine array_realloc(tab_mes)

      class(generic_measurement), dimension(:), allocatable, intent(inout) :: tab_mes

      !.. Local variables ..
      integer  :: current_size
      class(generic_measurement), dimension(:), allocatable :: tab_tmp

      current_size = size(tab_mes)
      !allocate(tab_tmp(current_size+100), source=tab_mes(1))
      !tab_tmp(1:current_size) = tab_mes(1:current_size)
      tab_tmp = reshape(tab_mes,[current_size+100],tab_mes)
      call move_alloc(tab_tmp, tab_mes)
	  
    end subroutine array_realloc

end module measurements_management



program polymorphic_realloc
  use measurements_management
  implicit none
  type(measurement_as_integers) :: mes_int
  class(generic_measurement), dimension(:), allocatable, target :: tab_mes_int
  class(generic_measurement), pointer :: ptr_measurement
  integer  :: i
  
  ! init buffer
  allocate(tab_mes_int(50), mold=mes_int)
  
  ! fill buffer
  do i=1,3
    ptr_measurement => tab_mes_int(i)
    select type (ptr_measurement)
       class is (measurement_as_integers)
         ptr_measurement%int_measurement = 10*i
    end select
  end do

  ! extend buffer
  call array_realloc(tab_mes_int)

  write(*,*)
  do i=1,3
    ptr_measurement => tab_mes_int(i)
    select type (ptr_measurement)
       class is (measurement_as_integers)
         write (*,*) i, ' =>', ptr_measurement%int_measurement
    end select
  end do
end program polymorphic_realloc

It has to be compiled using both MKL support and 64-bit integers to show the problem. Please note that in this reduced example, MKL is not used but the mere presence of -qmkl changes the result.

With ifort v2021.3.0 : the two lines in 23-24 compile ok (-mkl -i8) and yield the correct result. The new line in 25 compiles ok but gives a wrong result (see below).

With ifort v2021.9.0 : the two lines in 23-24 are rejected. The new line in 25 compiles ok (-qmkl-ilp64 -i8) but gives a wrong result :

                     1  =>                    10
                     2  =>       139931968097812
                     3  =>                    30

If compiled with ifx, the program gives the correct result.

Lastly, forcing int_measurement to be a 32-bit integer (integer*4) does also correct the problem.

andrew_4619 · ‎09-26-2023

"the two lines in 23-24 are rejected." what is the error message?

franck_reinquin · ‎09-26-2023

Sorry, I forgot to insert it again :

error #8304: In an intrinsic assignment statement, variable shall not be a non-allocatable polymorphic.

This was the subject of my first question, and line 25 is the solution proposed by Steve Lionel. That solution works well but reveals something that looks like a bug with the compiler.

Franck

FortranFan · ‎09-26-2023

@franck_reinquin ,

Try the option per the code snippet below with various compiler options of interest to you such as `qmkl` (a side-bar here: please consider shying away from the compiler options `-i8` / `-i4`, etc.; consider using defined kinds of integer types instead):

module measurements_management

   ! Measurements types
   !-------------------
   type :: generic_measurement
   end type generic_measurement
   
   type, extends(generic_measurement) :: measurement_as_integers
     integer  :: int_measurement
   end type measurement_as_integers

contains

   subroutine array_realloc(tab_mes)

      class(generic_measurement), dimension(:), allocatable, intent(inout) :: tab_mes

      !.. Local variables ..
      integer :: current_size

      select type ( mes => tab_mes )
         type is ( measurement_as_integers )
            blk_extend: block
               type(measurement_as_integers), dimension(:), allocatable :: tab_tmp
               current_size = size( mes )
               allocate( tab_tmp(current_size+100) )
               tab_tmp(1:current_size) = mes(1:current_size) !<-- Take note of this assignment
               call move_alloc(tab_tmp, tab_mes)
            end block blk_extend
         ! type is ( .. ) ! extend for more types as necessary
         class default
            ! error handling
      end select
	  
    end subroutine array_realloc

end module measurements_management

So with this suggestion, take note of the assignment on line #27. If the type ("class") of interest i.e., `measurement_as_integers` involves derived type components such that the intrinsic assignment, which is basically a copy instruction, can be costly, then consider a "home-brew" with a procedure that enables performant `move` instruction instead of the assignment `=` operator.

Note the advice you received with `reshape` intrinsic is questionable. It can work but there is nothing in the standard that definitively guarantees you performance or accuracy when the optional `PAD` argument is in effect. It's basically a hit-or-miss depending on the compiler you use. I would not rely on it.

franck_reinquin · ‎09-28-2023

Thanks for your help,

Unfortunately, removing -i8 is not an option. Our project holds more than 1000 files, some of which are legacy code we don't wish to dive in.

The current code, even if not totally efficient, is generic and avoids having dedicated code for our 26 dervied types (the example above exhibits the problem we encounter but is very far from the real code).

As for the change that was proposed, I trust Steve Lionel. He has an incredible grasp on Fortran and how the compiler implements the standards.

So I believe that there is a bug with ifort. I may be wrong of course.

Franck

FortranFan · ‎09-29-2023

@franck_reinquin wrote:
..
Unfortunately, removing -i8 is not an option. Our project holds more than 1000 files, some of which are legacy code we don't wish to dive in.
The current code, even if not totally efficient, is generic and avoids having dedicated code for our 26 dervied types (the example above exhibits the problem we encounter but is very far from the real code).
As for the change that was proposed, I trust ..

@franck_reinquin ,

This codebase will be served well by viewing the modules such as `measurement_management` as independent library code that are authored in independent manner and the data such as measurement values are represented with clear and defined kinds rather than via compiler options. With a bit of thought the derived types can be written such as the compiler options don't impact them. But I will leave it that.

As to `RESHAPE` intrinsic, there is absolutely no defensible basis whatsoever to recommend it with arrays of derived types generally. And one takes into account all the possible ways in which compiler implementations will fail to process the instructions, especially when polymorphism is involved and as you notice with Intel Fortran already and likely with gfortran also, it makes it further questionable to advise its use. Your comment comes across as a corollary converse to the logical fallacy argumentum ab auctoritate . This is a sad state for Fortran. Any simple test one might attempt to verify the use of RESHAPE in such scenarios is likely to fail, either outright in terms of functionality and/or performance. You can use the following as one such case, you may help your codebase by reviewing the following closely:

module measurements_management

   ! Measurements types
   !-------------------
   type, abstract :: generic_measurement
   contains
      procedure(Imove), deferred :: move
   end type generic_measurement

   abstract interface
      elemental subroutine Imove( from, to )
         import :: generic_measurement
         ! Argument list
         class(generic_measurement), intent(inout) :: from  
         class(generic_measurement), intent(inout) :: to  
      end subroutine 
   end interface
   
contains

   pure subroutine array_realloc_reshape(tab_mes)

      class(generic_measurement), dimension(:), allocatable, intent(inout) :: tab_mes

      !.. Local variables ..
      integer :: current_size
      class(generic_measurement), dimension(:), allocatable :: tab_tmp

      current_size = size( tab_mes )
      tab_tmp = reshape( tab_mes, [current_size+100], tab_mes )
      call move_alloc( tab_tmp, tab_mes )
	  
   end subroutine array_realloc_reshape

   pure subroutine array_realloc_move(tab_mes)

      class(generic_measurement), allocatable, intent(inout) :: tab_mes(:)

      ! Local variables
      integer :: i
      integer :: current_size
      class(generic_measurement), allocatable :: tab_tmp(:)

      current_size = size( tab_mes )
      allocate( tab_tmp(current_size+100), mold=tab_mes(1) )
      do concurrent ( i = 1:current_size )
         call tab_mes(i)%move( tab_tmp(i) )
      end do
      call move_alloc( tab_tmp, tab_mes )
	  
    end subroutine array_realloc_move

end module measurements_management

module int_measurements_management

   use measurements_management
   
   type, extends(generic_measurement) :: measurement_as_integers
      integer, allocatable :: int_measurement(:)
   contains
      procedure :: move => move_int_measurement
   end type measurement_as_integers

contains

   elemental subroutine set_measurments( this, size_mes )

      class(generic_measurement), intent(inout) :: this
      integer, intent(in) :: size_mes

      integer :: i

      select type ( this )
         type is ( measurement_as_integers )
            this%int_measurement = [( 42 + i - 1, i = 1, size_mes )]
         class default 
      end select

      return

   end subroutine set_measurments

   elemental subroutine move_int_measurement( from, to )

      class(measurement_as_integers), intent(inout) :: from
      class(generic_measurement), intent(inout)     :: to

      select type ( to )
         type is ( measurement_as_integers )
            call move_alloc( from=from%int_measurement, to=to%int_measurement )
         class default
            ! error handling elided
      end select

      return

   end subroutine move_int_measurement

end module int_measurements_management

program p

   use, intrinsic :: iso_fortran_env, only : I8 => int64, WP => real64

   use measurements_management
   use int_measurements_management

   implicit none

   integer, parameter :: N = 1000
   integer, parameter :: DATA_SIZE = 1000000
   class(generic_measurement), allocatable :: mes(:)
   real(WP) :: start_time
   real(WP) :: end_time

   print *, "Initial array size, N = ", N
   print *, "Size of data in each array element, DATA_SIZE = ", DATA_SIZE

   blk1: block

      ! Set up foo array
      allocate( measurement_as_integers :: mes(N) )
      call set_measurments( mes, DATA_SIZE )

      print *, "Block 1: Use MOVE method:"
      print *, "Before reallocation:"
      print *, "size(mes) = ", size(mes)
      select type ( mes )
         type is ( measurement_as_integers )
            print *, "mes(N)%int_measurement(1) = ", mes(N)%int_measurement(1)
            print *, "mes(N)%int_measurement(DATA_SIZE) = ", mes(N)%int_measurement(DATA_SIZE)
         class default 
      end select

      call my_cpu_time( start_time )
      call array_realloc_move( mes )
      call my_cpu_time( end_time )

      print *
      print "(a,g12.4,a)", "MOVE method: CPU time = ", (end_time - start_time), &
         " seconds."
      print *, "After reallocation:"
      print *, "size(mes) = ", size(mes)
      select type ( mes )
         type is ( measurement_as_integers )
            print *, "mes(N)%int_measurement(1) = ", mes(N)%int_measurement(1)
            print *, "mes(N)%int_measurement(DATA_SIZE) = ", mes(N)%int_measurement(DATA_SIZE)
         class default 
      end select

      deallocate( mes )

   end block blk1
   print *
   blk2: block

      ! Set up foo array
      allocate( measurement_as_integers :: mes(N) )
      call set_measurments( mes, DATA_SIZE )

      print *, "Block 2: Use RESHAPE:"
      print *, "Before reallocation:"
      print *, "size(mes) = ", size(mes)
      select type ( mes )
         type is ( measurement_as_integers )
            print *, "mes(N)%int_measurement(1) = ", mes(N)%int_measurement(1)
            print *, "mes(N)%int_measurement(DATA_SIZE) = ", mes(N)%int_measurement(DATA_SIZE)
         class default 
      end select

      call my_cpu_time( start_time )
      call array_realloc_reshape( mes )
      call my_cpu_time( end_time )

      print *
      print "(a,g12.4,a)", "Use RESHAPE: CPU time = ", (end_time - start_time), &
         " seconds."
      print *, "After reallocation:"
      print *, "size(mes) = ", size(mes)
      select type ( mes )
         type is ( measurement_as_integers )
            print *, "mes(N)%int_measurement(1) = ", mes(N)%int_measurement(1)
            print *, "mes(N)%int_measurement(DATA_SIZE) = ", mes(N)%int_measurement(DATA_SIZE)
         class default 
      end select

      deallocate( mes )

   end block blk2

   stop

contains

   subroutine my_cpu_time( time )

      ! Argument list
      real(WP), intent(inout) :: time

      ! Local variables
      integer(I8) :: tick
      integer(I8) :: rate

      call system_clock (tick, rate)

      time = real(tick, kind=kind(time) ) / real(rate, kind=kind(time) )

      return

   end subroutine my_cpu_time

end program p

C:\temp>ifort /standard-semantics /free /heap-arrays /traceback p.f
Intel(R) Fortran Intel(R) 64 Compiler Classic for applications running on Intel(R) 64, Version 2021.10.0 Build 20230609_000000
Copyright (C) 1985-2023 Intel Corporation.  All rights reserved.

Microsoft (R) Incremental Linker Version 14.36.32537.0
Copyright (C) Microsoft Corporation.  All rights reserved.

-out:p.exe
-subsystem:console
-incremental:no
p.obj

C:\temp>p.exe
 Initial array size, N =  1000
 Size of data in each array element, DATA_SIZE =  1000000
 Block 1: Use MOVE method:
 Before reallocation:
 size(mes) =  1000
 mes(N)%int_measurement(1) =  42
 mes(N)%int_measurement(DATA_SIZE) =  1000041

MOVE method: CPU time =    0.000     seconds.
 After reallocation:
 size(mes) =  1100
 mes(N)%int_measurement(1) =  42
 mes(N)%int_measurement(DATA_SIZE) =  1000041

 Block 2: Use RESHAPE:
 Before reallocation:
 size(mes) =  1000
 mes(N)%int_measurement(1) =  42
 mes(N)%int_measurement(DATA_SIZE) =  1000041

Use RESHAPE: CPU time =   0.4040     seconds.
 After reallocation:
 size(mes) =  1100
forrtl: severe (157): Program Exception - access violation
Image              PC                Routine            Line        Source
p.exe              00007FF697552D36  MAIN__                    183  p.f
p.exe              00007FF69759FE9B  Unknown               Unknown  Unknown
p.exe              00007FF6975A0240  Unknown               Unknown  Unknown
KERNEL32.DLL       00007FF827517344  Unknown               Unknown  Unknown
ntdll.dll          00007FF8294026B1  Unknown               Unknown  Unknown

JohnNichols · ‎09-30-2023

@FortranFan , your remove to Latin brought a smile to my face.

The classic example is

One example of the use of the appeal to authority in science dates to 1923,[29] when leading American zoologist Theophilus Painter declared, based on poor data and conflicting observations he had made,[30][31] that humans had 24 pairs of chromosomes. From the 1920s until 1956,[32] scientists propagated this "fact" based on Painter's authority,[33][34][31] despite subsequent counts totaling the correct number of 23.[30][35] Even textbooks[30] with photos showing 23 pairs incorrectly declared the number to be 24[35] based on the authority of the then-consensus of 24 pairs.[36]

taken from Wikipedia, it is fair use of the copyright material.

The second is

Manifesto of the Ninety-Three

and finally, the Common law of England and the USA is built on the absurdity of the authority.

The simplest example is the definition of a watercourse, it has bed banks and water, ok in England from the 1200's but apply it to Australia or Arizona in the modern age and you have a problem. As I explained to a Supreme Court Justice in NSW once, we can take your concept and define it scientifically, he called me lost in Fairyland.

I contend that a watercourse does not have to have banks or water all the time, merely a bed and a potential water stream.

polymorphism, mkl and integer*8

Manifesto of the Ninety-Three