gonum.org/v1/gonum@v0.14.0/lapack/gonum/dpbtf2.go (about)

     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  //
    18  //	A = Uᵀ * U  if ul == blas.Upper
    19  //	A = L * Lᵀ  if ul == blas.Lower
    20  //
    21  // ul also specifies the storage of ab. If ul == blas.Upper, then
    22  // ab is stored as an upper-triangular banded matrix with kd super-diagonals,
    23  // and if ul == blas.Lower, ab is stored as a lower-triangular banded matrix
    24  // with kd sub-diagonals. On exit, the banded matrix U or L is stored in-place
    25  // into ab depending on the value of ul. Dpbtf2 returns whether the factorization
    26  // was successfully completed.
    27  //
    28  // The band storage scheme is illustrated below when n = 6, and kd = 2.
    29  // The resulting Cholesky decomposition is stored in the same elements as the
    30  // input band matrix (a11 becomes u11 or l11, etc.).
    31  //
    32  //	ul = blas.Upper
    33  //	a11 a12 a13
    34  //	a22 a23 a24
    35  //	a33 a34 a35
    36  //	a44 a45 a46
    37  //	a55 a56  *
    38  //	a66  *   *
    39  //
    40  //	ul = blas.Lower
    41  //	 *   *  a11
    42  //	 *  a21 a22
    43  //	a31 a32 a33
    44  //	a42 a43 a44
    45  //	a53 a54 a55
    46  //	a64 a65 a66
    47  //
    48  // Dpbtf2 is the unblocked version of the algorithm, see Dpbtrf for the blocked
    49  // version.
    50  //
    51  // Dpbtf2 is an internal routine, exported for testing purposes.
    52  func (Implementation) Dpbtf2(uplo blas.Uplo, n, kd int, ab []float64, ldab int) (ok bool) {
    53  	switch {
    54  	case uplo != blas.Upper && uplo != blas.Lower:
    55  		panic(badUplo)
    56  	case n < 0:
    57  		panic(nLT0)
    58  	case kd < 0:
    59  		panic(kdLT0)
    60  	case ldab < kd+1:
    61  		panic(badLdA)
    62  	}
    63  
    64  	// Quick return if possible.
    65  	if n == 0 {
    66  		return true
    67  	}
    68  
    69  	if len(ab) < (n-1)*ldab+kd+1 {
    70  		panic(shortAB)
    71  	}
    72  
    73  	bi := blas64.Implementation()
    74  
    75  	kld := max(1, ldab-1)
    76  	if uplo == blas.Upper {
    77  		// Compute the Cholesky factorization A = Uᵀ * U.
    78  		for j := 0; j < n; j++ {
    79  			// Compute U(j,j) and test for non-positive-definiteness.
    80  			ajj := ab[j*ldab]
    81  			if ajj <= 0 {
    82  				return false
    83  			}
    84  			ajj = math.Sqrt(ajj)
    85  			ab[j*ldab] = ajj
    86  			// Compute elements j+1:j+kn of row j and update the trailing submatrix
    87  			// within the band.
    88  			kn := min(kd, n-j-1)
    89  			if kn > 0 {
    90  				bi.Dscal(kn, 1/ajj, ab[j*ldab+1:], 1)
    91  				bi.Dsyr(blas.Upper, kn, -1, ab[j*ldab+1:], 1, ab[(j+1)*ldab:], kld)
    92  			}
    93  		}
    94  		return true
    95  	}
    96  	// Compute the Cholesky factorization A = L * Lᵀ.
    97  	for j := 0; j < n; j++ {
    98  		// Compute L(j,j) and test for non-positive-definiteness.
    99  		ajj := ab[j*ldab+kd]
   100  		if ajj <= 0 {
   101  			return false
   102  		}
   103  		ajj = math.Sqrt(ajj)
   104  		ab[j*ldab+kd] = ajj
   105  		// Compute elements j+1:j+kn of column j and update the trailing submatrix
   106  		// within the band.
   107  		kn := min(kd, n-j-1)
   108  		if kn > 0 {
   109  			bi.Dscal(kn, 1/ajj, ab[(j+1)*ldab+kd-1:], kld)
   110  			bi.Dsyr(blas.Lower, kn, -1, ab[(j+1)*ldab+kd-1:], kld, ab[(j+1)*ldab+kd:], kld)
   111  		}
   112  	}
   113  	return true
   114  }