Compiler generates call to wrong dotp() routine with /64

Mysterious crashes, slow convergence and incorrect results were observed when a large program was compiled with /64 and run. This behaviour occurred with /64, /64 /opt and /64 /check, but not with /64 /checkmate, nor with any 32 bit run. I am using FTN95 V 8.51

The following short subroutine demonstrates the problem. It is, of necessity, not ready to run, but compiling and observing the /exp listing shows the problem.

MODULE mcm
  IMPLICIT NONE
  INTEGER, PARAMETER :: kdp = SELECTED_REAL_KIND(14,60)
  REAL(KIND=kdp), DIMENSION(:), ALLOCATABLE, TARGET :: rhs
  INTEGER :: nbn
  REAL(KIND=kdp) :: epsslv
  INTEGER, PARAMETER :: lrcgd1=19
END MODULE mcm

SUBROUTINE gcgris(ra,xx)
  USE mcm
  IMPLICIT NONE
  REAL(KIND=kdp), DIMENSION(:,:), INTENT(INOUT) :: ra
  REAL(KIND=kdp), DIMENSION(:), INTENT(INOUT), TARGET :: xx
  !
  INTEGER :: i, j
  REAL(KIND=kdp) :: r00
!
!...... code to calculate array rhs(), removed for this reproducer
!
  r00 = SQRT(DOT_PRODUCT(rhs,rhs))
  xx = 0._kdp
  IF (r00 <= epsslv) RETURN
  DO  i=1,nbn
     DO  j=1,lrcgd1
        ra(j,i) = 0.0_kdp
     END DO
  END DO
END SUBROUTINE gcgris

If you compile this code with /64 /exp, you will see the following call in the code listing:

R:\lang\FTN95\tbed>grep -i dot gcg.lis
   0019     r00 = SQRT(DOT_PRODUCT(rhs,rhs))                                                     AT 1fb
