Implementing Low rank update of block-LDLᵀ decomposition

Post here if you have a question about LAPACK or ScaLAPACK algorithm or data format

Implementing Low rank update of block-LDLᵀ decomposition

Postby pushpendre » Thu Dec 29, 2016 3:30 pm


I was wondering if any stable algorithm for efficiently updating the block-LDLᵀ decomposition of symmetric indefinite matrices was already available in the LAPACK library? I know that the dsytf2* functions compute this decomposition but I was not able to find anything for updating the decomposition.

Professor Danny Sorensen's thesis "UPDATING THE SYMMETRIC INDEFINITE FACTORIZATION WITH APPLICATIONS IN A MODIFIED NEWTON'S METHOD"( contains one such algorithm including fortran source code on page 140, but the thesis is available as OCR'd pdf therefore I thought of asking on the LAPACK forum before spending time manually typing something that might already exist.

Posts: 3
Joined: Wed Dec 28, 2016 4:07 am

Re: Implementing Low rank update of block-LDLᵀ decomposition

Postby Julien Langou » Mon Jan 02, 2017 7:32 am

I understand you want to do: Given symmetric indefinite n-by-n A. Given n-by-k B with k << n. Given P, L and D, the LDLᵀ decomposition of A. Compute (in a fast way) the LDLᵀ decomposition of A + BBᵀ (by updating somehow P, L and D).

There is no such low rank update of LDLᵀ decomposition in LAPACK.

Note: You mention DSYTF2 to obtain an LDLᵀ decomposition from LAPACK. This is Level 2 BLAS using Bunch-Kaufman diagonal pivoting. Note that you can use DSYTRF ( Level 3 BLAS, Bunch-Kaufman diagonal pivoting ) and DSYTRF_RK (since v3.7, Level 3 BLAS, bounded Bunch-Kaufman (rook) diagonal pivoting). There is also DSYTRF_AA which performs LTLᵀ where T is symmetric tridiagonal (since v3.7, Level 3 BLAS, Aasen's factorization).

Note: DSYTRF_ROOK (since v3.5) should be deprecated any time soon, so please do not use.

If you want to contribute some codes, please contact us. ;)

Cheers, Julien.
Julien Langou
Posts: 810
Joined: Thu Dec 09, 2004 12:32 pm
Location: Denver, CO, USA

Return to Algorithm / Data

Who is online

Users browsing this forum: No registered users and 2 guests