1// Copyright ©2019 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
5package gonum
6
7import (
8	"math"
9
10	"gonum.org/v1/gonum/blas"
11	"gonum.org/v1/gonum/blas/blas64"
12)
13
14// Dpbcon returns an estimate of the reciprocal of the condition number (in the
15// 1-norm) of an n×n symmetric positive definite band matrix using the Cholesky
16// factorization
17//  A = Uᵀ*U  if uplo == blas.Upper
18//  A = L*Lᵀ  if uplo == blas.Lower
19// computed by Dpbtrf. The estimate is obtained for norm(inv(A)), and the
20// reciprocal of the condition number is computed as
21//  rcond = 1 / (anorm * norm(inv(A))).
22//
23// The length of work must be at least 3*n and the length of iwork must be at
24// least n.
25func (impl Implementation) Dpbcon(uplo blas.Uplo, n, kd int, ab []float64, ldab int, anorm float64, work []float64, iwork []int) (rcond float64) {
26	switch {
27	case uplo != blas.Upper && uplo != blas.Lower:
28		panic(badUplo)
29	case n < 0:
30		panic(nLT0)
31	case kd < 0:
32		panic(kdLT0)
33	case ldab < kd+1:
34		panic(badLdA)
35	case anorm < 0:
36		panic(badNorm)
37	}
38
39	// Quick return if possible.
40	if n == 0 {
41		return 1
42	}
43
44	switch {
45	case len(ab) < (n-1)*ldab+kd+1:
46		panic(shortAB)
47	case len(work) < 3*n:
48		panic(shortWork)
49	case len(iwork) < n:
50		panic(shortIWork)
51	}
52
53	// Quick return if possible.
54	if anorm == 0 {
55		return 0
56	}
57
58	const smlnum = dlamchS
59
60	var (
61		ainvnm float64
62		kase   int
63		isave  [3]int
64		normin bool
65
66		// Denote work slices.
67		x     = work[:n]
68		v     = work[n : 2*n]
69		cnorm = work[2*n : 3*n]
70	)
71	// Estimate the 1-norm of the inverse.
72	bi := blas64.Implementation()
73	for {
74		ainvnm, kase = impl.Dlacn2(n, v, x, iwork, ainvnm, kase, &isave)
75		if kase == 0 {
76			break
77		}
78		var op1, op2 blas.Transpose
79		if uplo == blas.Upper {
80			// Multiply x by inv(Uᵀ),
81			op1 = blas.Trans
82			// then by inv(Uᵀ).
83			op2 = blas.NoTrans
84		} else {
85			// Multiply x by inv(L),
86			op1 = blas.NoTrans
87			// then by inv(Lᵀ).
88			op2 = blas.Trans
89		}
90		scaleL := impl.Dlatbs(uplo, op1, blas.NonUnit, normin, n, kd, ab, ldab, x, cnorm)
91		normin = true
92		scaleU := impl.Dlatbs(uplo, op2, blas.NonUnit, normin, n, kd, ab, ldab, x, cnorm)
93		// Multiply x by 1/scale if doing so will not cause overflow.
94		scale := scaleL * scaleU
95		if scale != 1 {
96			ix := bi.Idamax(n, x, 1)
97			if scale < math.Abs(x[ix])*smlnum || scale == 0 {
98				return 0
99			}
100			impl.Drscl(n, scale, x, 1)
101		}
102	}
103	if ainvnm == 0 {
104		return 0
105	}
106	// Return the estimate of the reciprocal condition number.
107	return (1 / ainvnm) / anorm
108}
109