forall index variables

In FORALL assignment, using a forall statement, or a forall construct, the index variable or variables should retain the defined/undefined status (and value if defined) they have at the start of the statement when the forall assignment is completed.

So in the following code, if i is undefined at the start it is undefined at the end, and if i = 4 (say) at the start then i = 4 at the end.

integer i
...

forall(i=1:10)
   x(i) = ...
end forall

This seems to be treated correctly. Essentially the index variables have a different scope to the local ones. This 'rule' means the following code, in which the forall is nested inside a DO loop with the same index, is perfectly legal.

program main

   integer i
   real b(5)
   real :: a(5) = (/1.0, 2.0, 3.0, 4.0, 5.0/)

   b = 0.0

   ! Should be able to nest forall inside do loop and use same index variable
   do i=1,2
     
      ! forall doesn't change defined/undefined status (or value if defined) of i
      forall(i=1:5)
         b(i) = b(i) + a(i)
      end forall
      
   end do
   
   ! should print 3
   print *, i

   ! should print 2.0, 4.0, 6.0, ...
   print *, b

end program main

However, the compiler generates **error 953 - i is already being used as a do (or implied do) index **. If I make the compiler ignore the error using /IGNORE 953, the code compiles and prints the correct value for i of 3 (it works properly even with /CHECKMATE turned on). This means the compiler will do the correct thing if the error message can be subverted.

Can forall indices be eliminated from the check for error 953? I realise this is not an important bug, just something you might see how to address easily at some time in the future.

Regards, David.

David,

Where do you go to learn these things? I've been reading books on post-Fortran77 all the way since Metcalf & Reid's 'Fortran 8x Explained' (which was what it was called until the 80s decade was well and truly over) - and never understood the half of it. Mind you, I struggled for a long time with FORMATs when converting from Algol-60 back in 1969, so wrestling with things written in Klingon is nothing new!

The error message you refer to is near to what a Fortran oldie would expect - so perhaps it merits a Warning at least.

It always seemed to me that the standards in Fortran 90 et seq were specified by people who hated Fortran as it was. A huge amount of effort was expended in doing things a different way, rather than (say) providing a standard method of accessing graphics, or for interacting with users - things that are done well by Clearwin, and which are far more useful than multiple ways of defining the precision of a variable, or of doing DO loops.

Eddie

Hi Eddie.

Well I'm one of those people who like to study the different standards and then try things out on different compilers. For a good translation of the standard I use 'The Fortran 2003 Handbook' but I also use one of Metcalf and Reids books plus the different documentation with the compilers. It helps that I have a good understanding of C++ as well, which helps to make a lot of the 'concepts' a bit less alien. (I started out writing Fortran IV then Fortran 77 and Pascal in the 80s).

Forall was added to Fortran 95 to provide opportunities for a programmer to hint at vectorisation on a compiler/processor that supports it (using SSE3 for example). Silverfrost FTN95 processes a Forall sequentially as far as I know, so it looks like a loop. However, technically it isn't a loop, but a fancy assignment.

A lot of the new things added from Fortran 90 can seem a bit unecessary for the Fortran 77 programmer. However, generally, the features provide ways of doing things that are essential for writing correct Fortran programs of any complexity (e.g. Explicit interfaces and Modules).

To give a specific example, the introduction of Derived Types (and Classes in Fortran 2003) makes coding much easier. But it also means that array elements of larger types can be made to 'persist' on multiple, different cache lines (avoiding 'false sharing'), which can provide a boost to efficiency when processes which act on them are made to run on multiple processors (using OpenMP for example, or Fortran 2008 coarrays).

David

Thanks for the bug report. I have logged it for investigation.

David,

Thanks for the explanation - most helpful. So it is a type of assignment statement, not a loop. A pity it looks a lot like a loop then. I've spent a bit of time trying to find out about it, and what I discover is that it doesn't seem to be any sort of improvement on a DO loop.

One answer at least to my general ignorance: FORALL isn't in Metcalf and Reid's 'Fortran 8X Explained'.

My last remark in the previous posting stands. Do you think that you are the first person to discover the bug because you are the first person to use FORALL in FTN95? (Or the first person to reuse the 'i' variable in that particular way?).

Eddie

Eddie,

FORALL was only introduced in Fortran 95, so it wouldn't be in 'Fortran 8x explained' or any book about Fortran 90. It was originally introduced in High Performance Fortran (HPF) and was then included in Fortran 95.

