github.com/jingcheng-WU/gonum@v0.9.1-0.20210323123734-f1a2a11a8f7b/lapack/gonum/dtrcon.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 gonum
     6  
     7  import (
     8  	"math"
     9  
    10  	"github.com/jingcheng-WU/gonum/blas"
    11  	"github.com/jingcheng-WU/gonum/blas/blas64"
    12  	"github.com/jingcheng-WU/gonum/lapack"
    13  )
    14  
    15  // Dtrcon estimates the reciprocal of the condition number of a triangular matrix A.
    16  // The condition number computed may be based on the 1-norm or the ∞-norm.
    17  //
    18  // work is a temporary data slice of length at least 3*n and Dtrcon will panic otherwise.
    19  //
    20  // iwork is a temporary data slice of length at least n and Dtrcon will panic otherwise.
    21  func (impl Implementation) Dtrcon(norm lapack.MatrixNorm, uplo blas.Uplo, diag blas.Diag, n int, a []float64, lda int, work []float64, iwork []int) float64 {
    22  	switch {
    23  	case norm != lapack.MaxColumnSum && norm != lapack.MaxRowSum:
    24  		panic(badNorm)
    25  	case uplo != blas.Upper && uplo != blas.Lower:
    26  		panic(badUplo)
    27  	case diag != blas.NonUnit && diag != blas.Unit:
    28  		panic(badDiag)
    29  	case n < 0:
    30  		panic(nLT0)
    31  	case lda < max(1, n):
    32  		panic(badLdA)
    33  	}
    34  
    35  	if n == 0 {
    36  		return 1
    37  	}
    38  
    39  	switch {
    40  	case len(a) < (n-1)*lda+n:
    41  		panic(shortA)
    42  	case len(work) < 3*n:
    43  		panic(shortWork)
    44  	case len(iwork) < n:
    45  		panic(shortIWork)
    46  	}
    47  
    48  	bi := blas64.Implementation()
    49  
    50  	var rcond float64
    51  	smlnum := dlamchS * float64(n)
    52  
    53  	anorm := impl.Dlantr(norm, uplo, diag, n, n, a, lda, work)
    54  
    55  	if anorm <= 0 {
    56  		return rcond
    57  	}
    58  	var ainvnm float64
    59  	var normin bool
    60  	kase1 := 2
    61  	if norm == lapack.MaxColumnSum {
    62  		kase1 = 1
    63  	}
    64  	var kase int
    65  	isave := new([3]int)
    66  	var scale float64
    67  	for {
    68  		ainvnm, kase = impl.Dlacn2(n, work[n:], work, iwork, ainvnm, kase, isave)
    69  		if kase == 0 {
    70  			if ainvnm != 0 {
    71  				rcond = (1 / anorm) / ainvnm
    72  			}
    73  			return rcond
    74  		}
    75  		if kase == kase1 {
    76  			scale = impl.Dlatrs(uplo, blas.NoTrans, diag, normin, n, a, lda, work, work[2*n:])
    77  		} else {
    78  			scale = impl.Dlatrs(uplo, blas.Trans, diag, normin, n, a, lda, work, work[2*n:])
    79  		}
    80  		normin = true
    81  		if scale != 1 {
    82  			ix := bi.Idamax(n, work, 1)
    83  			xnorm := math.Abs(work[ix])
    84  			if scale == 0 || scale < xnorm*smlnum {
    85  				return rcond
    86  			}
    87  			impl.Drscl(n, scale, work, 1)
    88  		}
    89  	}
    90  }