2

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

我想用 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
4

1 回答 1

2

关于什么：

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

（请注意，我更改了 matrix 的索引A，Fortran 是 column major！）

由于我们正在重复使用B结果，因此必须按升序执行此循环。我不确定是否可以使用forall（不是由编译器来选择如何进行吗？ - 请参阅Fortran forall 限制）。

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

forall ( i=1:N )
  C(i) = (B(i) - sum( A(i+1:M,i) * B(i+1:M) )) / A(i,i)
endforall
于 2013-10-09T21:35:43.457 回答