The FORALL construct looks like a loop, but the FORALL statement does not.

The following code to extract the diagonal of a matrix is an example of the FORALL statement (meaning: assign the diagonal elements of A to the DIAG array):

FORALL(i=1:10) DIAG(i) = A(i,i)

This cannot be done using array notation, DIAG(:)=A(:,:) since the left-hand side and right-hand side are not conformable (they have different shapes)

Another example of a FORALL statement is the following code which copies the diagonal to another array in reverse order:

FORALL(i=1:10) B(i) = DIAG(11-i)

This cannot be done using array notation, B(:) = DIAG(:), since there is not way to do the reversal.

The FORALL construct is just a way of simplifying two or more FORALL statements which have the same indices. The above two statements can therefore be combined into a FORALL construct:

FORALL(i=1:10) DIAG(i) = A(i,i)
FORALL(i=1:10) B(i) = DIAG(11-i)

is the same as

FORALL(i=1:10)
   DIAG(i) = A(i,i)
   B(i) = DIAG(11-i)
END FORALL

It is important to note that first line inside the construct sets up ALL of the values of DIAG (for i=1:10), then the second line copies the values of DIAG to B in the reverse order.

This is definitely [u:c087d2b317]not[/u:c087d2b317] the same as the following DO loop, which is incorrect code, since DIAG(11-i) is not defined on the first 5 passes.

DO i=1, 10
   DIAG(i) = A(i,i)
   B(i) = DIAG(11-i)  ! Error, since DIAG(11-i) is undefined when i=1
END DO

Two DO loops are needed to do the same operations as the FORALL construct correctly:

DO i=1, 10
   DIAG(i) = A(i,i)
END DO
DO i=1, 10
   B(i) = DIAG(11-i) 
END DO

I know a lot of Fortran programmers who use FORALL a lot. The BUG I have found only came to light because I have been studying the scope of FORALL indices. I was looking specifically at how my different Fortran compilers (NAG, Intel, Silverfrost, Gfortran, Open64, Oracle, Absoft) handle these variables and checking for conformance with the standard.

David.

David,

An interesting and helpful exposition. Thank you for going to the trouble of explaining it. The two DO loops are the way I would have tackled this problem.

Eddie

David,

I note your confidence with FORALL, but I would be more cautious.

FORALL(i=1:10) DIAG(i) = A(i,i) B(i) = DIAG(11-i) END FORALL

With this construct, can you be sure that all compilers would process the first line completely before B(i) =... ?

Also, your statement: 'I know a lot of Fortran programmers who use FORALL a lot. ' Unfortunately, I know few Fortran programmers and fewer of them know FORALL !

Best wishes for the new year and look forward to your comments in 2013.

John

Quoted from JohnCampbell

With this construct, can you be sure that all compilers would process the first line completely before B(i) =... ?

All Fortran 95 compilers should, otherwise it would be a serious bug.

The point of FORALL is that it expresses the parallelism and lack of loop carried dependencies to the programmer. If the following appears in code I am writing or maintaining, I know that DIAG is assigned completely before B.

FORALL(i=1:10)
   DIAG(i) = A(i,i)
   B(i) = DIAG(11-i)
END FORALL

However, with the following code, I might be tempted to fuse the loops together to get more efficiency (thereby introducing a bug).

DO i=1,10
   DIAG(i) = A(i,i)
END DO
DO i=1,10
   B(i) = DIAG(11-i)
END DO

Hmm ... how can one be sure an optimizing compiler won't fuse these two loops together anyway. 😉

Best wishes to all for 2013.

David and John,

For a few years from early 1970 to about early 1974, if called upon to write code for this, I would probably have done the 2 do loops. After reading a slim tome by Kreitzberg & Schneiderman, I would have used just one loop to save the loop overhead:

      DO 100 I=1,10 
      DIAG (I) = A (I, I)
      B (I)    = A (11-I, 11-I) 
 100  CONTINUE 

Back in the 70s, I was using non-optimising compilers, and so would probably have precalculated 11-I, especially if I was using it a number of times rather than just twice. After a while, I might have indexed B differently, realising that A(I,I) would only be fetched once even if referred to twice:

C...    Assumes that A is fully populated first
      DO 100 I=1,10 
      DIAG (I) = A (I, I)
      B (11-I) = A (I, I) 
 100  CONTINUE 

