- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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.
Link Copied
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
The first case escalated as DPD200381581. The ICE is DPD200381582.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
First case (DPD200381581) is fixed for the 17.0 product release. The ICE is still in progress.

- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page