### DSYTRF question

Posted:

**Tue Apr 22, 2008 4:45 pm**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?

Page **1** of **1**

Posted: **Tue Apr 22, 2008 4:45 pm**

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?

Posted: **Wed Apr 23, 2008 2:02 am**

Golub has the following example:

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

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

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

Posted: **Wed Apr 23, 2008 9:47 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

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.

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.

Posted: **Wed Apr 23, 2008 5:11 pm**

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.

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.

Posted: **Thu Apr 24, 2008 12:15 pm**

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.

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.

Posted: **Fri Apr 25, 2008 4:13 pm**

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?

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?

Posted: **Fri Apr 25, 2008 4:28 pm**

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:

Posted: **Fri Apr 25, 2008 4:52 pm**

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

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

Posted: **Sun Apr 27, 2008 4:00 pm**

You can use the result that

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

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

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.

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

Posted: **Wed Apr 30, 2008 1:07 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

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

Posted: **Mon Jul 16, 2012 6:39 pm**

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