gonum.org/v1/gonum@v0.15.1-0.20240517103525-f853624cb1bb/lapack/gonum/dgecon.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  	"gonum.org/v1/gonum/blas"
    11  	"gonum.org/v1/gonum/blas/blas64"
    12  	"gonum.org/v1/gonum/lapack"
    13  )
    14  
    15  // Dgecon estimates and returns the reciprocal of the condition number of the
    16  // n×n matrix A, in either the 1-norm or the ∞-norm, using the LU factorization
    17  // computed by Dgetrf.
    18  //
    19  // An estimate is obtained for norm(A⁻¹), and the reciprocal of the condition
    20  // number rcond is computed as
    21  //
    22  //	rcond 1 / ( norm(A) * norm(A⁻¹) ).
    23  //
    24  // If n is zero, rcond is always 1.
    25  //
    26  // anorm is the 1-norm or the ∞-norm of the original matrix A. anorm must be
    27  // non-negative, otherwise Dgecon will panic. If anorm is 0 or infinity, Dgecon
    28  // returns 0. If anorm is NaN, Dgecon returns NaN.
    29  //
    30  // work must have length at least 4*n and iwork must have length at least n,
    31  // otherwise Dgecon will panic.
    32  func (impl Implementation) Dgecon(norm lapack.MatrixNorm, n int, a []float64, lda int, anorm float64, work []float64, iwork []int) float64 {
    33  	switch {
    34  	case norm != lapack.MaxColumnSum && norm != lapack.MaxRowSum:
    35  		panic(badNorm)
    36  	case n < 0:
    37  		panic(nLT0)
    38  	case lda < max(1, n):
    39  		panic(badLdA)
    40  	case anorm < 0:
    41  		panic(negANorm)
    42  	}
    43  
    44  	// Quick return if possible.
    45  	if n == 0 {
    46  		return 1
    47  	}
    48  
    49  	switch {
    50  	case len(a) < (n-1)*lda+n:
    51  		panic(shortA)
    52  	case len(work) < 4*n:
    53  		panic(shortWork)
    54  	case len(iwork) < n:
    55  		panic(shortIWork)
    56  	}
    57  
    58  	// Quick return if possible.
    59  	switch {
    60  	case anorm == 0:
    61  		return 0
    62  	case math.IsNaN(anorm):
    63  		// Propagate NaN.
    64  		return anorm
    65  	case math.IsInf(anorm, 1):
    66  		return 0
    67  	}
    68  
    69  	bi := blas64.Implementation()
    70  	var rcond, ainvnm float64
    71  	var kase int
    72  	var normin bool
    73  	isave := new([3]int)
    74  	onenrm := norm == lapack.MaxColumnSum
    75  	smlnum := dlamchS
    76  	kase1 := 2
    77  	if onenrm {
    78  		kase1 = 1
    79  	}
    80  	for {
    81  		ainvnm, kase = impl.Dlacn2(n, work[n:], work, iwork, ainvnm, kase, isave)
    82  		if kase == 0 {
    83  			if ainvnm != 0 {
    84  				rcond = (1 / ainvnm) / anorm
    85  			}
    86  			return rcond
    87  		}
    88  		var sl, su float64
    89  		if kase == kase1 {
    90  			sl = impl.Dlatrs(blas.Lower, blas.NoTrans, blas.Unit, normin, n, a, lda, work, work[2*n:])
    91  			su = impl.Dlatrs(blas.Upper, blas.NoTrans, blas.NonUnit, normin, n, a, lda, work, work[3*n:])
    92  		} else {
    93  			su = impl.Dlatrs(blas.Upper, blas.Trans, blas.NonUnit, normin, n, a, lda, work, work[3*n:])
    94  			sl = impl.Dlatrs(blas.Lower, blas.Trans, blas.Unit, normin, n, a, lda, work, work[2*n:])
    95  		}
    96  		scale := sl * su
    97  		normin = true
    98  		if scale != 1 {
    99  			ix := bi.Idamax(n, work, 1)
   100  			if scale == 0 || scale < math.Abs(work[ix])*smlnum {
   101  				return rcond
   102  			}
   103  			impl.Drscl(n, scale, work, 1)
   104  		}
   105  	}
   106  }