gonum.org/v1/gonum@v0.15.1-0.20240517103525-f853624cb1bb/lapack/gonum/dpttrf.go (about)

     1  // Copyright ©2023 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  // Dpttrf computes the L*D*Lᵀ factorization of an n×n symmetric positive
     8  // definite tridiagonal matrix A and returns whether the factorization was
     9  // successful.
    10  //
    11  // On entry, d and e contain the n diagonal and (n-1) subdiagonal elements,
    12  // respectively, of A.
    13  //
    14  // On return, d contains the n diagonal elements of the diagonal matrix D and e
    15  // contains the (n-1) subdiagonal elements of the unit bidiagonal matrix L.
    16  func (impl Implementation) Dpttrf(n int, d, e []float64) (ok bool) {
    17  	if n < 0 {
    18  		panic(nLT0)
    19  	}
    20  
    21  	if n == 0 {
    22  		return true
    23  	}
    24  
    25  	switch {
    26  	case len(d) < n:
    27  		panic(shortD)
    28  	case len(e) < n-1:
    29  		panic(shortE)
    30  	}
    31  
    32  	// Compute the L*D*Lᵀ (or Uᵀ*D*U) factorization of A.
    33  	i4 := (n - 1) % 4
    34  	for i := 0; i < i4; i++ {
    35  		if d[i] <= 0 {
    36  			return false
    37  		}
    38  		ei := e[i]
    39  		e[i] /= d[i]
    40  		d[i+1] -= e[i] * ei
    41  	}
    42  	for i := i4; i < n-4; i += 4 {
    43  		// Drop out of the loop if d[i] <= 0: the matrix is not positive
    44  		// definite.
    45  		if d[i] <= 0 {
    46  			return false
    47  		}
    48  
    49  		// Solve for e[i] and d[i+1].
    50  		ei := e[i]
    51  		e[i] /= d[i]
    52  		d[i+1] -= e[i] * ei
    53  		if d[i+1] <= 0 {
    54  			return false
    55  		}
    56  
    57  		// Solve for e[i+1] and d[i+2].
    58  		ei = e[i+1]
    59  		e[i+1] /= d[i+1]
    60  		d[i+2] -= e[i+1] * ei
    61  		if d[i+2] <= 0 {
    62  			return false
    63  		}
    64  
    65  		// Solve for e[i+2] and d[i+3].
    66  		ei = e[i+2]
    67  		e[i+2] /= d[i+2]
    68  		d[i+3] -= e[i+2] * ei
    69  		if d[i+3] <= 0 {
    70  			return false
    71  		}
    72  
    73  		// Solve for e[i+3] and d[i+4].
    74  		ei = e[i+3]
    75  		e[i+3] /= d[i+3]
    76  		d[i+4] -= e[i+3] * ei
    77  	}
    78  	// Check d[n-1] for positive definiteness.
    79  	return d[n-1] > 0
    80  }