Forall loop confusion

Accelerated Computing HPC Compilers Legacy PGI Compilers
1

The following function describes a map between two vectors of dimension n to a vector of dimension n^2. Here, vector is a type containing a rank-1 pointer, “values”, and a logical scalar, “tmp”.

module c_algebra

  type vector
    private
    complex(dp), dimension(:), pointer :: values => null()
    logical                            :: tmp    =  .false.
  end type vector

  function vec_dir_vec(v1,v2 ) result(v3)
    implicit none
    type(vector), intent(in)  :: v1,v2
    type(vector), intent(out) :: v3
    integer                   :: n1, n2, i, j
    
    n1=size(v1); n2=size(v2)
    allocate(v3%values(n1 * n2))
    
!    forall ( i = 1:n1, j = 1:n2 ) 
    do i=1, n1
       do j=1, n2
          v3%values(2*(i-1)+j) = v1%values(i) * v2%values(j)
       end do
    end do
!    end forall

    v3%tmp=.true.
    call cleanup( v1 , .true. )
    call cleanup( v2 , .true. )
  end function vec_dir_vec
end module c_algebra

Compiling and running with portland 6.1-5 on a x86-64 machine using RHEL produces the following results for v1%values=(/cmplx(1,0),cmplx(2,0/) and v2%values=(/(3,0),(4,0)/). With a forall loop, v3%values=((2,0),(4,0),(2,0),(4,0)), whereas with the nested do loops, the correct result of ((1,0),(2,0),(2,0),(4,0)) is produced.

Now I’m pretty sure the two should be identical. Can anybody help me here?

2

Further, should the following work, where a1 and a2 are integer arrays, larray is a logical array (all size n), and flag is a logical scalar?

larray=a1.eq.a2
    flag=.true.
    forall ( i=1:n )
       if (.not.larray(i)) then
          flag=.false.
          return
       end if
    end forall

Or is there perhaps an inbuilt function to discover whether arrays are identical?

3

Your thinking of

forall( i = 1:n1, j = 1:n2) 
   v3(i, j) = expression
end forall

mapping to:

do i = 1:n1
   do j = 1:n2
      v3(i,j) = expression
    end do
end do

All right hand side expressions get evaluated before the assignments, and may be calculated out of order or in parallel. The left hand side array indices are also calculated in any order, and must not overlap. Many-one assignments (as your code has) are not allowed.

The bottom line is that you need to use the do loops here.

Hope this helps,
Mat

4

Just to clear up any confusion. There was a typo here. The forall loop was supposed to be

forall ( i = 1:n1; j = 1:n2 )
   v3%values(n2*(i-1)+j) = v1%values(i) * v2%values(j)
end forall

… which produces the same spurious answer.

As you can see, this is actually one-to-one, in the sense that any pair of indices (i,j) in (1:n1, 1:n2) will produce a unique index on the LHS. Whilst I can understand that it may be difficult for the compiler to distinguish this from a many-one assignment, it is not such, so I think the standard rather expects it to work. Certainly had I not tested this carefully, I’d not have spotted it. At the very least, should not the compiler flag such things? I’m pretty sure the code is legal, and it’s a problem if it simply doesn’t work.

5

Hi Che,

Thanks for the clarification. I’m able to reproduce the error now and agree that it’s compiler problem. Note that it seems to only occur when a derived type’s allocatable array is used on the left-hand side of the assignment. I have filed a technical problem report (TPR#4045) and will have our compiler team investigate.

  • Mat
6

Many thanks indeed.