Intel® Fortran Compiler
Build applications that can scale for the future with optimized code designed for Intel® Xeon® and compatible processors.

What does #6541 mean?

JVanB
Valued Contributor II
1,094 Views

There was a question in stack overflow about implementing Matlab's diff function in Fortran. First attempt:

! diff.f90
module M
   use ISO_FORTRAN_ENV
   implicit none
   private
   public diff
   interface diff
      module procedure diff2
   end interface diff
   contains
recursive function diff2(A,n,dim) result(B)
   real(REAL64)  A(:,:)
   real(REAL64), allocatable :: B(:,:)
   integer, optional :: n, dim
   integer Ndim
   integer shB(rank(A))
   integer shA(rank(A))
   integer peB(rank(A))
   integer i

   if(present(n)) then
      if(n <= 0) then
!         allocate(B,source=A)
         B = A
         return
      else if(n > 1) then
!         allocate(B,source=diff(diff(A,n-1,dim),1,dim))
         B = diff(diff(A,n-1,dim),1,dim)
         return
      end if
   end if
   if(present(dim)) then
      if(dim <= 0 .OR. dim > rank(A)) then
!         allocate(B,source=A)
         B = A
         return
      else
         Ndim = dim
      end if
   else
      Ndim = 1
   end if
   shB = shape(A)
   shB(Ndim) = shB(Ndim)-1
   if(any(shB <= 0)) then
      allocate(B(shB(1),shB(2)))
      return
   end if
   peB = [(i,i=1,rank(A))]
   if(Ndim /= rank(A)) then
      peB([Ndim,rank(A)]) = peB([rank(A),Ndim])
   end if
   shA = shape(A)
!   allocate(B,source=reshape(reshape(eoshift(A,1,dim=Ndim)-A,shB1,order=peB),shB,order=peB))
   B = reshape(reshape(eoshift(A,1,dim=Ndim)-A,shA(peB),order=peB),shB,order=peB)
end function diff2
end module M
! main.f90
program P
   use M
   use ISO_FORTRAN_ENV
   implicit none
   real(REAL64), allocatable :: A(:,:), B(:,:)
   integer M1, N
   character(20) fmt
   integer i,j
   integer dim

   do dim = 1, rank(A)
      M1 = 3
      N = 4
      write(*,'(a)') 'A ='
      A = reshape([((10*i+j,i=1,M1),j=1,N)],[M1,N])
      write(fmt,'(*(g0))') '(',size(A,2),'(g0:1x))'
      write(*,fmt) transpose(A)
      B = diff(A,dim=dim)
      write(fmt,'(*(g0))') '(',size(B,2),'(g0:1x))'
      write(*,'(a)') 'B ='
      write(*,fmt) transpose(B)
   end do
end program P

Worked OK with gfortran, but with ifort I got:

diff.f90(55): error #6541: This element is not yet supported in this context.
[SHA]
   B = reshape(reshape(eoshift(A,1,dim=Ndim)-A,shA(peB),order=peB),shB,order=peB
)
-----------------------------------------------^
diff.f90(55): error #6361: An array-valued argument is required in this context.
   [RESHAPE]
   B = reshape(reshape(eoshift(A,1,dim=Ndim)-A,shA(peB),order=peB),shB,order=peB
)
---------------^
compilation aborted for diff.f90 (code 1)

And there is that error #6541, so what does it mean?

I figured the way to go was to create a rank-invariant expression for diff(A,1,dim=Ndim) so that it could be used in a function or subroutine with an assumed-rank input. Is this possible? Will it be possible when SELECT RANK is available? Actually a test program that incrementally wrote out a rank-invariant expression seemed to work OK, both on gfortran and ifort:

program Q
   use ISO_FORTRAN_ENV
   implicit none
   real(REAL64), allocatable :: A(:,:,:)
   integer M, N, P
   integer i, j, k
   character(30) fmt
   integer Ndim

   M = 3
   N = 4
   P = 5
do Ndim = 1, rank(A)
   allocate(A(M,N,P))
   write(fmt,'(a,i0,a)') '(',size(A,2),'(g0:1x))'
   A = reshape([(((100*i**2+10*j**2+k**2,i=1,M),j=1,N),k=1,P)],shape(A))
   write(*,'(a)') 'A ='
do i = 1, P
   write(*,fmt) transpose(A(:,:,i))
   write(*,'()')