On my first introduction to Fortran I was aghast at statement numbers: I was used to Algol-60's BEGIN ... END constructs, and used indenting to wander code across all 132 lineprinter columns. You got a sense of the flow of code in this way. I hated the 72 column business with its lack of space, and continuation codes in column 6. Eventually, I realised that a statement number is an 'outdent', and a DO loop with a statement number tells you where you are going before you get there! Moreover, since late 1983 I have had printers that only print 80 columns readably, and the lack of space became less critical, as I realised that I could only take in relatively simple formulae at a glance, or understand them decades (or even weeks!) after they were written. My Algol-60 codes had usually used much less than the 132 columns anyway by the time multi-level indenting was taken into account.

I'd go even further than John, and say that you, David, are the first Fortran programmer I have ever come across who even knew that FORALL existed, let alone used it and mastered it to the point of debating compile time errors (not that one comes across many Fortran programmers these days, particularly in the shed at the end of my garden).

In the early PC days, I collected compilers much as David does, but mainly used Microsoft's. It was something of a delight not to be stuck with just one compiler on any computer, as was the case in mainframe days (the CDC 6xxx series had MNF and FTN). The advantages of one were always offset by a disadvantage elsewhere. In the case of FTN95, having Clearwin is the trump card, but even in DBOS days, one had graphics, whole libraries of useful functions and primitive windowing - as well as the ability to access lots of RAM on a 386 or 486 (without which it wouldn't run anyway).

For this old dog, is FORALL really a new trick worth learning?

Eddie

Interesting points Eddie and John.

It's only worth it if you want or need the speed-up that vectorisation can bring and you are using a compiler that supports it. (4X double precision or 8X single precision, on the new processors which support AVX instructions.)

In the following code for in place LU factorisation, the FORALL version is about 8 times faster for moderate values of n (using the Intel compiler). Notice how most of the loops and overhead have disappeared from the FORALL version.

I doubt you would see much difference with FTN95 though.

subroutine lu_f77(a)
   integer i,j,k,n
   real a(n,n)
   do k = 1, n-1
      do i = k+1, n
         a(i, k) = a(i, k) / a(k, k)
      end do
      do i = k+1, n                  
         do j = k+1, n
            a(i, j) = a(i, j) - a(i, k)*a(k, j)
         end do
      end do
   end do
end subroutine lu_f77

subroutine lu_f95(a,n)
   integer i,j,k,n
   real a(n,n)
   do k = 1, n-1
      a(k+1:n, k) = a(k+1:n, k) / a(k, k)
      forall (i=k+1:n, j=k+1:n)
         a(i, j) = a(i, j) - a(i, k)*a(k, j)
      end forall
   end do      
end subroutine lu_f95

Since Fortran 95 must support both ways of doing things, I hope my bug report proves useful for Silverfrost to provide a compliant FORALL for the, possible minority, of programmers who utilise it (even if that minority is only me 😃).

David and Eddie,

An interesting twist to the tale of FORALL is that Salford FTN95 had some bugs with it's implementation of FORALL for a number of years. As a consequence I have been slow to use it's syntax. My impression is that FORALL was implemented in Fortran 95 as a pseudo 'parallel' code syntax, but few compilers provided any efficiency gains from it. Even today, other optimising compilers do not automatically parallelise FORALL structures, or show any performance benefits over DO codings for the same compiler options.

FORALL has the important restriction in that the code must be 'pure' and the order of calculation can not be assumed, which the compiler is to use when optimising.

Most FORALL can be easily replaced by a 1-line DO loop, such as: DO I=1,10 ; DIAG (I) = A (I, I) ; END DO

It can be an elegent syntax, but provides little practical benefit over DO.

John

Sometimes FORALL is clearer code.

Look at these ways of shuffling the contents of an array. The FORALL is much clearer and closer to the mathematical syntax; the DO loop needs to work sequentially in reverse order.

program shuffle

   real a(5)

   ! Shuffle with forall   
   a = (/(real(i), i=1,5)/)
   
   forall(i=1:4) a(i+1) = a(i)
      
   print *, a

   ! Incorrect implementation of shuffle with do
   a = (/(real(i), i=1,5)/)

   do i=1,4
      a(i+1) = a(i)
   end do
   
   print *, a

   ! Correct shuffle with do (need to go backwards!!!)
   a = (/(real(i), i=1,5)/)
   
   do i=4,1,-1
      a(i+1) = a(i)
   end do
   
   print *, a

end program shuffle

Several points from this discussion, firstly that a compiler ought to conform to the standard, even if it seems to anyone that the particular facility in a standard is pointless. David's findings are therefore of great use.

