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