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