algorithm - 公式中的智能求和)翻译)

Question

好吧，今天我带着这个疑问来到这里......

我想用 Fortran 写这个等式:

当然我可以采用“经典”的方式，这样写:

 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

但由于我是用 Fortran 编写的，所以我想用一种看起来“更像原始”的表达方式来编写它，因此，紧凑。我在想:

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

看起来更像原来的表达。但事实是，我很难找到一种方法来以简单紧凑的方式编写求和表达式。当然，我可以编写一个函数“real function summation(<expression>,<lower bound>, <upper bound> )”，但由于我们正在谈论 Fortran，我认为应该有一种简单的(也许是内在的(？))编写它的方式。

所以有一种紧凑的方式来编写该表达式，还是我必须采用更丑陋的方式(两个显式循环)？

编辑:在实际代码中 x是一个二维数组，每列有一个解。所以使用内部函数 sum到这里为止，这似乎是个好主意(如@alexander-vogt 在他的回答中所示)导致代码几乎相同的“紧凑性”:

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

Answer 1

关于:

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

(请注意，我更改了矩阵 A 的索引，Fortran 是主要列!)

因为我们重新使用 B 作为结果，所以这个循环必须按升序进行。我不确定是否可以使用 forall 来做到这一点(是否由编译器来选择如何继续？ - 请参阅 Fortran forall restrictions )。

将结果写入新向量C 不会覆盖B 并且可以按任何顺序执行:

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