Secondly, the standard doesn't say anything about 'warnings', and these can be as helpful as possible to the majority of users (FTN95 scores well here). But, the 'errors' should be 'standard conforming'.

Thirdly, FORALL doesn't actually do anything that you can't do in a DO loop - all these discussions end up by showing how the DO loop would look like. Sometimes, the DO loop(s) are a bit contrived. A useful comment that the loop needs to be 'PURE'. I always believed that FUNCTIONs should be 'pure', but some programming styles lead to 'impure' functions, and I suspect that the changes in scoping rules with MODULEs makes it more likely that routines will not be 'PURE'. I would concede that there will be cases where 90 and 95 constructs might be more elegant than 77 constructs, but I have yet to encounter a situation in my own programming that needs them.

As well as missing the opportunity to add graphics and GUI facilities to Fortran, the Standard Committee could have made in incumbent on compiler writers to generate the most efficient code ...

Eddie

Eddie,

Your points are well made. However F77 had its problems. F90 standard was a significant improvement, although I regret that I did not provide input into the standard committee when it was being developed. The standard committee has lacked the input from Fortran users, being more dominated by compiler developers and computer scientists who use other languages. F77 had a number of poor compromises and lack of definition, due to operating system limitations of (a) major hardware manufacturer. There was discussion in 1980's that a significant source of coding errors could be addressed by improving the linking process, but like GUI, none of this became included in the standard. The complexity of KIND is typical of how to make something simple more complex. It implies you can 'adjust' the precision, when there are very few hardware options available. REAL*8 is much more effective for a programmer, but never got to the main table. Post F95, Fortran has been taken over by C programmers, with most of the effort to turn Fortran into C, but little attention to improving the needs of numerical calculations. Improved functionality for ALLOCATE has not happened, nor has addressing GUI's or the linking process. We now have inter-operability and more obscure allocated derived types, few of which are needed for FORmula TRANslation computations. I wrote indexed list structures back in 70's and do not see any great leap forward with F2008. All I see is a standard that is so complex, that few compiler developers would be able to pay for developing a stable conforming compiler. While I might dismiss your holding on to F77, I am exactly the same as you as I hold onto F95.

Unfortunately, the code committee has long dismissed us programmers and our comments and persists with a new Fortran version of C that fewer require. Any improvements are coming from compiler developers and hardware, such as Checkmate from FTN95 and OpenMP for parallel computation, outside the standard. F90 was developed with the intention of improving portability, which is an aim that is being forgotten.

My apologies David, if these old views differ from your experience, but we must all embrace the diversity.

Thanks for your discussions in 2012 and I look forward to your comments and assistance in 2013,

John

David, being an oldie, I hardly ever if at all use a forall statement. However, I have just created a slightly different version of your program where the outer do loop variable is used as an array index.

program main 

   integer i 
   real b(5,2) 
   real :: a(5) = (/1.0, 2.0, 3.0, 4.0, 5.0/) 

   b = 0.0 

   ! Should be able to nest forall inside do loop and use same index variable 
   do j=1,2 
      
      ! forall doesn't change defined/undefined status (or value if defined) of i 
      forall(i=1:5) 
         b(i,j) = b(i,j) + a(i) 
      end forall 
      
   end do 
    
   ! should print 3 
   print *, i 

   ! should print 2.0, 4.0, 6.0, ... 
   print *, b 

end program main 

Note this is hypothetical and not intended to give the same result as yours. I have changed the do loop variable to j to differentiate it from the i of the forall construct, but I think that this is a perfectly valid program structure that may be required. If the do loop index remained as i, then there would be a conflict just as in any other nested do loops.

The error provided by FTN95 is therefore perfectly reasonable, but feel free to shout at me if you wish.

Ian

Ian,

You are an oldie, just like me. The critical thing (pointed out by David) is that FORALL is not a loop construct, it is a parallel assignment. Hence, it is like the array assignment:

      DIMENSION A(5), B(5), C(5)
      A = B + C

Now, you or I could write that in a loop:

     DO 10 I=1,5
     A(I) = B(I) + C(I)
  10   CONTINUE

and we would be happy to do so. We would also know that at the end of the loop, I would be undefined, and on the basis of experience and a number of compilers, we would know that I was perhaps 5, or more probably 0, or perhaps even -1. Whatever, we would know that the loop messed around with the value of I. What we would be absolutely sure of is that I would not contain any value it had before the loop was entered, except by chance.

