Confused by dsytrf documentation

Jason_Jones · ‎07-22-2014

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:

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?

mecej4 · ‎07-23-2014

Most users probably follow up a call to ?SYTRF with a call to ?SYTRS in order to solve a symmetric indefinite system. It is when we wish to intercept and interpret the intermediate results (the factors and the permutation vector) that the complexities confront us. Do you really need to look at the intermediate results? If so, please read on ...

The example that you ran gives a purely diagonal factor for D -- there are no blocks larger than 1 X 1 involved, and all the elements of the permutation vector are positive. You may find it helpful to compare the MKL results with the results from Matlab. If you do

[L,D,p]=ldl(M,'vector')

after entering your matrix as M in Matlab, you will get

L =

    1.0000         0         0         0
   -1.2548    1.0000         0         0
   -0.4004    0.5735    1.0000         0
   -0.1930    0.5097   -0.2054    1.0000


D =

    1.4527         0         0         0
         0   12.4156         0         0
         0         0    0.3765         0
         0         0         0    0.0004


p =

     1     3     2     4

Note that the encodings of the permutations are slightly different in Matlab and MKL. Matlab generates the permuted factorization A(p,p) = L D L^T. Compare these results with the MKL results to discover how the results are packed into the input/output 2-D array.

Jason_Jones · ‎07-23-2014

Thanks for that answer! I have dug in a little more and worked through some examples with 2x2 blocks and got that to work. As you suggested, I did the calculation in Matlab. When I call:

[ L, D, P ] = ldl(A);
P*L*D*P' - A;

I get an answer of zero (close to machine precision). The P and D appear to be the same from dsytrf and from Matlab, however, I am seeing a strange difference in the ordering of the elements of L. In Matlab, I get:

    1.0000         0         0         0
   -1.2548**  1.0000         0         0
   -0.4004**  0.5735    1.0000         0
   -0.1930    0.5097   -0.2054    1.0000

However, from 'dsytrf', I get L to be:

    1.0000         0         0         0 
   -0.4004**  1.0000         0         0 
   -1.2548**  0.5735    1.0000         0 
   -0.1930    0.5097   -0.2054    1.0000

If I multiply Matlab's P*L, I get:

    1.0000         0         0         0
   -0.4004**  0.5735    1.0000         0
   -1.2548**  1.0000         0         0
   -0.1930    0.5097   -0.2054    1.0000

Whereas if I multiply 'dsytrf's P*L, I get:

   1.0000        0        0        0 
  -1.2548** 0.5735   1.0000        0 
  -0.4004** 1.0000        0        0 
  -0.1930   0.5097  -0.2054   1.0000

It looks to me like 'dsytrf' is giving me the wrong answer, but I assume it is either a bug in my part of the code (though I have looked and looked) or I am not understanding something. Is it possible that P * L * D * L' * P' does not equal A, but I have the correct answer?

For the record, 'ipiv' is returning:

which I am assuming indicates a permutation matrix that looks like this:

1   0   0   0
0   0   1   0
0   1   0   0
0   0   0   1

When I calculate P * L * D * L' * P' - A, from 'dsytrf', I get:

       0  -1.2412   1.2412        0 
 -1.2412   2.0545        0   0.2395 
  1.2412        0  -2.0545  -0.2395 
       0   0.2395  -0.2395        0

mecej4 · ‎07-23-2014

You will have to read carefully as to the details of how the permutations are encoded in the MKL manual. To undo a sequence of row/column interchanges, the interchanges need to be applied in reverse order, since it is possible that a given row/column has been exchanged more than once. For your matrix, the permutation array IPIV is [1,3,3,4], which means that only one exchange was needed. Again, please read the MKL manual to ascertain what entities were interchanged.