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