Smart summation in FOR(mula) TRAN(slation)

Question

Well, today I'm here with this doubt...

I want to write this equation in Fortran:

enter image description here

Of course I can take the "classic" approach and write it like this:

 do i=1,N
        ac=0.0
        do j=i+1,M
            ac=ac+A(i,j)*B(j)
        enddo
        B(i)=(B(i)-ac)/A(i,i)
  enddo

But since I'm writing in Fortran, I would like to write it with an expression that looks like "more like the original" and, thus, compact. I was thinking in something like:

forall(i=1:N,j=i+1:M)
    B(i)=(B(i)- <MySummationExpression(A(i,j)*B(j))> )/A(i,i)
endforall

Expression that looks much more like the original. But the thruth is that I'm having a hard time finding a way to write the summation expression in a simple and compact way. Of course that I can write a function "real function summation(<expression>,<lower bound>, <upper bound>)", but since we're talking about Fortran I thought that there should be a simple (maybe intrinsic(?)) way of writing it.

So there is a compact way to write that expression, or I have to take the uglier way (two explicit loops)?

EDIT: In the actual code x is a two dimensional array, with one solution in each column. So using the intrinsic function sum that, up to here, seems like a great idea (as @alexander-vogt shown in his answer) leads to almost the same "compacticity" of the code:

do j=1,size(B,2)
        do i=nA,1,-1
            B(i,j)=(B(i,j)-sum(A(i,i+1:nA)*B(i+1:nA,j)))/A(i,i)
        enddo
 enddo

A valid effort, but "smarter" written code is often less readable and may come with free bugs. Write code in the way that is easiest for you. — milancurcic, Oct 09 '13 at 22:27

score 2 · Accepted Answer · edited May 23 '17 at 12:30

2

What about:

do i=1,N
  B(i) = (B(i) - sum( A(i+1:M,i) * B(i+1:M) )) / A(i,i)
enddo

(Please notice that I changed the indexing of matrix A, Fortran is column major! )

As we are re-using B for the result, this loop has to be done in ascending order. I'm not sure it's possible to do this using forall (isn't it up to the compiler to choose how to proceed? - see Fortran forall restrictions).

Writing the result into a new vector C does not overwrite B and can be executed in any order:

forall ( i=1:N )
  C(i) = (B(i) - sum( A(i+1:M,i) * B(i+1:M) )) / A(i,i)
endforall

edited May 23 '17 at 12:30

Community

1
1

answered Oct 09 '13 at 21:35

Alexander Vogt

17,879
13
52
68

I see that you had changed the indexing, but in that way the result is wrong. I'll investigate further why. – Alfonso Santiago Oct 09 '13 at 22:21
2

Well, just changing to `A(i,i+1:M)` I get the same results as with the code in the original question... The matrix `A` has to be transposed to work with the solution I provided ;-) – Alexander Vogt Oct 09 '13 at 22:25

Smart summation in FOR(mula) TRAN(slation)

1 Answers1