DSYTRF question

[kdallen]

Hi, I'm trying to get L and D from DSYTRF but can't decipher the explanation in the comments within the source code. Is there somewhere else that explains how L and D can be extracted?

[kdallen]

Golub has the following example:

Code: Select all

( 10 20  30 )   ( 1 0 0 ) (10 0 0 ) ( 1 2 3 )
( 20 45  80 ) = ( 2 1 0 ) ( 0 5 0 ) ( 0 1 4 )
( 30 80 171 )   ( 3 4 1 ) ( 0 0 1 ) ( 0 0 1 )


And DSYTRF yields the following factorization (0.25 and 0.062 are rounded for display):

Code: Select all

( 10    0     0 )
(  2   81     0 )
(  3 0.25 0.062 )

( 1 3 3 ) = IPIV

Can somone explain how to get the above 3 matrices from the DSYTRF factorization? Thanks.

[sven]

The short answer is that you cannot produce your three matrices from DSYTRF. This is because DSYTRF has performed a permutation at the second step, whereas your example is based upon no permutations.

In exact arithmetic DSYTRF would produce the factorization

P A P^T = L D L^T,

where

Code: Select all

   A = ( 10  20  30  )
       ( 20  45  80  )
       ( 30  80 171  )

   P = ( 1  0  0 ),  L = ( 1   0     0 ),  D = ( 10   0  0    ).
       ( 0  0  1 )       ( 3   1     0 )       (  0  81  0    )
       ( 0  1  0 )       ( 2  20/81  1 )       (  0   0  5/81 )


The vector IPIV = ( 1 3 3 ) means that no interchange (permutation) took place at steps 1 and 3, but at step 2, a permutation in the (2,3) plane took place.

Unless you wish to do something special, you generally should not need the explicit factors as routines such as DSYTRS take the output from DSYTRF.

You can see more information in Section 4.4 of the third edition of Golub and Van Loan. Hope this helps.

Sven Hammarling.

[kdallen]

Sven,

Thanks, I understand it now based on your explanation. Now the hard part. How to construct D when IPIV has negative values?

The reason I'm trying to construct L, D, (and P) is for handling constraints in orthogonal least squares problems. In ordinary LS problems w/ constraints A = J^T J, so things work out nicely. In orthogonal least squares, A is symmetric, sometimes positive definite, but not always, so LDL^T factorization comes in real handy.

[sven]

Two consecutive negative values in IPIV simply mean that D has a two by two block instead of a one by one block. The actual element in L that is used to store the off-diagonal element of D is zero.

If you have access to MATLAB, the a good way to see what is going on is to use a command such as

[L,D,P] = ldl(A)

which returns L, D and P such that

P^T A P = L D L^T.

You can compare the output from MATLAB with that from LAPACK.

Sven.

[kdallen]

Sven,

Thanks, almost there. No access to MATLAB, unfortunately, but an example in Golub w/ pivoting helped. I managed to factor a 12x12 matrix (w/ all 2x2 block diagonals for D), but 2 factorizations were required, one to determine P and the 2nd to factor the permutation.

Is it possible to determine P w/o factoring the original matrix?

[Julien Langou]

Is it possible to determine P w/o factoring the original matrix?

Nope, the easiest (known) way to find P is to factor the original matrix.
If you find a less expensive way to get P, let us know !!!!! :wink:

[kdallen]

Julien,

Thanks. I was afraid of that. I did manage to find a link that explains the gory details:

http://www-curri.u-strasbg.fr/documenta ... #m86Plchri

Example 3-2. Symmetric indefinite matrix factorization

But, the method looks REAL painful. I guess the choices are:

1. Perform two factorizations. Straightforward, but 2x computations needed for the factorization.
2. Use the L = P1*L1*P2*L2, etc. method, though L is really PL. Looks very ugly
3. Forget about using LDL^T and just calculate A^-1

[sven]

You can use the result that

Code: Select all

L(r) P(p,q) = P(p,q) M(r),  r < p,q


where M(r) is the unit lower triangular matrix L(r) with elements l(p,r) and l(q,r) interchanged. Moving the permutations matrices in this way is what allows us to express the factorization as

Code: Select all

P^T A P = L D L^T


The same idea is used to express Gaussian elimination with partial pivoting as

Code: Select all

P A = L U


Of course here we have the minor complication of two by two blocks, but I hope this helps you to see how to readily extract the factorization from the information returned by LAPACK.

Good luck,

Sven.

[kdallen]

Sven,

Thanks for all your help. I was able to come up with a concise algorithm for forming P, L, D based on element swapping of DSYTRF output. :D

[holyjoly]

I was able to come up with a concise algorithm for forming P, L, D based on element swapping of DSYTRF output.

Hi,

I realise that it's been a while since you did this, but would you mind outlining the algorithm you came up with?

Many thanks,
Joe