github.com/gonum/lapack@v0.0.0-20181123203213-e4cdc5a0bff9/native/dpotrf.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  	"github.com/gonum/blas"
     9  	"github.com/gonum/blas/blas64"
    10  )
    11  
    12  // Dpotrf computes the Cholesky decomposition of the symmetric positive definite
    13  // matrix a. If ul == blas.Upper, then a is stored as an upper-triangular matrix,
    14  // and a = U^T U is stored in place into a. If ul == blas.Lower, then a = L L^T
    15  // is computed and stored in-place into a. If a is not positive definite, false
    16  // is returned. This is the blocked version of the algorithm.
    17  func (impl Implementation) Dpotrf(ul blas.Uplo, n int, a []float64, lda int) (ok bool) {
    18  	if ul != blas.Upper && ul != blas.Lower {
    19  		panic(badUplo)
    20  	}
    21  	checkMatrix(n, n, a, lda)
    22  
    23  	if n == 0 {
    24  		return true
    25  	}
    26  
    27  	nb := impl.Ilaenv(1, "DPOTRF", string(ul), n, -1, -1, -1)
    28  	if nb <= 1 || n <= nb {
    29  		return impl.Dpotf2(ul, n, a, lda)
    30  	}
    31  	bi := blas64.Implementation()
    32  	if ul == blas.Upper {
    33  		for j := 0; j < n; j += nb {
    34  			jb := min(nb, n-j)
    35  			bi.Dsyrk(blas.Upper, blas.Trans, jb, j,
    36  				-1, a[j:], lda,
    37  				1, a[j*lda+j:], lda)
    38  			ok = impl.Dpotf2(blas.Upper, jb, a[j*lda+j:], lda)
    39  			if !ok {
    40  				return ok
    41  			}
    42  			if j+jb < n {
    43  				bi.Dgemm(blas.Trans, blas.NoTrans, jb, n-j-jb, j,
    44  					-1, a[j:], lda, a[j+jb:], lda,
    45  					1, a[j*lda+j+jb:], lda)
    46  				bi.Dtrsm(blas.Left, blas.Upper, blas.Trans, blas.NonUnit, jb, n-j-jb,
    47  					1, a[j*lda+j:], lda,
    48  					a[j*lda+j+jb:], lda)
    49  			}
    50  		}
    51  		return true
    52  	}
    53  	for j := 0; j < n; j += nb {
    54  		jb := min(nb, n-j)
    55  		bi.Dsyrk(blas.Lower, blas.NoTrans, jb, j,
    56  			-1, a[j*lda:], lda,
    57  			1, a[j*lda+j:], lda)
    58  		ok := impl.Dpotf2(blas.Lower, jb, a[j*lda+j:], lda)
    59  		if !ok {
    60  			return ok
    61  		}
    62  		if j+jb < n {
    63  			bi.Dgemm(blas.NoTrans, blas.Trans, n-j-jb, jb, j,
    64  				-1, a[(j+jb)*lda:], lda, a[j*lda:], lda,
    65  				1, a[(j+jb)*lda+j:], lda)
    66  			bi.Dtrsm(blas.Right, blas.Lower, blas.Trans, blas.NonUnit, n-j-jb, jb,
    67  				1, a[j*lda+j:], lda,
    68  				a[(j+jb)*lda+j:], lda)
    69  		}
    70  	}
    71  	return true
    72  }