Fortran: array syntax

[Telemachus]

Well, I'm trying to move a bit into modern Fortran, and I am more used to Fortran77 syntax.

So, for example, when I wanted to make an array whose two of its components equal another array of another dimension, I used to do something like this:

     do k=1,N
      do j=1,N
       do i=1,N        
        f(i,j,k)=f(i,j)
       enddo
      enddo
     enddo

So that would do the job for what I need. My doubt is if with modern Fortran I could do it this way:

     do k=1,N
        f(:,:,k)=f(:,:)
     enddo

Thanks in advance.

 

[jedishrfu]

I think which ever way you do it if it works it will work equally fast. FORTRAN is pretty good at optimizing code.

My preference would be the first way as it it’s clear to anyone reading your code. One thing though you need different names for the f array like f and g or f1 and f2.

 

[Telemachus]

Right. Thanks.

 

[DrClaude]

An even better solution in modern Fortran would be

forall(k=1:N) f2(:,:,k) = f1(:,:)

That way, the compiler is not forced to do the assignment following a strict increment of k from 1 to N, allowing for further optimisation.

 

[Telemachus]

An even better solution in modern Fortran would be

forall(k=1:N) f2(:,:,k) = f1(:,:)

That way, the compiler is not forced to do the assignment following a strict increment of k from 1 to N, allowing for further optimisation.

Can a subroutine be called inside the forall statement?

 

[Mark44]

Can a subroutine be called inside the forall statement?

Inside the body of the forall statement? Sure.

Later edit: No, subroutine calls are not allowed in a forall statement.

 

[Telemachus]

It gives me an error message, when I use the do loop I have no trouble with it.
program.f:817:72: Error: Unexpected CALL statement in FORALL block at (1)

I am using the gfortran compiler. BTW, the subroutine is using the index in the forall (k in here) when the call is made.

 

[Mark44]

It gives me an error message, when I use the do loop I have no trouble with it.
program.f:817:72: Error: Unexpected CALL statement in FORALL block at (1)

I am using the gfortran compiler. BTW, the subroutine is using the index in the forall (k in here) when the call is made.

Please show us the code that generates the error.

 

[Telemachus]

       forall(ja=1:N)
        call calcr(xc(ja),yc(ja),ja,rIp)
         fxp(:,:,ja)=fx(:,:)
         fyp(:,:,ja)=fy(:,:)
       end forall

If instead I use:

      do ja=1,2*mang
       call calcr(xc(ja),yc(ja),ja,rIp)
        fxp(:,:,ja)=fx(:,:)
        fyp(:,:,ja)=fy(:,:)
      enddo

I have no error messages and everything goes ok. I actually wanted to see if there was any performance improvement with the forall statement.

 

[Mark44]

I'm going to change my previous answer. As it turns out, forall is actually just an assignment-type statement, and is not a loop. Instead of using forall, you can use a do loop, as you already mentioned.

 

[Telemachus]

But should it work? or there is a problem with it?

 

[Mark44]

But should it work? or there is a problem with it?

From what I've read, it shouldn't work as you are trying to use it. It's really nothing more than a shortcut way to assign values to an array.

See https://en.wikipedia.org/wiki/Fortran_95_language_features#The_FORALL_Statement_and_Construct

BTW, if you get a compiler error on some code, your first instinct should be to check the documentation to see the correct usage.

 

[DrClaude]

      do ja=1,2*mang
       call calcr(xc(ja),yc(ja),ja,rIp)
        fxp(:,:,ja)=fx(:,:)
        fyp(:,:,ja)=fy(:,:)
      enddo

This is not efficient. I assume that there isn't a link between the function call and the arrays fx and fy. Therefore, you should have a loop for the call and a forall for the assignments.

 

[anorlunda]

But should it work? or there is a problem with it?

I had to upgrade my own FORTRAN knowledge in this case (FORTRAN IV and FORTRAN 77 were my experiences.)

Your error is thinking that the FORALL is a loop. That's not true. The Wikipedia article that @Mark44 linked says this.

When a DO construct is executed, each successive iteration is performed in order and one after the other—an impediment to optimization on a parallel processor.

