From b27ae13fdd03d83a4be8e309a5e8427e159ee35f Mon Sep 17 00:00:00 2001
From: Vladimir Chalupecky <vladimir@chalupecky.dev>
Date: Sat, 25 Nov 2023 09:36:39 +0100
Subject: [PATCH] lapack/gonum: add Dptcon

---
 lapack/gonum/dptcon.go      | 99 +++++++++++++++++++++++++++++++++++++
 lapack/gonum/lapack_test.go |  5 ++
 lapack/testlapack/dptcon.go | 81 ++++++++++++++++++++++++++++++
 3 files changed, 185 insertions(+)
 create mode 100644 lapack/gonum/dptcon.go
 create mode 100644 lapack/testlapack/dptcon.go

diff --git a/lapack/gonum/dptcon.go b/lapack/gonum/dptcon.go
new file mode 100644
index 00000000..cd41e317
--- /dev/null
+++ b/lapack/gonum/dptcon.go
@@ -0,0 +1,99 @@
+// Copyright ©2023 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/blas64"
+)
+
+// Dptcon computes and returns the reciprocal of the condition number (in the
+// 1-norm) of a symmetric positive definite tridiagonal matrix A using the
+// factorization A = L*D*Lᵀ or A = Uᵀ*D*U computed by Dpttrf.
+//
+// The reciprocal of the condition number is computed as
+//
+//	rcond = 1 / (anorm * ‖A⁻¹‖)
+//
+// and ‖A⁻¹‖ is computed by a direct method.
+//
+// d and e contain, respectively, the n diagonal elements of the diagonal matrix
+// D and the (n-1) off-diagonal elements of the unit bidiagonal factor U or L
+// from the factorization of A, as computed by Dpttrf.
+//
+// anorm is the 1-norm of the original matrix A.
+//
+// work must have length n, otherwise Dptcon will panic.
+func (impl Implementation) Dptcon(n int, d, e []float64, anorm float64, work []float64) (rcond float64) {
+	switch {
+	case n < 0:
+		panic(nLT0)
+	case anorm < 0:
+		panic(badNorm)
+	}
+
+	// Quick return if possible.
+	if n == 0 {
+		return 1
+	}
+
+	switch {
+	case len(d) < n:
+		panic(shortD)
+	case len(e) < n-1:
+		panic(shortE)
+	case len(work) < n:
+		panic(shortWork)
+	}
+
+	// Quick return if possible.
+	switch {
+	case anorm == 0:
+		return 0
+	case math.IsNaN(anorm):
+		// Propagate NaN.
+		return anorm
+	case math.IsInf(anorm, 1):
+		return 0
+	}
+
+	// Check that d[0:n] is positive.
+	for _, di := range d[:n] {
+		if di <= 0 {
+			return 0
+		}
+	}
+
+	// Solve M(A) * x = e, where M(A) = (m[i,j]) is given by
+	//
+	// 	m[i,j] =  abs(A[i,j]), i == j,
+	// 	m[i,j] = -abs(A[i,j]), i != j,
+	//
+	// and e = [1,1,...,1]ᵀ. Note M(A) = M(L)*D*M(L)ᵀ.
+
+	// Solve M(L) * b = e.
+	work[0] = 1
+	for i := 1; i < n; i++ {
+		work[i] = 1 + work[i-1]*math.Abs(e[i-1])
+	}
+
+	// Solve D * M(L)ᵀ * x = b.
+	work[n-1] /= d[n-1]
+	for i := n - 2; i >= 0; i-- {
+		work[i] = work[i]/d[i] + work[i+1]*math.Abs(e[i])
+	}
+
+	// Compute ainvnm = max(x[i]), 0<=i<n.
+	bi := blas64.Implementation()
+	ix := bi.Idamax(n, work, 1)
+	ainvnm := math.Abs(work[ix])
+	if ainvnm == 0 {
+		return 0
+	}
+
+	// Compute the reciprocal condition number.
+	return 1 / ainvnm / anorm
+}
diff --git a/lapack/gonum/lapack_test.go b/lapack/gonum/lapack_test.go
index d5b9b1dc..879bdd53 100644
--- a/lapack/gonum/lapack_test.go
+++ b/lapack/gonum/lapack_test.go
@@ -593,6 +593,11 @@ func TestDptsv(t *testing.T) {
 	testlapack.DptsvTest(t, impl)
 }
 