0000025c(#78,117,10):      CALL      __dotp4

The routine being called is for evaluating the DOT_PRODUCT of two single precision arrays. The correct routine to call for double precision arrays is __dotp8.

The bug does not occur with /64 /checkmate or 32-bit compilations because in those cases the code for evaluating the scalar product is generated inline. Nor does the bug occur with /64 when the code contains the scalar product calculation of just local arrays.

I know this reply has little to do with the error, but I have been calculating RMS vibrations for the last week !
Could I suggest considering this alternative solution when using FTN95 /64. (I would encourage the awareness and use of DOT_PRODUCT8@. I also use DOT_PRODUCT8@ in my ddotp function.)

MODULE mcm 
  IMPLICIT NONE 
  INTEGER, PARAMETER :: kdp = SELECTED_REAL_KIND(14,60) 
  REAL(KIND=kdp), DIMENSION(:), ALLOCATABLE, TARGET :: rhs 
  INTEGER :: nbn 
  REAL(KIND=kdp) :: epsslv 
  INTEGER, PARAMETER :: lrcgd1=19 

  contains
    real (kind=kdp) function rms ( xx, n )
      REAL(KIND=kdp) :: XX(*)
      INTEGER*4 :: N
!
      real*8, external :: DOT_PRODUCT8@ 
      real(KIND=kdp) :: ss
      INTEGER*8 :: lx
      lx = n
      ss = DOT_PRODUCT8@ (xx, xx, lx ) 
      if ( ss > 0 ) ss = sqrt (ss)  !  rms should be ss = sqrt ( ss/n )
      rms = ss
    end function rms
  
END MODULE mcm 

SUBROUTINE gcgris(ra,xx) 
  USE mcm 
  IMPLICIT NONE 
  REAL(KIND=kdp), DIMENSION(:,:), INTENT(INOUT) :: ra 
  REAL(KIND=kdp), DIMENSION(:), INTENT(INOUT), TARGET :: xx 
  ! 
  INTEGER :: i, j, n 
  REAL(KIND=kdp) :: r00 
! 
!...... code to calculate array rhs(), removed for this reproducer 
! 
  n = size(rhs)          !   size(rhs) is of unknown kind
  r00 = rms ( rhs, n )   !   SQRT(DOT_PRODUCT(rhs,rhs)) 
  xx = 0._kdp 
  IF (r00 <= epsslv) RETURN 
  DO  i=1,nbn 
     DO  j=1,lrcgd1 
        ra(j,i) = 0.0_kdp 
     END DO 
  END DO 
END SUBROUTINE gcgris

Also, what is '0.0_kdp' ? Can't you just use '0' ?

I was able to reproduce the problem reported.

Thanks for the tip, John.

However, I think that the people developing the software may object to using DOT_PRODUCT8@ and other such functions and subroutines that are specific to a particular compiler. They would need to put directives in the source code and then run it through a preprocessor to use the standard DOT_PRODUCT for other compilers. Secondly, if they wish to prepare a single precision version of the EXE, they would again have to modify the source code, instead of changing the definition of KDP (Ken's Default Precision?) in just one place.

The choice of '@' as part of the name creates a bit of a problem since this character is not one of the characters approved in standard Fortran for use in a variable name -- '@' is a 'special character'.

I used /64 /exp /xref to see what is the type of RHS

In MCM it describes: DOUBLE PRECISION, POINTER, ALLOCATABLE, POINTER, TARGET, DIMENSION(:) :: RHS (from module MCM)

However in SUBROUTINE GCGRIS it describes: Pointer, POINTER, ALLOCATABLE, POINTER, TARGET, DIMENSION(:) :: RHS (from module MCM) and DOUBLE PRECISION :: FUNCTION DOT_PRODUCT

Not sure this suggests any indication of the problem, as 'DOUBLE PRECISION :: FUNCTION DOT_PRODUCT' suggests __dotp4 should not be used ? /64 /checkmate and 32 bit don't adopt __dotp8, but appear to replace by loop code. (I am using Ver 8.40.0)

What is the status of 'POINTER, TARGET' for type flexibility ? (after 30 years of F90, these are still not in my F95 subset !) Finally patched a run version, without POINTER, TARGET, but still uses __dotp4 and gets the wrong result. A definate bug.

Tried the following 2 changes, but still calls --dotp4: REAL*8, DIMENSION(:), ALLOCATABLE :: rhs r00 = DOT_PRODUCT(rhs,rhs)

finally changed to the following r00 = SUM(rhsrhs) which used real8 routine 000002a7(#197,141,43): CALL DSUM@

Definite bug of selecting __dotp4 rather than __dotp8 for DOT_PRODUCT

Regarding use of using DOT_PRODUCT8@; being contained in the same module as where kdp is defined does weaken your argument. Also a utility library for compiler specific routines is fairly common, although mine typically use F77 style wrappers, with no interface requirements.

The XREF listing seems to follow this pattern:

1. When it sees USE statements, the type and other attributes of variables in the named module are listed.

2. When the module variables are used in the body of the subprogram, the corresponding line numbers are listed.

That the attribute POINTER is listed three times instead of merging the three instances is a bit odd. Perhaps the code that generates the XREF listings should be cleaned up a bit.

Thus, we see module variables listed twice. The first time, the type information is given. The subsequent listing(s) of the same variable will not show the type, but will show other attributes.

R:\HST216\src>grep -i rhs gcgris.XRF
INTEGER(KIND=4) :: RHS_B:extent:1) (from module (MCM)
INTEGER(KIND=4) :: RHS_R:extent:1) (from module (MCM)
Pointer, POINTER, ALLOCATABLE, POINTER, TARGET, DIMENSION(:) :: RHS (from module MCM)
Pointer, POINTER, ALLOCATABLE, POINTER, TARGET, DIMENSION(:) :: RHSSBC (from module MCM)
Pointer, POINTER, POINTER, DIMENSION(:) :: RHS_B (from module MCM)
Pointer, POINTER, POINTER, DIMENSION(:) :: RHS_R (from module MCM)

Note that RHS_B and RHS_R have no type shown in their second appearances in the listing. Both are of type INTEGER, as shown in their first appearances.

As to the uses of POINTER and TARGET: version 2 of HST3D was released in 1998. At that time, many compilers did not allow allocatable arrays to be passed as actual arguments. FTN95 7.2 did not, either. Pointers gave some of the same functionality that allocatable arguments did, so they got used in quite a few packages of the time.

In HST3D, these pointers are used quite a bit to pass unit-stride segments of arrays as arguments. FTN95 8.51, 64 bit seems to generate rather inefficient code for such usage. The HST3D test problem ELDER_HEAT, for example, runs in about 2 seconds with gFortran and 32-bit FTN95, but with 64-bit FTN95 the run time jumps up to 117 seconds, 97 percent of that in EFACT.

P.S.: I find that the run time falls back from 117 s to about 2 seconds if I replace

     ipenvv => ipenv(ifirst:neqn+1)
     envutv => envut(iipenv:ipenv(neqn+1)) 
     envlv => envl(iipenv:ipenv(neqn+1)) 
     CALL el1slv(iband,ipenvv,envl,envutv)

by

     CALL el1slv(iband,ipenv(ifirst:),envl,envut(iipenv:))

and a similar call to ELSLV the same way. Thus, at least with the FTN95 compiler, one has to suffer a huge loss of performance in return for the improved readability of code that one gains by using pointers to represent array sections.

Quoted from mecej4 FTN95 7.2 did not, either. Pointers gave some of the same functionality that allocatable arguments did, so they got used in quite a few packages of the time.

I am not familiar with this problem. In older versions of FTN95 ( pre Ver7 ) I could use allocatable arrays as arguments, but had trouble using them with /check. (I mainly used with F77 style calls) Certainly, using allocatable arrays in modules was a lot easier. I have rarely used pointers, but allocatable arrays have always been an extremely useful part of FTN95

I think the use of 'Pointer,' in .xrf is to imply allocatable, ie without a pre-defined address ?

It is the attribute ALLOCATABLE of a **dummy ** argument that is the issue. If the callee is not going to change or query the allocation status of a dummy argument, then an actual argument, whether allocatable or not, can be passed in the usual way, and in the callee there is no need to declare the argument as allocatable.

The following modified code does use __dotp8, so there is something about the module and subroutine interface that is causing the problem. !MODULE mcm IMPLICIT NONE INTEGER, PARAMETER :: kdp = SELECTED_REAL_KIND(14,60) REAL(KIND=kdp), DIMENSION(:), ALLOCATABLE :: rhs ! REAL*8, DIMENSION(:), ALLOCATABLE :: rhs INTEGER :: nbn REAL(KIND=kdp) :: epsslv = 0.1d0 INTEGER, PARAMETER :: lrcgd1=19 !END MODULE mcm

!program run_test
  real(kind=kdp) :: ra(lrcgd1,8), xx(8)
  REAL(KIND=kdp) :: r00 
  INTEGER :: i, j 
! 
!...... code to calculate array rhs(), removed for this reproducer 
! 
  nbn = 8

  allocate ( rhs(nbn) )
  rhs = 1
!
  r00 = SQRT(DOT_PRODUCT(rhs,rhs)) 
!  r00 = DOT_PRODUCT(rhs,rhs)
!  r00 = SUM(rhs*rhs)
  xx = 0._kdp 
  write (*,*) 'r00=',r00
  IF (r00 <= epsslv) stop 
  DO  i=1,nbn 
     DO  j=1,lrcgd1 
        ra(j,i) = 0.0_kdp 
     END DO 
  END DO 
END  ! SUBROUTINE gcgris

Quite so, I tried to create a self-contained reproducer without a module, and failed.

mecej4

Many thanks for the bug report. This has now been fixed for the next release of FTN95.

Thanks for the rapid resolution of the problem.

The choice of '@' as part of the name creates a bit of a problem since this character is not one of the characters approved in standard Fortran for use in a variable name -- '@' is a 'special character'.

I suspect that was the point of it in the first place.

As to whether or not it is a problem, when you run a source code containing @ through another compiler, the lines that use it are highlighted.

Of course, if you want easy-to-use graphics, and an easy-to-use GUI builder system, then you have really no choice other than FTN95, in which case the use of @ is no problem at all.

In their wisdom, the Fortran committee have introduced over the years things that perhaps a lot of people asked for, but which I could easily do without, and failed to provide an alternative for Clearwin+, which is fundamental to what I need. That's where it's @!

Eddie

Paul,

Thanks again for the quick work. Can you indicate what restricted type of code caused the problem ?

John

I suspect that it was the TARGET attribute (on rhs) that caused the problem.

In the example code that I provided, the TARGET attribute is not needed. I removed it and recompiled. The incorrect call to __dotp4 was again generated.

OK sorry. It turns out to be related to the ALLOCATABLE attribute. FTN95 missed the extra indirection needed and got the type of the array wrong as a result.