documentation bug: DSYTRF with UPLO='U'

[jaganmn]

The Bunch-Kaufman factorization accessing the upper triangular part of A yields A = U*D*U**T not A = U**T*D*U. There are three places in SRC/dsytrf.f where the factorization is stated incorrectly: lines 42, 145, and 262.

I'm not familiar with Aesen's algorithm. It might be worth checking whether SRC/dsytrf_aa.f is correct in stating the factorization as A = U**T*T*U. Unfortunately, LAWN 294 does not clarify ...

[langou]

Maybe this is a jargon thing but, in the upper case, LAPACK SYTRF returns a factorization of the type A = U**T*D*U. Maybe we should not call this a Bunch-Kaufman factorization, but in any case the comment that LAPACK SYTRF returns A = U**T*D*U is correct. If you have a reference that refers to a Bunch-Kaufman factorization being of the type A = U*D*U**T, let us know. This is certainly possible. We would "just" need to start at the lower right corner instead of starting at the upper left. LAPACK algorithms such as LU, Cholesky, LDLT always return lower times upper. (No matter the UPLO flag.) Returning upper times lower is certainly possible but we do not have algorithms who do this.

[ closing this issue, feel free to reopen if you have more comments ]

[jaganmn]

I've not seen a reference discussing the details of a Bunch-Kaufman factorization other than A = L*D*L**T, so I can't comment on the language choice. But the motivation of the report was my noticing that dsytri.f and dsytrs.f have:

https://github.com/Reference-LAPACK/lapack/blob/a83d8d28214f1a3760943d56e49c58a4e55ef13d/SRC/dsytri.f#L51-L52

https://github.com/Reference-LAPACK/lapack/blob/a83d8d28214f1a3760943d56e49c58a4e55ef13d/SRC/dsytrs.f#L51-L52

directly contradicting dsytrf.f. So I tested using R package Matrix, which now has methods (calling DSYTRF) for computing and expanding the Bunch-Kaufman factorization. Computing U*D*U**T reproduced A, computing U**T*D*U did not:

## zzz.R
## install.packages("Matrix", repos = "http://R-Forge.R-project.org")
library(Matrix)
(A <- matrix(c(1, 2, 2, 3), 2L, 2L))
(bk.A <- BunchKaufman(A, uplo = "U"))
(e.bk.A <- expand2(bk.A))
(U <- as(Reduce(`%*%`, e.bk.A[1:4]), "matrix"))
(D <- as(e.bk.A[[5L]], "matrix"))
norm(A - U %*% D %*% t(U), "2") / norm(A, "2")
norm(A - t(U) %*% D %*% U, "2") / norm(A, "2")

$ R -f zzz.R 
> ## zzz.R
> ## install.packages("Matrix", repos = "http://R-Forge.R-project.org")
> library(Matrix)
> (A <- matrix(c(1, 2, 2, 3), 2L, 2L))
     [,1] [,2]
[1,]    1    2
[2,]    2    3
> (bk.A <- BunchKaufman(A, uplo = "U"))
Bunch-Kaufman factorization of Formal class 'BunchKaufman' [package "Matrix"] with 5 slots
  ..@ uplo    : chr "U"
  ..@ x       : num [1:4] -0.333 0 0.667 3
  ..@ perm    : int [1:2] 1 2
  ..@ Dim     : int [1:2] 2 2
  ..@ Dimnames:List of 2
  .. ..$ : NULL
  .. ..$ : NULL
> (e.bk.A <- expand2(bk.A))
[[1]]
2 x 2 sparse Matrix of class "pMatrix"
        
[1,] | .
[2,] . |

[[2]]
2 x 2 sparse Matrix of class "dtCMatrix" (unitriangular)
                
[1,] I 0.6666667
[2,] . I        

[[3]]
2 x 2 sparse Matrix of class "pMatrix"
        
[1,] | .
[2,] . |

[[4]]
2 x 2 sparse Matrix of class "dtCMatrix" (unitriangular)
        
[1,] I .
[2,] . I

[[5]]
2 x 2 sparse Matrix of class "dsCMatrix"
                 
[1,] -0.3333333 .
[2,]  .         3

[[6]]
2 x 2 sparse Matrix of class "dtCMatrix" (unitriangular)
        
[1,] I .
[2,] . I

[[7]]
2 x 2 sparse Matrix of class "pMatrix"
        
[1,] | .
[2,] . |

[[8]]
2 x 2 sparse Matrix of class "dtCMatrix" (unitriangular)
                
[1,] I         .
[2,] 0.6666667 I

[[9]]
2 x 2 sparse Matrix of class "pMatrix"
        
[1,] | .
[2,] . |

> (U <- as(Reduce(`%*%`, e.bk.A[1:4]), "matrix"))
     [,1]      [,2]
[1,]    1 0.6666667
[2,]    0 1.0000000
> (D <- as(e.bk.A[[5L]], "matrix"))
           [,1] [,2]
[1,] -0.3333333    0
[2,]  0.0000000    3
> norm(A - U %*% D %*% t(U), "2") / norm(A, "2")
[1] 0
> norm(A - t(U) %*% D %*% U, "2") / norm(A, "2")
[1] 0.7177927
>

[jaganmn]

I should add that dsptrf.f has

https://github.com/Reference-LAPACK/lapack/blob/a83d8d28214f1a3760943d56e49c58a4e55ef13d/SRC/dsptrf.f#L117

and is consistent with dsptri.f and dsptrs.f. And my tests again indicate that U*D*U**T is correct there. So the language in dsytrf.f is really exceptional here ...