vendor/gonum.org/v1/gonum/lapack/gonum/dpocon.go

   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 )
  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         checkMatrix(n, n, a, lda)
  25         if uplo != blas.Upper && uplo != blas.Lower {
  26                 panic(badUplo)
  27         }
  28         if len(work) < 3*n {
  29                 panic(badWork)
  30         }
  31         if len(iwork) < n {
  32                 panic(badWork)
  33         }
  34         var rcond float64
  35         if n == 0 {
  36                 return 1
  37         }
  38         if anorm == 0 {
  39                 return rcond
  40         }
  41
  42         bi := blas64.Implementation()
  43         var ainvnm float64
  44         smlnum := dlamchS
  45         upper := uplo == blas.Upper
  46         var kase int
  47         var normin bool
  48         isave := new([3]int)
  49         var sl, su float64
  50         for {
  51                 ainvnm, kase = impl.Dlacn2(n, work[n:], work, iwork, ainvnm, kase, isave)
  52                 if kase == 0 {
  53                         if ainvnm != 0 {
  54                                 rcond = (1 / ainvnm) / anorm
  55                         }
  56                         return rcond
  57                 }
  58                 if upper {
  59                         sl = impl.Dlatrs(blas.Upper, blas.Trans, blas.NonUnit, normin, n, a, lda, work, work[2*n:])
  60                         normin = true
  61                         su = impl.Dlatrs(blas.Upper, blas.NoTrans, blas.NonUnit, normin, n, a, lda, work, work[2*n:])
  62                 } else {
  63                         sl = impl.Dlatrs(blas.Lower, blas.NoTrans, blas.NonUnit, normin, n, a, lda, work, work[2*n:])
  64                         normin = true
  65                         su = impl.Dlatrs(blas.Lower, blas.Trans, blas.NonUnit, normin, n, a, lda, work, work[2*n:])
  66                 }
  67                 scale := sl * su
  68                 if scale != 1 {
  69                         ix := bi.Idamax(n, work, 1)
  70                         if scale == 0 || scale < math.Abs(work[ix])*smlnum {
  71                                 return rcond
  72                         }
  73                         impl.Drscl(n, scale, work, 1)
  74                 }
  75         }
  76 }