github.com/gopherd/gonum@v0.0.4/lapack/gonum/dpotf2.go (about)

     1  // Copyright ©2015 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  	"github.com/gopherd/gonum/blas"
    11  	"github.com/gopherd/gonum/blas/blas64"
    12  )
    13  
    14  // Dpotf2 computes the Cholesky decomposition of the symmetric positive definite
    15  // matrix a. If ul == blas.Upper, then a is stored as an upper-triangular matrix,
    16  // and a = Uᵀ U is stored in place into a. If ul == blas.Lower, then a = L Lᵀ
    17  // is computed and stored in-place into a. If a is not positive definite, false
    18  // is returned. This is the unblocked version of the algorithm.
    19  //
    20  // Dpotf2 is an internal routine. It is exported for testing purposes.
    21  func (Implementation) Dpotf2(ul blas.Uplo, n int, a []float64, lda int) (ok bool) {
    22  	switch {
    23  	case ul != blas.Upper && ul != blas.Lower:
    24  		panic(badUplo)
    25  	case n < 0:
    26  		panic(nLT0)
    27  	case lda < max(1, n):
    28  		panic(badLdA)
    29  	}
    30  
    31  	// Quick return if possible.
    32  	if n == 0 {
    33  		return true
    34  	}
    35  
    36  	if len(a) < (n-1)*lda+n {
    37  		panic(shortA)
    38  	}
    39  
    40  	bi := blas64.Implementation()
    41  
    42  	if ul == blas.Upper {
    43  		for j := 0; j < n; j++ {
    44  			ajj := a[j*lda+j]
    45  			if j != 0 {
    46  				ajj -= bi.Ddot(j, a[j:], lda, a[j:], lda)
    47  			}
    48  			if ajj <= 0 || math.IsNaN(ajj) {
    49  				a[j*lda+j] = ajj
    50  				return false
    51  			}
    52  			ajj = math.Sqrt(ajj)
    53  			a[j*lda+j] = ajj
    54  			if j < n-1 {
    55  				bi.Dgemv(blas.Trans, j, n-j-1,
    56  					-1, a[j+1:], lda, a[j:], lda,
    57  					1, a[j*lda+j+1:], 1)
    58  				bi.Dscal(n-j-1, 1/ajj, a[j*lda+j+1:], 1)
    59  			}
    60  		}
    61  		return true
    62  	}
    63  	for j := 0; j < n; j++ {
    64  		ajj := a[j*lda+j]
    65  		if j != 0 {
    66  			ajj -= bi.Ddot(j, a[j*lda:], 1, a[j*lda:], 1)
    67  		}
    68  		if ajj <= 0 || math.IsNaN(ajj) {
    69  			a[j*lda+j] = ajj
    70  			return false
    71  		}
    72  		ajj = math.Sqrt(ajj)
    73  		a[j*lda+j] = ajj
    74  		if j < n-1 {
    75  			bi.Dgemv(blas.NoTrans, n-j-1, j,
    76  				-1, a[(j+1)*lda:], lda, a[j*lda:], 1,
    77  				1, a[(j+1)*lda+j:], lda)
    78  			bi.Dscal(n-j-1, 1/ajj, a[(j+1)*lda+j:], lda)
    79  		}
    80  	}
    81  	return true
    82  }