- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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 dii 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?
Link Copied
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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 LT. Compare these results with the MKL results to discover how the results are packed into the input/output 2-D array.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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:
1 3 3 4
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
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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.

- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page