github.com/gonum/lapack@v0.0.0-20181123203213-e4cdc5a0bff9/native/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 native
     6  
     7  import (
     8  	"math"
     9  
    10  	"github.com/gonum/blas"
    11  	"github.com/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^T U is stored in place into a. If ul == blas.Lower, then a = L L^T
    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  	if ul != blas.Upper && ul != blas.Lower {
    23  		panic(badUplo)
    24  	}
    25  	checkMatrix(n, n, a, lda)
    26  
    27  	if n == 0 {
    28  		return true
    29  	}
    30  
    31  	bi := blas64.Implementation()
    32  	if ul == blas.Upper {
    33  		for j := 0; j < n; j++ {
    34  			ajj := a[j*lda+j]
    35  			if j != 0 {
    36  				ajj -= bi.Ddot(j, a[j:], lda, a[j:], lda)
    37  			}
    38  			if ajj <= 0 || math.IsNaN(ajj) {
    39  				a[j*lda+j] = ajj
    40  				return false
    41  			}
    42  			ajj = math.Sqrt(ajj)
    43  			a[j*lda+j] = ajj
    44  			if j < n-1 {
    45  				bi.Dgemv(blas.Trans, j, n-j-1,
    46  					-1, a[j+1:], lda, a[j:], lda,
    47  					1, a[j*lda+j+1:], 1)
    48  				bi.Dscal(n-j-1, 1/ajj, a[j*lda+j+1:], 1)
    49  			}
    50  		}
    51  		return true
    52  	}
    53  	for j := 0; j < n; j++ {
    54  		ajj := a[j*lda+j]
    55  		if j != 0 {
    56  			ajj -= bi.Ddot(j, a[j*lda:], 1, a[j*lda:], 1)
    57  		}
    58  		if ajj <= 0 || math.IsNaN(ajj) {
    59  			a[j*lda+j] = ajj
    60  			return false
    61  		}
    62  		ajj = math.Sqrt(ajj)
    63  		a[j*lda+j] = ajj
    64  		if j < n-1 {
    65  			bi.Dgemv(blas.NoTrans, n-j-1, j,
    66  				-1, a[(j+1)*lda:], lda, a[j*lda:], 1,
    67  				1, a[(j+1)*lda+j:], lda)
    68  			bi.Dscal(n-j-1, 1/ajj, a[(j+1)*lda+j:], lda)
    69  		}
    70  	}
    71  	return true
    72  }