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

Question

更新:这是一个纯粹的 Fortran 问题现在 ; I put the maths stuff on M.SE .

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

我找到了this paper 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) , 根据 this Fortran Primer 这似乎是可能的?

Answer 1

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

显式形状:显式声明尺寸

假定形状:没有给出尺寸，只有冒号表示数组的等级

假定大小:最后一个维度是星号，前导维度显式声明

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

所以，在你的情况下，虽然 Y(1)工作，因为你可以越界，这是非常糟糕的，因为当你用 -fcheck=bounds 编译它时程序会失败.应该编写有效的 Fortran 77:

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

或者，更好:

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