From fa306f215a0186dfb3c1cad95190d7c9a249eb15 Mon Sep 17 00:00:00 2001
From: Vladimir Chalupecky <vladimir@chalupecky.dev>
Date: Mon, 13 Nov 2023 09:45:04 +0100
Subject: [PATCH] lapack/gonum: handle NaN and Inf input to Dgecon

See https://github.com/Reference-LAPACK/lapack/pull/926
---
 lapack/gonum/dgecon.go      | 30 ++++++++++++++++++------------
 lapack/testlapack/dgecon.go | 21 +++++++++++++++++++++
 2 files changed, 39 insertions(+), 12 deletions(-)

diff --git a/lapack/gonum/dgecon.go b/lapack/gonum/dgecon.go
index 618d7fdb..1d046441 100644
--- a/lapack/gonum/dgecon.go
+++ b/lapack/gonum/dgecon.go
@@ -12,14 +12,20 @@ import (
 	"gonum.org/v1/gonum/lapack"
 )
 
-// Dgecon estimates the reciprocal of the condition number of the n×n matrix A
-// given the LU decomposition of the matrix. The condition number computed may
-// be based on the 1-norm or the ∞-norm.
+// Dgecon estimates and returns the reciprocal of the condition number of the
+// n×n matrix A, in either the 1-norm or the ∞-norm, using the LU factorization
+// computed by Dgetrf.
 //
-// The slice a contains the result of the LU decomposition of A as computed by Dgetrf.
+// An estimate is obtained for norm(A⁻¹), and the reciprocal of the condition
+// number rcond is computed as
 //
-// anorm is the corresponding 1-norm or ∞-norm of the original matrix A. anorm
-// must be non-negative.
+//	rcond 1 / ( norm(A) * norm(A⁻¹) ).
+//
+// If n is zero, rcond is always 1.
+//
+// anorm is the 1-norm or the ∞-norm of the original matrix A. anorm must be
+// non-negative, otherwise Dgecon will panic. If anorm is 0 or infinity, Dgecon
+// returns 0. If anorm is NaN, Dgecon returns NaN.
 //
 // work must have length at least 4*n and iwork must have length at least n,
 // otherwise Dgecon will panic.
@@ -50,14 +56,14 @@ func (impl Implementation) Dgecon(norm lapack.MatrixNorm, n int, a []float64, ld
 	}
 
 	// Quick return if possible.
-	if anorm == 0 {
+	switch {
+	case anorm == 0:
 		return 0
-	}
-	if math.IsNaN(anorm) {
-		// The reference implementation treats the NaN anorm as invalid and
-		// returns an error code. Our error handling is to panic which seems too
-		// harsh for a runtime condition, so we just propagate the NaN instead.
+	case math.IsNaN(anorm):
+		// Propagate NaN.
 		return anorm
+	case math.IsInf(anorm, 1):
+		return 0
 	}
 
 	bi := blas64.Implementation()
diff --git a/lapack/testlapack/dgecon.go b/lapack/testlapack/dgecon.go
index c5007578..a6679d3a 100644
--- a/lapack/testlapack/dgecon.go
+++ b/lapack/testlapack/dgecon.go
@@ -6,6 +6,7 @@ package testlapack
 
 import (
 	"fmt"
+	"math"
 	"testing"
 
 	"golang.org/x/exp/rand"
@@ -87,5 +88,25 @@ func dgeconTest(t *testing.T, impl Dgeconer, rnd *rand.Rand, n, lda int) {
 			t.Errorf("%v: unexpected value of rcond; got=%v, want=%v (ratio=%v)",
 				name, rcondGot, rcondWant, ratio)
 		}
+
+		// Check for corner-case values of anorm.
+		for _, anorm := range []float64{0, math.Inf(1), math.NaN()} {
+			rcondGot = impl.Dgecon(norm, n, aFac, lda, anorm, work, iwork)
+			if n == 0 {
+				if rcondGot != 1 {
+					t.Errorf("%v: unexpected rcond when anorm=%v: got=%v, want=1", name, anorm, rcondGot)
+				}
+				continue
+			}
+			if math.IsNaN(anorm) {
+				if !math.IsNaN(rcondGot) {
+					t.Errorf("%v: NaN not propagated when anorm=NaN: got=%v", name, rcondGot)
+				}
+				continue
+			}
+			if rcondGot != 0 {
+				t.Errorf("%v: unexpected rcond when anorm=%v: got=%v, want=0", name, anorm, rcondGot)
+			}
+		}
 	}
 }