end do
   write(*,'(*(g0:1x))') 'shape(A) =',shape(A)
   write(*,'(*(g0:1x))') '[(i,i=1,rank(A))] =',[(i,i=1,rank(A))]
   write(*,'(*(g0:1x))') 'eoshift([rank(A)-Ndim,(0,i=2,rank(A))],1-Ndim) =', &
      eoshift([rank(A)-Ndim,(0,i=2,rank(A))],1-Ndim)
   write(*,'(*(g0:1x))') 'eoshift([Ndim-rank(A),(0,i=2,rank(A))],1-rank(A)) =', &
      eoshift([Ndim-rank(A),(0,i=2,rank(A))],1-rank(A))
   write(*,'(*(g0:1x))') '[(i,i=1,rank(A))] + ' // &
      'eoshift([rank(A)-Ndim,(0,i=2,rank(A))],1-Ndim) + ' // &
      'eoshift([Ndim-rank(A),(0,i=2,rank(A))],1-rank(A)) =', &
      [(i,i=1,rank(A))] + &
      eoshift([rank(A)-Ndim,(0,i=2,rank(A))],1-Ndim) + &
      eoshift([Ndim-rank(A),(0,i=2,rank(A))],1-rank(A))
   write(fmt,'(a,i0,a)') '(',rank(A),'(g0:1x))'
   write(*,'(a)') 'spread([1,(0,i=2,rank(A))],2,rank(A)) ='
   write(*,fmt) transpose(spread([1,(0,i=2,rank(A))],2,rank(A)))
   write(*,'(a)') 'eoshift(' // &
      'spread([1,(0,i=2,rank(A))],2,rank(A)), 1-(' // &
      '[(i,i=1,rank(A))] +' // &
      'eoshift([rank(A)-Ndim,(0,i=2,rank(A))],1-Ndim) +' // &
      'eoshift([Ndim-rank(A),(0,i=2,rank(A))],1-rank(A)))) ='
   write(*,fmt) eoshift( &
      spread([1,(0,i=2,rank(A))],2,rank(A)), 1-( &
      [(i,i=1,rank(A))] + &
      eoshift([rank(A)-Ndim,(0,i=2,rank(A))],1-Ndim) + &
      eoshift([Ndim-rank(A),(0,i=2,rank(A))],1-rank(A))))
   write(*,'(a)') 'matmul(eoshift(' // &
      'spread([1,(0,i=2,rank(A))],2,rank(A)), 1-(' // &
      '[(i,i=1,rank(A))] +' // &
      'eoshift([rank(A)-Ndim,(0,i=2,rank(A))],1-Ndim) +' // &
      'eoshift([Ndim-rank(A),(0,i=2,rank(A))],1-rank(A)))),' // &
      'shape(A)) ='
   write(*,'(*(g0:1x))') matmul(eoshift( &
      spread([1,(0,i=2,rank(A))],2,rank(A)), 1-( &
      [(i,i=1,rank(A))] + &
      eoshift([rank(A)-Ndim,(0,i=2,rank(A))],1-Ndim) + &
      eoshift([Ndim-rank(A),(0,i=2,rank(A))],1-rank(A)))), &
      shape(A))
   write(fmt,'(*(g0))') '(',M,'(',N,'(g0:1x)/))'
   write(*,'(a)') 'eoshift(A,1,dim=Ndim)-A ='
   write(*,fmt) reshape(eoshift(A,1,dim=Ndim)-A,[N,M,P],order=[2,1,3])
   write(fmt,'("(",i0,"(",i0,"(g0:1x)/))")') matmul(reshape([1,0,0, &
                                                            0,1,0], &
      [2,3],order=[2,1]), &
      matmul(eoshift( &
      spread([1,(0,i=2,rank(A))],2,rank(A)), 1-( &
      [(i,i=1,rank(A))] + &
      eoshift([rank(A)-Ndim,(0,i=2,rank(A))],1-Ndim) + &
      eoshift([Ndim-rank(A),(0,i=2,rank(A))],1-rank(A)))), &
      shape(A)))
   write(*,'(a)') 'reshape(' // &
      'eoshift(A,1,dim=Ndim)-A,' // &
      'matmul(eoshift(' // &
      'spread([1,(0,i=2,rank(A))],2,rank(A)), 1-(' // &
      '[(i,i=1,rank(A))] +' // &
      'eoshift([rank(A)-Ndim,(0,i=2,rank(A))],1-Ndim) +' // &
      'eoshift([Ndim-rank(A),(0,i=2,rank(A))],1-rank(A)))),' // &
      'shape(A)),' // &
      'order =' // &
      '[(i,i=1,rank(A))] +' // &
      'eoshift([rank(A)-Ndim,(0,i=2,rank(A))],1-Ndim) +' // &
      'eoshift([Ndim-rank(A),(0,i=2,rank(A))],1-rank(A)))'
   write(*,fmt) reshape(reshape( &
      eoshift(A,1,dim=Ndim)-A, &
      matmul(eoshift( &
      spread([1,(0,i=2,rank(A))],2,rank(A)), 1-( &
      [(i,i=1,rank(A))] + &
      eoshift([rank(A)-Ndim,(0,i=2,rank(A))],1-Ndim) + &
      eoshift([Ndim-rank(A),(0,i=2,rank(A))],1-rank(A)))), &
      shape(A)), &
      order = &
      [(i,i=1,rank(A))] + &
      eoshift([rank(A)-Ndim,(0,i=2,rank(A))],1-Ndim) + &
      eoshift([Ndim-rank(A),(0,i=2,rank(A))],1-rank(A))), &
      matmul(eoshift( &
      spread([1,(0,i=2,rank(A))],2,rank(A)), 1-( &
      [(i,i=1,rank(A))] + &
      eoshift([rank(A)-Ndim,(0,i=2,rank(A))],1-Ndim) + &
      eoshift([Ndim-rank(A),(0,i=2,rank(A))],1-rank(A)))), &
      shape(A)), &
      order=[2,1,3])
   write(*,'(a)') 'shape(A)-eoshift([1,(0,i=2,rank(A))],1-Ndim) ='
   write(*,'(*(g0:1x))') shape(A)-eoshift([1,(0,i=2,rank(A))],1-Ndim)
   write(fmt,'(*(g0))') '(',M-merge(1,0,Ndim==1),'(',N-merge(1,0,Ndim==2),'(g0:1x)/))'
   write(*,'(a)') 'reshape(' // &
      'reshape(' // &
      'eoshift(A,1,dim=Ndim)-A,' // &
      'matmul(eoshift(' // &
      'spread([1,(0,i=2,rank(A))],2,rank(A)), 1-(' // &
      '[(i,i=1,rank(A))] +' // &
      'eoshift([rank(A)-Ndim,(0,i=2,rank(A))],1-Ndim) +' // &
      'eoshift([Ndim-rank(A),(0,i=2,rank(A))],1-rank(A)))),' // &
      'shape(A)),' // &
      'order =' // &
      '[(i,i=1,rank(A))] +' // &
      'eoshift([rank(A)-Ndim,(0,i=2,rank(A))],1-Ndim) +' // &
      'eoshift([Ndim-rank(A),(0,i=2,rank(A))],1-rank(A))),' // &
      'shape(A)-eoshift([1,(0,i=2,rank(A))],1-Ndim),' // &
      'order =' // &
      '[(i,i=1,rank(A))] +' // &
      'eoshift([rank(A)-Ndim,(0,i=2,rank(A))],1-Ndim) +' // &
      'eoshift([Ndim-rank(A),(0,i=2,rank(A))],1-rank(A)))'
   write(*,fmt) reshape(reshape( &
      reshape( &
      eoshift(A,1,dim=Ndim)-A, &
      matmul(eoshift( &
      spread([1,(0,i=2,rank(A))],2,rank(A)), 1-( &
      [(i,i=1,rank(A))] + &
      eoshift([rank(A)-Ndim,(0,i=2,rank(A))],1-Ndim) + &
      eoshift([Ndim-rank(A),(0,i=2,rank(A))],1-rank(A)))), &
      shape(A)), &
      order = &
      [(i,i=1,rank(A))] + &
      eoshift([rank(A)-Ndim,(0,i=2,rank(A))],1-Ndim) + &
      eoshift([Ndim-rank(A),(0,i=2,rank(A))],1-rank(A))), &
      shape(A)-eoshift([1,(0,i=2,rank(A))],1-Ndim), &
      order = &
      [(i,i=1,rank(A))] + &
      eoshift([rank(A)-Ndim,(0,i=2,rank(A))],1-Ndim) + &
      eoshift([Ndim-rank(A),(0,i=2,rank(A))],1-rank(A))), &
      [N-merge(1,0,Ndim==2),M-merge(1,0,Ndim==1),P-merge(1,0,Ndim==3)],order=[2,1,3])
   deallocate(A)
