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

     1  // Copyright ©2019 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  	"gonum.org/v1/gonum/blas"
     9  	"gonum.org/v1/gonum/blas/blas64"
    10  )
    11  
    12  // Dpbtrs solves a system of linear equations A*X = B with an n×n symmetric
    13  // positive definite band matrix A using the Cholesky factorization
    14  //
    15  //	A = Uᵀ * U  if uplo == blas.Upper
    16  //	A = L * Lᵀ  if uplo == blas.Lower
    17  //
    18  // computed by Dpbtrf. kd is the number of super- or sub-diagonals of A. See the
    19  // documentation for Dpbtrf for a description of the band storage format of A.
    20  //
    21  // On entry, b contains the n×nrhs right hand side matrix B. On return, it is
    22  // overwritten with the solution matrix X.
    23  func (Implementation) Dpbtrs(uplo blas.Uplo, n, kd, nrhs int, ab []float64, ldab int, b []float64, ldb int) {
    24  	switch {
    25  	case uplo != blas.Upper && uplo != blas.Lower:
    26  		panic(badUplo)
    27  	case n < 0:
    28  		panic(nLT0)
    29  	case kd < 0:
    30  		panic(kdLT0)
    31  	case nrhs < 0:
    32  		panic(nrhsLT0)
    33  	case ldab < kd+1:
    34  		panic(badLdA)
    35  	case ldb < max(1, nrhs):
    36  		panic(badLdB)
    37  	}
    38  
    39  	// Quick return if possible.
    40  	if n == 0 || nrhs == 0 {
    41  		return
    42  	}
    43  
    44  	if len(ab) < (n-1)*ldab+kd+1 {
    45  		panic(shortAB)
    46  	}
    47  	if len(b) < (n-1)*ldb+nrhs {
    48  		panic(shortB)
    49  	}
    50  
    51  	bi := blas64.Implementation()
    52  	if uplo == blas.Upper {
    53  		// Solve A*X = B where A = Uᵀ*U.
    54  		for j := 0; j < nrhs; j++ {
    55  			// Solve Uᵀ*Y = B, overwriting B with Y.
    56  			bi.Dtbsv(blas.Upper, blas.Trans, blas.NonUnit, n, kd, ab, ldab, b[j:], ldb)
    57  			// Solve U*X = Y, overwriting Y with X.
    58  			bi.Dtbsv(blas.Upper, blas.NoTrans, blas.NonUnit, n, kd, ab, ldab, b[j:], ldb)
    59  		}
    60  	} else {
    61  		// Solve A*X = B where A = L*Lᵀ.
    62  		for j := 0; j < nrhs; j++ {
    63  			// Solve L*Y = B, overwriting B with Y.
    64  			bi.Dtbsv(blas.Lower, blas.NoTrans, blas.NonUnit, n, kd, ab, ldab, b[j:], ldb)
    65  			// Solve Lᵀ*X = Y, overwriting Y with X.
    66  			bi.Dtbsv(blas.Lower, blas.Trans, blas.NonUnit, n, kd, ab, ldab, b[j:], ldb)
    67  		}
    68  	}
    69  }