Finding diagonal elements in a Fortran array using C? Easy?

[susan_in_california]

Hello, all! I know this will be another easy one, but I just want to be absolutely sure I am doing the right thing. The LAPACK Cholesky factorization function (DPOTRF) documentation says the following:

"If UPLO = 'U', the leading N-by-N upper triangular part of A contains the upper triangular part of the matrix A, and the strictly lower triangular part of A is not referenced."

Now, I know that Fortran arrays are stored in column-major order (meaning the column entries are adjacent in memory), and I would like the diagonal elements accessing the array from my C++ code.

Therefore, clearly A[0] = A[0,0] (the first diagonal element).

However, if the lower triangular part of the matrix is not referenced, what about A[1,1]? Is A[(N + 1] = A[1,1] and A[N] = A[1,0] = 0? Or, since the lower triangular are not referenced, is A[N] = A[1,1]? And then A[2N - 1] = A[2,2]?

Basically, does "not referenced" mean "not present", and thus the entries are offset accordingly, preserving memory? Or does "not referenced" mean contains junk data, but is present for a place holder?

Thank you for any insight!

Best,
Susan

[Julien Langou]

Both of your answers are correct in the framework of LAPACK. All depends of the
storage format you choose for your coefficient matrix. For symmetric matrices, they are
basically two basic storage format: full storage format and packed storage format.

The full storage format corresponds to LAPACK routines that begins with DPO (for SPD)
or DSY (for symmetric). (For example DPOTRF.) The full storage format is the first one
you have described. You allocate a full matrix n-by-n (that is n^2 data value) and only
reference the upper part (in your case) (that is n(n+1)/2 date value) and yes LAPACK is
in Fortran, so that means that you need to deal with the column major by yourself. So
A[0] is A(0,0), then comes A(0,1) that is A[n], then comes A(1,1) that is A[n+1], then
comes A(0,2) that is A[2n], then comes A(1,2) that is A[2n+1], then comes A(2,2) that is
A[2n+2], etc. So the diagonal entries in this format are A[0], A[n+1], A[2n+2], etc. What
happens to the lower part of the matrix? It is unreferenced, this means that it is not
referenced by the LAPACK code. Whatever has been there, remains there.

The packed storage format corresponds to LAPACK routines that begins with DPP (for
SPD) or DSP (for symmetric). (For example DPPTRF.) The packed storage is the
second one you have described. You only allocate the n(n+1)/2 entry to represent your
matrix A in memory. So in case of the upper case, A[0] is A(0,0), A[1] is A(0,1), A[2] is
A(1,1), etc. So the diagonal of A is in memory as A[0], A[2], A[5], A[9], A[14], etc.

The packed storage format offers an obvious advantage in term of memory requirement
as opposed to the full storage format. (You only need to allocate n^2/2 instead of n^2.)
However, for reasons that I'll skip, the packed storage format leads to extremely poor
performance (say 10 times slower or even more). So practically few people use
nowadays the packed storage.

There are some current research effort from Gustavson, Wasniewski et al. to introduce
new data storage in LAPACK that will lead to advantage of both (high performance +
optimal storage).

-julien

[susan_in_california]

Again, thank you! That was *very* clear, and it worked perfectly.

Best,
Susan