## DSYTRF question

Open discussion regarding features, bugs, issues, vendors, etc.

### DSYTRF question

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

Posts: 39
Joined: Sun Jan 28, 2007 12:48 am

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.
kdallen

Posts: 39
Joined: Sun Jan 28, 2007 12:48 am

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.
sven

Posts: 146
Joined: Wed Dec 22, 2004 4:28 am

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.
kdallen

Posts: 39
Joined: Sun Jan 28, 2007 12:48 am

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.
sven

Posts: 146
Joined: Wed Dec 22, 2004 4:28 am

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?
kdallen

Posts: 39
Joined: Sun Jan 28, 2007 12:48 am

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:
Julien Langou

Posts: 835
Joined: Thu Dec 09, 2004 12:32 pm
Location: Denver, CO, USA

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
kdallen

Posts: 39
Joined: Sun Jan 28, 2007 12:48 am

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.
sven

Posts: 146
Joined: Wed Dec 22, 2004 4:28 am

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
kdallen

Posts: 39
Joined: Sun Jan 28, 2007 12:48 am

### Re:

kdallen wrote: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
holyjoly

Posts: 1
Joined: Mon Jul 16, 2012 5:24 pm