vendor/gonum.org/v1/gonum/lapack/gonum/dpbtf2.go

   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.
   4
   5 package gonum
   6
   7 import (
   8         "math"
   9
  10         "gonum.org/v1/gonum/blas"
  11         "gonum.org/v1/gonum/blas/blas64"
  12 )
  13
  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.
  25 //
  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.).
  29 //
  30 //  ul = blas.Upper
  31 //  a11 a12 a13
  32 //  a22 a23 a24
  33 //  a33 a34 a35
  34 //  a44 a45 a46
  35 //  a55 a56  *
  36 //  a66  *   *
  37 //
  38 //  ul = blas.Lower
  39 //   *   *  a11
  40 //   *  a21 a22
  41 //  a31 a32 a33
  42 //  a42 a43 a44
  43 //  a53 a54 a55
  44 //  a64 a65 a66
  45 //
  46 // Dpbtf2 is the unblocked version of the algorithm, see Dpbtrf for the blocked
  47 // version.
  48 //
  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 {
  52                 panic(badUplo)
  53         }
  54         checkSymBanded(ab, n, kd, ldab)
  55         if n == 0 {
  56                 return
  57         }
  58         bi := blas64.Implementation()
  59         kld := max(1, ldab-1)
  60         if ul == blas.Upper {
  61                 for j := 0; j < n; j++ {
  62                         // Compute U(J,J) and test for non positive-definiteness.
  63                         ajj := ab[j*ldab]
  64                         if ajj <= 0 {
  65                                 return false
  66                         }
  67                         ajj = math.Sqrt(ajj)
  68                         ab[j*ldab] = ajj
  69                         // Compute elements j+1:j+kn of row J and update the trailing submatrix
  70                         // within the band.
  71                         kn := min(kd, n-j-1)
  72                         if kn > 0 {
  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)
  75                         }
  76                 }
  77                 return true
  78         }
  79         for j := 0; j < n; j++ {
  80                 // Compute L(J,J) and test for non positive-definiteness.
  81                 ajj := ab[j*ldab+kd]
  82                 if ajj <= 0 {
  83                         return false
  84                 }
  85                 ajj = math.Sqrt(ajj)
  86                 ab[j*ldab+kd] = ajj
  87
  88                 // Compute elements J+1:J+KN of column J and update the trailing submatrix
  89                 // within the band.
  90                 kn := min(kd, n-j-1)
  91                 if kn > 0 {
  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)
  94                 }
  95         }
  96         return true
  97 }