FORALL(i = 1:n) a(i, i) = x(i)

where the individual assignments may be carried out in any order, and even simultaneously. The FORALL may be considered to be an array assignment expressed with the help of indices.

In other words, the compiler is free to implement the assignment in any way, not necessarily with a loop. It might have been less confusing if the feature was named ARRAY_ASSIGN rather than FORALL.

You can have arithmetic operations in the expression, but if you did have a function call within a FORALL (including, it might violate this.

Assignment in a FORALL is like an array assignment: as if all the expressions were evaluated in any order, held in temporary storage, then all the assignments performed in any order.

That is rather subtle, and it goes against the sequential execution of statements visualized by old fashioned FORTRAN programmers.

 

[Telemachus]

Great. Thanks. I think I can rearrange things to use forall. There was an slight improvement in the code after changing some do loops for this statement.

This was the running time after changing some do loops by forall:
2797.387 seconds.
And the time without those changes: 2964.344 seconds.

It is not a great improvement, but I think I can modify the code, using more arrays I should be able to include more forall's, at the expense of using more memory. But I would like to know how much faster could it be running. Also, maybe I can think some way to modify the subroutine call to include a forall in those situations.

This is not efficient. I assume that there isn't a link between the function call and the arrays fx and fy. Therefore, you should have a loop for the call and a forall for the assignments.

Sadly, there is a link between the call and those arrays. So I think I have to keep the do statement in that case.

One last question. The improvement comes from the execution in parallel of the assignments for the arrays? I haven't specified any parallel environment, so I think it shouldn't be that.

Thank you all.

 

[DrClaude]

One last question. The improvement comes from the execution in parallel of the assignments for the arrays? I haven't specified any parallel environment, so I think it shouldn't be that.

The speed up is not due to parallel execution, but due to memory management. The compiler is free to do the assignment any which way, so it can take advantage of the cache and avoid pagefaults. I guess that in some cases it can also reduce the overhead of the for loop.

 

[anorlunda]

One last question. The improvement comes from the execution in parallel of the assignments for the arrays? I haven't specified any parallel environment, so I think it shouldn't be that.

Where's your ego? You should expect that your program will be used by millions of people, over and over in all kinds of environments, for eternity.

Seriously, if it did everything strictly sequentially, then there would be no speed gain compared to a traditional DO loop. There are behind the curtain things that happen with memory access and cache memory that are completely invisible to you.

Edit: I see that @DrClaude beat me to it.

 

[Telemachus]

Great. Thank you. How do you know all of these things about how the compiler does the optimization? is there any reference, a book or something that explains all this?

 

[Mark44]

How do you know all of these things about how the compiler does the optimization?

By having a deeper level of understanding of the computer's architecture, including instruction pipeline, memory access vs. cache access vs. register access, multi-threading, to name a few. As for optimization, being able to see the code that the compiler actually emits makes it possible to optimize the sections of code that are taking the most time.

is there any reference, a book or something that explains all this?

Michael Abrash published "Zen of Assembly Language," 1990, and "Zen of Code Optimization," 1994.
A more current resource that I use is Agner Fog, a Danish computer scientist (https://agner.org/optimize/). The first manual listed on this page is C++ software optimization.

 

[harborsparrow]

Personally I don't recommend using a difficult-to-read (but shorter) syntax in code. When some future programmer looks at it, or you yourself 2 years later, it will be more helpful to see the longer but quite explicit version. So I like that first code you showed that you are trying to make more concise, and I don't like the more concise versions whatsoever because it is simply trickier to read.

 

[DrClaude]

Personally I don't recommend using a difficult-to-read (but shorter) syntax in code. When some future programmer looks at it, or you yourself 2 years later, it will be more helpful to see the longer but quite explicit version. So I like that first code you showed that you are trying to make more concise, and I don't like the more concise versions whatsoever because it is simply trickier to read.

We are not talking about obfuscating code here. We are talking about using the features of a high-level language to get simpler code that can be more easily optimised by a compiler. Come to think of it, it is actually the opposite of obfuscation. The idea is to code what you want the program to do, without having to detail exactly how to do it.