The rules are different for FORALL, and different to any Fortran construct you ever came across up to Fortran-77. The I in a FORALL is not the same as a variable of the same name elsewhere in the (sub)program unit: it is a new variable whose scope is local to the FORALL. Essentially, FORALL is a potential subset operation on an array assignment, so that one might have (for example) a couple of FORALLs so that even and odd I values might be B+C or B-C etc.

I will never use a FORALL, as the modification of the scoping rules upsets my view of what is correct and proper. David points out that the FORALL construct hints at the potential for parallelisation, but this is only possible if the FORALL is PURE (i.e. has no side effects). In the loop example above, that would be pure, but

     DO 10 I=1,5
     A(I) = A(I-1) + B(I) + C(I)
  10   CONTINUE

perhaps isn't (showing my ignorance here!) and if done using a FORALL instead of a DO, might conceivably give a different answer. It seems to me from the earlier discussion that a FORALL that used I-1 ought to use A(I-1) on entry to the FORALL, whereas a DO should use the A(I-1) from the previous loop step.

An implicit DO loop in a READ or WRITE statement follows what I consider to be classic scoping rules, and therefore has the effect I am used to. Mind you, Fortran 90 et seq. has interfered with the scoping rules through MODULE and CONTAINS to such an extent that all the old certainties have gone, unless one either embraces the new paradigm wholeheartedly or steadfastly refuses to have anything to do with it.

Eddie

Quoted from IanLambley

the error provided by FTN95 is therefore perfectly reasonable, but feel free to shout at me if you wish.

There is no error in the original program I posted though. 😄 It is not perfectly reasonable for a compiler that claims to be compatible with Fortran 95 to issue an error in this case.

Its a bug in the compiler and I hope it can be fixed in due course. There may be an underlying issue with FORALL indices scope as well, which I have made Paul aware of in another thread, and I am sure he will give it some attention, but I expect he has other fish to fry at the moment (64 bit stuff, etc.)

It doesn't help thinking about FORALL as a loop. As Eddie says above, it is really a more general form of array assignment. The indices are there to illustrate what gets assigned.

Its not even close to being a loop. If you run the following program you will see the forall and the do loop give completely different results.

program shuffle

   real a(5)

   a = (/(real(i), i=1,5)/)
   
   ! forall assignment

   forall(i=1:4)
      a(i+1) = a(i)
   end forall
     
   print *, a

   a = (/(real(i), i=1,5)/)

  ! do loop

   do i=1,4
      a(i+1) = a(i)
   end do
   
   print *, a

end program shuffle 

As I said, array assignment and FORALL is really only a hint to the compiler that something can be parallelised (i.e. SIMD parallelisation). An optimizing compilers can sometimes figure this out from a DO loop though, but the FORALL makes this easier.

There are other developments going on in Fortran 2008, such as DO CONCURRENT which really is a loop, and would be more suited to SPMD parallelisation.

In my codes (some work, some hobby stuff) I make extensive use of SPMD (using OpenMP) and SIMD parallelisation and I find I use FORALL quite a lot. I also find that codes I have to maintain also use (or misuse sometimes) FORALL, so others are using it and it is something I need to know about.

There's a lot of new things in Fortran these days. You don't have to use them, however. But a good compiler (which FTN95 definitely is) needs to support as many of the new things it can to encourage 'new adopters of Fortran' as well as keeping us 'oldies' happy.

I hope in a small way my reports here on 'corner cases' will help to improve FTN95 over time.

David,

I'm not the authority on what FORALL is or does - I rely on your knowledge here!

Having thought about it (a lot) I can maybe see a possible benefit of it. But the scoping rules for many of the things in Fortran-90 break what I consider to be sacrosanct rules, and seem to me to have the potential to cripple my ability to see certain types of bugs 'at a glance'.

Most of my programming life has been done with subset compilers, and I have no difficulty with the subset of FTN95 that I use being not even the whole of FTN77! Just so long as it has Clearwin+

Do you get the same behaviour with .NET?

Eddie

Good evening Eddie. How is the snow there 😃

When you use a Statement Function in Fortran 77 the arguments have a scope that is different to the scope of the containing function. Well, the use of a separate scope for forall is similar. What are the other Fortran 90 scoping rules you mentioned that 'break the rules'? Can you give an example?

When programming in Fortran I find the best approach is to stick with the parts of the language I know very well. FORALL happens to be one of those areas for me, but I accept its a personal choice.

I have not tried compiling under .NET. I don't use it much.