给定一个巨大的对称正定矩阵,如何计算其逆的几个对角元素？

Question

更新:这是一个纯粹的Fortran问题 ; 我把数学的东西放在M.SE上.

考虑Px P对称和正定矩阵A(P = 70000,即A使用8字节双精度大约为40 GB).我们想要计算逆矩阵的前三个对角元素inv(A)[1,1],inv(A)[2,2]和inv(A)[3,3].

我发现 James R. Bunch的这篇论文似乎解决了这个问题而没有计算完全逆inv(A); 不幸的是他使用了Fortran和LINPACK,这两者都是我从未使用过的.

我正在尝试理解这个功能:

    SUBROUTINE SOLVEJ(A,LDA,P,Y,J)
    INTEGER LDA,P,J
    REAL A(LDA,1),Y(1)
C
    INTEGER K
    Y(J) = 1/A(J,J)
    DO 10 K = J + 1,P
    Y(K) = - SDOT(K - J,A(J,K),1,Y(J),1)/A(K,K)
    10 CONTINUE
    RETURN
    END

其中A是大小为LDA x P的矩阵,Y是一个长度为P的向量.

你能解释为什么他Y(1)在函数头中定义然后分配给Y(J)？Fortran不关心已定义数组的大小并允许您访问超出其结束的范围吗？为什么不定义Y(P),这似乎可能根据这个Fortran入门？

Answer 1

首先，您应该了解不同的 Fortran 版本，尤其是 77 VS 90/95 及更高版本，而且您确实可以（通常）像在 C 中一样越界。Fortran 中的数组可能会引起很多混乱，我会说这有点乱。为了将讨论限制在您的特定情况下，我们可以使用以下事实：这是关于虚拟数组，它是出现在过程的虚拟参数列表中的数组。对于虚拟数组，我们可以有 3 种类型：

1. 显式形状：显式声明尺寸
2. 假定形状：没有给出尺寸，只有冒号来表示数组的等级
3. 假定大小：最后一个维度是星号，前导维度已明确声明

使事情变得复杂的是，(3) 可以与(1) 分组，并且(2) 通常与延迟形状数组分组，例如可分配数组。deferred-shape 和假定形状仅适用于 Fortran 90/95 及更高版本，如果您想将它们用作虚拟参数，则需要显式接口，因此它通常在模块中使用。

因此，就您的情况而言，虽然Y(1)可以工作，因为您可能会超出范围，但这是非常糟糕的，因为当您使用-fcheck=bounds. 应编写有效的 Fortran 77 之一：

REAL A(LDA,*),Y(*)

或者，更好：

REAL A(LDA,P),Y(P)