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