Given a huge symmetric positive definite matrix, how to calculate a few diagonal elements of its inverse?

Update: This is a pure Fortran question now; I put the maths stuff on M.SE.

Consider a PxP symmetric and positive definite matrix A (P=70000, i.e. A is roughly 40 GB using 8-byte doubles). We want to calculate the first three diagonal elements of the inverse matrix inv(A)[1,1], inv(A)[2,2] and inv(A)[3,3].

I have found this paper by James R. Bunch who seems to solve this exact problem without calculating the full inverse inv(A); unfortunately he uses Fortran and LINPACK, both of which I've never used.

I'm trying to understand this function:

 开发者_运维知识库   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

where A is a matrix of size LDA x P and Y is a vector of length P.

Can you explain why he defines Y(1) in the function head but then assigns to Y(J)? Does Fortran just not care about the size of the defined array and lets you access beyond its end? Why not define Y(P), which seems possible according to this Fortran Primer?

First, you should be aware of the different Fortran versions, especially 77 VS 90/95 and beyond, and that indeed you can (normally) go out of bounds just like in C. Arrays in fortran can cause a lot of confusion, and I would say that it's a bit of a mess. To limit the discussion to your specific case, we can use the fact that this is about a dummy array, which is an array that appears in the dummy argument list of a procedure. For dummy arrays, we can have 3 types:

1. explicit shape: dimensions are explicitly declared
2. assumed-shape: no dimensions given, only colons to denote the rank of the array
3. assumed-size: last dimension is an asterisk, leading dimensions are explicitly declared

To complicate things, (3) can be grouped with (1), and (2) is usually grouped with deferred-shape arrays, such as e.g. allocatable arrays. The deferred-shape and assumed-shape is only for Fortran 90/95 and beyond and requires an explicit interface if you want to use them as dummy arguments, so it's typically used in a module.

So, in your case, while Y(1) works because you can go out of bounds, it's very bad since the program will fail when you would compile it with -fcheck=bounds. One should write either the valid Fortran 77:

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

or, much better:

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