github.com/gonum/lapack@v0.0.0-20181123203213-e4cdc5a0bff9/native/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 native
     6  
     7  import (
     8  	"math"
     9  
    10  	"github.com/gonum/blas"
    11  	"github.com/gonum/blas/blas64"
    12  	"github.com/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  	if norm != lapack.MaxColumnSum && norm != lapack.MaxRowSum {
    23  		panic(badNorm)
    24  	}
    25  	if uplo != blas.Upper && uplo != blas.Lower {
    26  		panic(badUplo)
    27  	}
    28  	if diag != blas.NonUnit && diag != blas.Unit {
    29  		panic(badDiag)
    30  	}
    31  	if len(work) < 3*n {
    32  		panic(badWork)
    33  	}
    34  	if len(iwork) < n {
    35  		panic(badWork)
    36  	}
    37  	if n == 0 {
    38  		return 1
    39  	}
    40  	bi := blas64.Implementation()
    41  
    42  	var rcond float64
    43  	smlnum := dlamchS * float64(n)
    44  
    45  	anorm := impl.Dlantr(norm, uplo, diag, n, n, a, lda, work)
    46  
    47  	if anorm <= 0 {
    48  		return rcond
    49  	}
    50  	var ainvnm float64
    51  	var normin bool
    52  	kase1 := 2
    53  	if norm == lapack.MaxColumnSum {
    54  		kase1 = 1
    55  	}
    56  	var kase int
    57  	isave := new([3]int)
    58  	var scale float64
    59  	for {
    60  		ainvnm, kase = impl.Dlacn2(n, work[n:], work, iwork, ainvnm, kase, isave)
    61  		if kase == 0 {
    62  			if ainvnm != 0 {
    63  				rcond = (1 / anorm) / ainvnm
    64  			}
    65  			return rcond
    66  		}
    67  		if kase == kase1 {
    68  			scale = impl.Dlatrs(uplo, blas.NoTrans, diag, normin, n, a, lda, work, work[2*n:])
    69  		} else {
    70  			scale = impl.Dlatrs(uplo, blas.Trans, diag, normin, n, a, lda, work, work[2*n:])
    71  		}
    72  		normin = true
    73  		if scale != 1 {
    74  			ix := bi.Idamax(n, work, 1)
    75  			xnorm := math.Abs(work[ix])
    76  			if scale == 0 || scale < xnorm*smlnum {
    77  				return rcond
    78  			}
    79  			impl.Drscl(n, scale, work, 1)
    80  		}
    81  	}
    82  }