end do
end program Q

(ifort needs /assume:realloc_lhs for the above.) Encouraged by this result, I started on a version with a constant rank-invariant expression:

program P
   use ISO_FORTRAN_ENV
   implicit none
   integer, parameter :: M = 3
   integer, parameter :: N = 4
   real(REAL64) :: A(M,N)
   integer i, j
   parameter(A = reshape([((10*i**2+j**2,i=1,M),j=1,N)],shape(A)))
   integer, parameter :: Ndim1 = 1
   real(REAL64) :: B1(M-merge(1,0,Ndim1==1),N-merge(1,0,Ndim1==2))
   parameter(B1 = reshape( &
      reshape( &
      eoshift(A,1,dim=Ndim1)-A, &
      matmul(eoshift( &
      spread([1,(0,i=2,rank(A))],2,rank(A)), 1-( &
      [(i,i=1,rank(A))] + &
      eoshift([rank(A)-Ndim1,(0,i=2,rank(A))],1-Ndim1) + &
      eoshift([Ndim1-rank(A),(0,i=2,rank(A))],1-rank(A)))), &
      shape(A)), &
      order = &
      [(i,i=1,rank(A))] + &
      eoshift([rank(A)-Ndim1,(0,i=2,rank(A))],1-Ndim1) + &
      eoshift([Ndim1-rank(A),(0,i=2,rank(A))],1-rank(A))), &
      shape(A)-eoshift([1,(0,i=2,rank(A))],1-Ndim1), &
      order = &
      [(i,i=1,rank(A))] + &
      eoshift([rank(A)-Ndim1,(0,i=2,rank(A))],1-Ndim1) + &
      eoshift([Ndim1-rank(A),(0,i=2,rank(A))],1-rank(A))))
   integer, parameter :: Ndim2 = 2
   real(REAL64) :: B2(M-merge(1,0,Ndim2==1),N-merge(1,0,Ndim2==2))
   parameter(B2 = reshape( &
      reshape( &
      eoshift(A,1,dim=Ndim2)-A, &
      matmul(eoshift( &
      spread([1,(0,i=2,rank(A))],2,rank(A)), 1-( &
      [(i,i=1,rank(A))] + &
      eoshift([rank(A)-Ndim2,(0,i=2,rank(A))],1-Ndim2) + &
      eoshift([Ndim2-rank(A),(0,i=2,rank(A))],1-rank(A)))), &
      shape(A)), &
      order = &
      [(i,i=1,rank(A))] + &
      eoshift([rank(A)-Ndim2,(0,i=2,rank(A))],1-Ndim2) + &
      eoshift([Ndim2-rank(A),(0,i=2,rank(A))],1-rank(A))), &
      shape(A)-eoshift([1,(0,i=2,rank(A))],1-Ndim2), &
      order = &
      [(i,i=1,rank(A))] + &
      eoshift([rank(A)-Ndim2,(0,i=2,rank(A))],1-Ndim2) + &
      eoshift([Ndim2-rank(A),(0,i=2,rank(A))],1-rank(A))))
   character(30) fmt

   write(fmt,'(*(g0))') '(',size(A,2),'(g0:1x))'
   write(*,'(a)') 'A ='
   write(*,fmt) transpose(A)
