I have created a wrapper/helper function to call 'dsytrf'. It works (i.e l*d*l' == A) in all cases in which 'ipiv' is simply sequential and all the "blocks" of the "block diagonal" 'D' are 1x1. However, the other case I cannot figure out. I am using 'dystrf' to factorize this matrix, 'A':

1.4527167 -0.58173740 -1.8229176 -0.28031340 -0.58173740 4.6931840 7.8504984 3.6645482 -1.8229176 7.8504984 14.703079 6.6805245 -0.28031340 3.6645482 6.6805245 3.2963768

The return that I get from 'dsytrf' is:

1.452717 -0.400448 -1.254834 -0.192958 -0.581737 12.415621 0.573513 0.509743 -1.822918 7.850498 0.376524 -0.205396 -0.280313 3.664548 6.680525 0.000352

and the 'ipiv' vector contains:

1 3 3 4

The remarks for the 'ipiv' output say that:

If ipiv(i) = k >0, then d

_{ii}is a 1-by-1 block, and the i-th row and column of A was interchanged with the k-th row and column.

I have taken that to mean that if ipiv(1) = 1, ipiv(2) = 2, ipiv(3) = 3, etc for all ipiv, then all of D's diagonal blocks are 1x1 and the permutation matrix is simply the identity matrix so that L*D*L' should equal A. However, as I read further on, I see:

If uplo = 'L' and ipiv(i) =ipiv(i+1) = -m < 0, then D has a 2-by-2 block in rows/columns i and i+1, and (i+1)-th row and column of A was interchanged with the m-th row and column.

I have taken this to mean that since the second and third elements of 'ipiv' both equal 3, then D has a 2x2 block at D(2:3, 2:3) and that I need to use a permutation matrix. The instructions for filling in L and D are further down in the document:

If s = 1, D(k) overwrites A(k,k), and v overwrites A(k+1:n,k).

If s = 2, the lower triangle of D(k) overwrites A(k,k), A(k+1,k), and A(k+1,k+1), and v overwrites A(k+2:n,k:k+1).

This sounds to me as if in the s=1 case, D is the diagonal of the output A from dsytrf and L is the lower triangular portion. This has been working for me. I assume that the s=2 case is the case where D has 2x2 blocks. And specifically in this case, from the dystrf output A: D(2,2) = A(2,2), D(3,2) = A(3,2), L(4, 3) = A(4,3). It does not mention what the upper triangular portion of the 2x2 block in D is. Is it supposed to be the same as the lower triangular portion i.e. is D symmetric, or is it supposed to be 0, i.e. is D a lower-triangular-block-matrix? Here are my best guess for L and D with A for reference:

returned 'A' = 1.452717 -0.400448 -1.254834 -0.192958 -0.581737 12.415621 0.573513 0.509743 -1.822918 7.850498 0.376524 -0.205396 -0.280313 3.664548 6.680525 0.000352

best guess D = 1.4527167 0.0 0.0 0.0 0.0 12.415621 0.0 0.0 0.0 0.57351258 0.37652362 0.0 0.0 0.0 0.0 0.00035225881

best guess L = 1.0 0.0 0.0 0.0 -0.40044793 1.0 0.0 0.0 -1.2548335 0.0 1.0 0.0 -0.19295806 0.50974316 -0.20539592 1.0

What am I doing wrong?