1 // Copyright ©2017 The Gonum Authors. All rights reserved.
2 // Use of this source code is governed by a BSD-style
3 // license that can be found in the LICENSE file.
10 "gonum.org/v1/gonum/blas"
11 "gonum.org/v1/gonum/blas/blas64"
14 // Dpbtf2 computes the Cholesky factorization of a symmetric positive banded
15 // matrix ab. The matrix ab is n×n with kd diagonal bands. The Cholesky
16 // factorization computed is
17 // A = U^T * U if ul == blas.Upper
18 // A = L * L^T if ul == blas.Lower
19 // ul also specifies the storage of ab. If ul == blas.Upper, then
20 // ab is stored as an upper-triangular banded matrix with kd super-diagonals,
21 // and if ul == blas.Lower, ab is stored as a lower-triangular banded matrix
22 // with kd sub-diagonals. On exit, the banded matrix U or L is stored in-place
23 // into ab depending on the value of ul. Dpbtf2 returns whether the factorization
24 // was successfully completed.
26 // The band storage scheme is illustrated below when n = 6, and kd = 2.
27 // The resulting Cholesky decomposition is stored in the same elements as the
28 // input band matrix (a11 becomes u11 or l11, etc.).
46 // Dpbtf2 is the unblocked version of the algorithm, see Dpbtrf for the blocked
49 // Dpbtf2 is an internal routine, exported for testing purposes.
50 func (Implementation) Dpbtf2(ul blas.Uplo, n, kd int, ab []float64, ldab int) (ok bool) {
51 if ul != blas.Upper && ul != blas.Lower {
54 checkSymBanded(ab, n, kd, ldab)
58 bi := blas64.Implementation()
61 for j := 0; j < n; j++ {
62 // Compute U(J,J) and test for non positive-definiteness.
69 // Compute elements j+1:j+kn of row J and update the trailing submatrix
73 bi.Dscal(kn, 1/ajj, ab[j*ldab+1:], 1)
74 bi.Dsyr(blas.Upper, kn, -1, ab[j*ldab+1:], 1, ab[(j+1)*ldab:], kld)
79 for j := 0; j < n; j++ {
80 // Compute L(J,J) and test for non positive-definiteness.
88 // Compute elements J+1:J+KN of column J and update the trailing submatrix
92 bi.Dscal(kn, 1/ajj, ab[(j+1)*ldab+kd-1:], kld)
93 bi.Dsyr(blas.Lower, kn, -1, ab[(j+1)*ldab+kd-1:], kld, ab[(j+1)*ldab+kd:], kld)