end program P

But that one got a C0000005 with ifort.

 

0 Kudos
6 Replies
Steven_L_Intel1
Employee
1,094 Views

I have no idea what that means. If this were a specification or initialization expression, there might be some missing functionality but this is just an assignment. I will ask the developers.

0 Kudos
JVanB
Valued Contributor II
1,094 Views

Well it's sort of an initialization expression in that

20 SHAPE shall be a rank-one integer array. SIZE (x), where x is the actual argument corresponding to

21 SHAPE, shall be a constant expression whose value is positive and less than 16. It shall not have

22 an element whose value is negative.

0 Kudos
Steven_L_Intel1
Employee
1,094 Views

Hmm, I see. SHA(Peb) is not a constant expression but it is required to be. Maybe that's what the compiler tried to say. Still a rather unhelpful message.

0 Kudos
JVanB
Valued Contributor II
1,094 Views

Looking over the standard, shA(peB) is a variable in that shA is a variable and this is an array section of it with a vector subscript. Then the question is whether SIZE(shA(peB)) is a constant expression, and that flies because SIZE is an array inquiry function, so it's a specification inquiry and the result comes ultimately from SIZE(peB), which is given by RANK(A), and that is OK as a constant expression via the specification inquiry clause, all the way from F95 at least (if F95 had a RANK intrinsic).

We can confirm that SIZE(shA(peB)) is a constant expression by running it through our handy-dandy constant expression tester:

module M
   use ISO_FORTRAN_ENV
   implicit none
   contains
      recursive subroutine S(A)
         real(REAL64) A(:,:)
         integer peB(rank(A))
         integer shA(rank(A))
         type T
            integer :: x = 0
         end type T
         type(T) x(SIZE(shA(peB)))
         save
         print '(i0)', x(1)%x
         x(1)%x = 1
      end subroutine S
end module M

program P
   use M
   implicit none
   real(REAL64) A(2,2)
   A = 7
   call S(A)
   call S(A)
end program P

And both ifort and gfortran print out:

0
1

So it is a constant expression.

 

0 Kudos
Steven_L_Intel1
Employee
1,094 Views

The first case escalated as DPD200381581. The ICE is DPD200381582.

0 Kudos
Steven_L_Intel1
Employee
1,094 Views

First case (DPD200381581) is fixed for the 17.0 product release. The ICE is still in progress.

0 Kudos
Reply