+func TestDptcon(t *testing.T) {
+	t.Parallel()
+	testlapack.DptconTest(t, impl)
+}
+
 func TestDrscl(t *testing.T) {
 	t.Parallel()
 	testlapack.DrsclTest(t, impl)
diff --git a/lapack/testlapack/dptcon.go b/lapack/testlapack/dptcon.go
new file mode 100644
index 00000000..d25ee50a
--- /dev/null
+++ b/lapack/testlapack/dptcon.go
@@ -0,0 +1,81 @@
+// Copyright ©2023 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 testlapack
+
+import (
+	"fmt"
+	"math"
+	"testing"
+
+	"golang.org/x/exp/rand"
+
+	"gonum.org/v1/gonum/floats"
+	"gonum.org/v1/gonum/lapack"
+)
+
+type Dptconer interface {
+	Dptcon(n int, d, e []float64, anorm float64, work []float64) (rcond float64)
+
+	Dpttrf(n int, d, e []float64) (ok bool)
+	Dpttrs(n, nrhs int, d, e []float64, b []float64, ldb int)
+}
+
+func DptconTest(t *testing.T, impl Dptconer) {
+	rnd := rand.New(rand.NewSource(1))
+	for _, n := range []int{0, 1, 2, 3, 4, 5, 10, 20, 50, 51, 52, 53, 54, 100} {
+		dptconTest(t, impl, rnd, n)
+	}
+}
+
+func dptconTest(t *testing.T, impl Dptconer, rnd *rand.Rand, n int) {
+	const tol = 1e-15
+
+	name := fmt.Sprintf("n=%v", n)
+
+	// Generate a random diagonally dominant symmetric tridiagonal matrix A.
+	d, e := newRandomSymTridiag(n, rnd)
+	aNorm := dlanst(lapack.MaxColumnSum, n, d, e)
+
+	// Compute the Cholesky factorization of A.
+	ok := impl.Dpttrf(n, d, e)
+	if !ok {
+		t.Errorf("%v: bad test matrix, Dpttrf failed", name)
+		return
+	}
+
+	// Compute the reciprocal of the condition number of A.
+	dCopy := make([]float64, len(d))
+	copy(dCopy, d)
+	eCopy := make([]float64, len(e))
+	copy(eCopy, e)
+	work := make([]float64, 3*n)
+	rcondGot := impl.Dptcon(n, d, e, aNorm, work)
+
+	// Check that Dptcon didn't modify d and e.
+	if !floats.Equal(d, dCopy) {
+		t.Errorf("%v: unexpected modification of d", name)
+	}
+	if !floats.Equal(e, eCopy) {
+		t.Errorf("%v: unexpected modification of e", name)
+	}
+
+	// Compute the norm of A⁻¹.
+	aInv, lda := make([]float64, n*n), max(1, n)
+	for i := 0; i < n; i++ {
+		aInv[i*lda+i] = 1
+	}
+	impl.Dpttrs(n, n, d, e, aInv, lda)
+	aInvNorm := dlange(lapack.MaxColumnSum, n, n, aInv, lda)
+
+	rcondWant := 1.0
+	if aNorm > 0 && aInvNorm > 0 {
+		rcondWant = 1 / aNorm / aInvNorm
+	}
+
+	diff := math.Abs(rcondGot - rcondWant)
+	if diff > tol {
+		t.Errorf("%v: unexpected value of rcond. got=%v, want=%v (diff=%v)", name, rcondGot, rcondWant, diff)
+	}
+}