blob: f55a596c347ec001f1074a3ae83c2b22267c8bb1 [file] [log] [blame]
 // Copyright ©2019 The Gonum Authors. All rights reserved. // Use of this source code is governed by a BSD-style // license that can be found in the LICENSE file. package gonum import ( "math" "gonum.org/v1/gonum/blas" "gonum.org/v1/gonum/blas/blas64" ) // Dpbcon returns an estimate of the reciprocal of the condition number (in the // 1-norm) of an n×n symmetric positive definite band matrix using the Cholesky // factorization // A = Uᵀ*U if uplo == blas.Upper // A = L*Lᵀ if uplo == blas.Lower // computed by Dpbtrf. The estimate is obtained for norm(inv(A)), and the // reciprocal of the condition number is computed as // rcond = 1 / (anorm * norm(inv(A))). // // The length of work must be at least 3*n and the length of iwork must be at // least n. func (impl Implementation) Dpbcon(uplo blas.Uplo, n, kd int, ab []float64, ldab int, anorm float64, work []float64, iwork []int) (rcond float64) { switch { case uplo != blas.Upper && uplo != blas.Lower: panic(badUplo) case n < 0: panic(nLT0) case kd < 0: panic(kdLT0) case ldab < kd+1: panic(badLdA) case anorm < 0: panic(badNorm) } // Quick return if possible. if n == 0 { return 1 } switch { case len(ab) < (n-1)*ldab+kd+1: panic(shortAB) case len(work) < 3*n: panic(shortWork) case len(iwork) < n: panic(shortIWork) } // Quick return if possible. if anorm == 0 { return 0 } const smlnum = dlamchS var ( ainvnm float64 kase int isave [3]int normin bool // Denote work slices. x = work[:n] v = work[n : 2*n] cnorm = work[2*n : 3*n] ) // Estimate the 1-norm of the inverse. bi := blas64.Implementation() for { ainvnm, kase = impl.Dlacn2(n, v, x, iwork, ainvnm, kase, &isave) if kase == 0 { break } var op1, op2 blas.Transpose if uplo == blas.Upper { // Multiply x by inv(Uᵀ), op1 = blas.Trans // then by inv(Uᵀ). op2 = blas.NoTrans } else { // Multiply x by inv(L), op1 = blas.NoTrans // then by inv(Lᵀ). op2 = blas.Trans } scaleL := impl.Dlatbs(uplo, op1, blas.NonUnit, normin, n, kd, ab, ldab, x, cnorm) normin = true scaleU := impl.Dlatbs(uplo, op2, blas.NonUnit, normin, n, kd, ab, ldab, x, cnorm) // Multiply x by 1/scale if doing so will not cause overflow. scale := scaleL * scaleU if scale != 1 { ix := bi.Idamax(n, x, 1) if scale < math.Abs(x[ix])*smlnum || scale == 0 { return 0 } impl.Drscl(n, scale, x, 1) } } if ainvnm == 0 { return 0 } // Return the estimate of the reciprocal condition number. return (1 / ainvnm) / anorm }