github.com/jingcheng-WU/gonum@v0.9.1-0.20210323123734-f1a2a11a8f7b/lapack/gonum/dpbcon.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  	"math"
     9  
    10  	"github.com/jingcheng-WU/gonum/blas"
    11  	"github.com/jingcheng-WU/gonum/blas/blas64"
    12  )
    13  
    14  // Dpbcon returns an estimate of the reciprocal of the condition number (in the
    15  // 1-norm) of an n×n symmetric positive definite band matrix using the Cholesky
    16  // factorization
    17  //  A = Uᵀ*U  if uplo == blas.Upper
    18  //  A = L*Lᵀ  if uplo == blas.Lower
    19  // computed by Dpbtrf. The estimate is obtained for norm(inv(A)), and the
    20  // reciprocal of the condition number is computed as
    21  //  rcond = 1 / (anorm * norm(inv(A))).
    22  //
    23  // The length of work must be at least 3*n and the length of iwork must be at
    24  // least n.
    25  func (impl Implementation) Dpbcon(uplo blas.Uplo, n, kd int, ab []float64, ldab int, anorm float64, work []float64, iwork []int) (rcond float64) {
    26  	switch {
    27  	case uplo != blas.Upper && uplo != blas.Lower:
    28  		panic(badUplo)
    29  	case n < 0:
    30  		panic(nLT0)
    31  	case kd < 0:
    32  		panic(kdLT0)
    33  	case ldab < kd+1:
    34  		panic(badLdA)
    35  	case anorm < 0:
    36  		panic(badNorm)
    37  	}
    38  
    39  	// Quick return if possible.
    40  	if n == 0 {
    41  		return 1
    42  	}
    43  
    44  	switch {
    45  	case len(ab) < (n-1)*ldab+kd+1:
    46  		panic(shortAB)
    47  	case len(work) < 3*n:
    48  		panic(shortWork)
    49  	case len(iwork) < n:
    50  		panic(shortIWork)
    51  	}
    52  
    53  	// Quick return if possible.
    54  	if anorm == 0 {
    55  		return 0
    56  	}
    57  
    58  	const smlnum = dlamchS
    59  
    60  	var (
    61  		ainvnm float64
    62  		kase   int
    63  		isave  [3]int
    64  		normin bool
    65  
    66  		// Denote work slices.
    67  		x     = work[:n]
    68  		v     = work[n : 2*n]
    69  		cnorm = work[2*n : 3*n]
    70  	)
    71  	// Estimate the 1-norm of the inverse.
    72  	bi := blas64.Implementation()
    73  	for {
    74  		ainvnm, kase = impl.Dlacn2(n, v, x, iwork, ainvnm, kase, &isave)
    75  		if kase == 0 {
    76  			break
    77  		}
    78  		var op1, op2 blas.Transpose
    79  		if uplo == blas.Upper {
    80  			// Multiply x by inv(Uᵀ),
    81  			op1 = blas.Trans
    82  			// then by inv(Uᵀ).
    83  			op2 = blas.NoTrans
    84  		} else {
    85  			// Multiply x by inv(L),
    86  			op1 = blas.NoTrans
    87  			// then by inv(Lᵀ).
    88  			op2 = blas.Trans
    89  		}
    90  		scaleL := impl.Dlatbs(uplo, op1, blas.NonUnit, normin, n, kd, ab, ldab, x, cnorm)
    91  		normin = true
    92  		scaleU := impl.Dlatbs(uplo, op2, blas.NonUnit, normin, n, kd, ab, ldab, x, cnorm)
    93  		// Multiply x by 1/scale if doing so will not cause overflow.
    94  		scale := scaleL * scaleU
    95  		if scale != 1 {
    96  			ix := bi.Idamax(n, x, 1)
    97  			if scale < math.Abs(x[ix])*smlnum || scale == 0 {
    98  				return 0
    99  			}
   100  			impl.Drscl(n, scale, x, 1)
   101  		}
   102  	}
   103  	if ainvnm == 0 {
   104  		return 0
   105  	}
   106  	// Return the estimate of the reciprocal condition number.
   107  	return (1 / ainvnm) / anorm
   108  }