Legacy Software and gfortran-10


I don't follow any Fortran blogs. In fact, I'm not sure if any of them exist. If there are such blogs though, I might expect them to be ablaze right now about gfortran version 10 and all the legacy software is breaking. What is the problem you ask? Well consider the following program, which wraps up computing the eigenvectors of a matrix using LAPACK.

SUBROUTINE Routine(MatA, MatV)
  DOUBLE PRECISION, DIMENSION(:,:), INTENT(IN) :: MatA
  DOUBLE PRECISION, DIMENSION(:,:), INTENT(INOUT) :: MatV
  CHARACTER, PARAMETER :: job = 'V', uplo = 'U'
  INTEGER :: N, LDA
  DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: W
  DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: WORK
  DOUBLE PRECISION :: WORKTEMP
  INTEGER :: LWORK
  INTEGER, DIMENSION(:), ALLOCATABLE :: IWORK
  INTEGER :: IWORKTEMP
  INTEGER :: LIWORK
  INTEGER :: INFO
  INTEGER :: II

  MatV = MatA

  N = SIZE(MatA,DIM=1)
  LDA = N

  !! Allocations
  ALLOCATE(W(N))

  !! Determine the scratch space size
  LWORK = -1
  CALL DSYEVD(JOB, UPLO, N, MatA, LDA, W, &
       & WORKTEMP, LWORK, &
       & IWORKTEMP, LIWORK, &
       & INFO)
  N = LDA
  LWORK = INT(WORKTEMP)
  ALLOCATE(WORK(LWORK))
  LIWORK = INT(IWORKTEMP)
  ALLOCATE(IWORK(LIWORK))

  !! Run Lapack For Real
  CALL DSYEVD(JOB, UPLO, N, MatV, LDA, W, &
       & WORK, LWORK, &
       & IWORK, LIWORK, &
       & INFO)

  !! Cleanup
  DEALLOCATE(W)
  DEALLOCATE(Work)
  DEALLOCATE(IWork)

END SUBROUTINE Routine

The important thing to note here is the two calls to DSYEVD. We perform two calls because the first one tells us the required scratch space size, and the second performs the actual decomposition. There is a subtle point about these calls. In the first call we pass WORKTEMP and IWORKTEMP which are double precision and integer variables. In the second, we pass the actual scratch space arrays. If we compile this code with gfortran version 9, there are no problems. But with version 10 we get the following errors:

   27 |        & WORKTEMP, LWORK, &
      |         2
......
   38 |        & WORK, LWORK, &
      |         1
Error: Rank mismatch between actual argument at (1) and actual argument at (2) (scalar and rank-1)

The solution is simple: change the variables WORKTEMP and IWORKTEMP to be arrays of length one, and now the arguments match. Another time I've seen this come up is in old code that does something like this:

CALL DCOPY(NX, ZERO, 0, ARRAY1(1:NX), 1)
...
CALL DCOPY(NX, ARRAY1, 1, ARRAY2, 1)

In the first case, we're using dcopy to initialize some memory to zero. In the second, we're using it as a memcopy. I'm not sure why other programmers didn't just use the assignment operator to initialize the memory, though I wouldn't be surprised if older compilers weren't always optimal.

One final place I've encountered these kinds of errors is in calls to MPI. Fortunately, placing the USE MPI statement in your code can resolve this (and should have been there in the first place to prevent bugs).

These errors can be resolved with a flag like --std=legacy, but this is for sure a good time to bring legacy code up to spec. Version 10 is already available on homebrew, so if you have Mac users now is the time to check your code!