lapack/gonum: propagate NaN from Dgecon

lapack/gonum/dgecon.go

@@ -18,11 +18,11 @@ import (
 //
 // The slice a contains the result of the LU decomposition of A as computed by Dgetrf.
 //
-// anorm is the corresponding 1-norm or ∞-norm of the original matrix A.
+// anorm is the corresponding 1-norm or ∞-norm of the original matrix A. anorm
+// must be non-negative.
 //
-// work is a temporary data slice of length at least 4*n and Dgecon will panic otherwise.
-//
-// iwork is a temporary data slice of length at least n and Dgecon will panic otherwise.
+// work must have length at least 4*n and iwork must have length at least n,
+// otherwise Dgecon will panic.
 func (impl Implementation) Dgecon(norm lapack.MatrixNorm, n int, a []float64, lda int, anorm float64, work []float64, iwork []int) float64 {
 	switch {
 	case norm != lapack.MaxColumnSum && norm != lapack.MaxRowSum:
@@ -31,6 +31,8 @@ func (impl Implementation) Dgecon(norm lapack.MatrixNorm, n int, a []float64, ld
 		panic(nLT0)
 	case lda < max(1, n):
 		panic(badLdA)
+	case anorm < 0:
+		panic(negANorm)
 	}
 
 	// Quick return if possible.
@@ -51,6 +53,12 @@ func (impl Implementation) Dgecon(norm lapack.MatrixNorm, n int, a []float64, ld
 	if 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.
+		return anorm
+	}
 
 	bi := blas64.Implementation()
 	var rcond, ainvnm float64