Учитывая огромную симметричную положительную определённую матрицу, как вычислить несколько диагональных элементов её инверсии?

Обновление: Это чистый Фортран вопрос сейчас ; Я ставлю математические вещи на M.SE.

Рассмотрим PxP симметричную и положительную определенную матрицу A (P=70000, т.е. A составляет примерно 40 Гб с использованием 8-байтовых дубликатов). Мы хотим вычислить первые три диагональных элемента обратной матрицы inv(A)[1,1], inv(A)[2,2] и inv(A)[3,3].

Я нашел эту статью Джеймса Р. Банча, который, кажется, решил именно эту проблему, не вычислив полную обратную inv(A); , к сожалению, он использует Фортран и 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 Primer?