Merge 2ec9487 into 16abd5a

lapack/netlib/lapack.go

@@ -532,6 +532,134 @@ func (impl Implementation) Dlange(norm lapack.MatrixNorm, m, n int, a []float64,
 	return lapacke.Dlange(byte(norm), m, n, a, lda, work)
 }
 
+// Dlansb returns the given norm of an n×n symmetric band matrix with kd
+// super-diagonals.
+//
+// When norm is lapack.MaxColumnSum or lapack.MaxRowSum, the length of work must
+// be at least n.
+func (impl Implementation) Dlansb(norm lapack.MatrixNorm, uplo blas.Uplo, n, kd int, ab []float64, ldab int, work []float64) float64 {
+	// This is a pure Go implementation of Dlansb copied verbatim from the gonum
+	// repository. Both implementations should be kept in sync. It must exist
+	// here because LAPACKE doesn't provide it but Implementation must implement
+	// lapack.Float64.
+
+	switch {
+	case norm != lapack.MaxAbs && norm != lapack.MaxRowSum && norm != lapack.MaxColumnSum && norm != lapack.Frobenius:
+		panic(badNorm)
+	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)
+	}
+
+	// Quick return if possible.
+	if n == 0 {
+		return 0
+	}
+
+	switch {
+	case len(ab) < (n-1)*ldab+kd+1:
+		panic(shortAB)
+	case len(work) < n && (norm == lapack.MaxColumnSum || norm == lapack.MaxRowSum):
+		panic(shortWork)
+	}
+
+	var value float64
+	switch norm {
+	case lapack.MaxAbs:
+		if uplo == blas.Upper {
+			for i := 0; i < n; i++ {
+				for j := 0; j < min(n-i, kd+1); j++ {
+					aij := math.Abs(ab[i*ldab+j])
+					if aij > value || math.IsNaN(aij) {
+						value = aij
+					}
+				}
+			}
+		} else {
+			for i := 0; i < n; i++ {
+				for j := max(0, kd-i); j < kd+1; j++ {
+					aij := math.Abs(ab[i*ldab+j])
+					if aij > value || math.IsNaN(aij) {
+						value = aij
+					}
+				}
+			}
+		}
+	case lapack.MaxColumnSum, lapack.MaxRowSum:
+		work = work[:n]
+		var sum float64
+		if uplo == blas.Upper {
+			for i := range work {
+				work[i] = 0
+			}
+			for i := 0; i < n; i++ {
+				sum := work[i] + math.Abs(ab[i*ldab])
+				for j := i + 1; j < min(i+kd+1, n); j++ {
+					aij := math.Abs(ab[i*ldab+j-i])
+					sum += aij
+					work[j] += aij
+				}
+				if sum > value || math.IsNaN(sum) {
+					value = sum
+				}
+			}
+		} else {
+			for i := 0; i < n; i++ {
+				sum = 0
+				for j := max(0, i-kd); j < i; j++ {
+					aij := math.Abs(ab[i*ldab+kd+j-i])
+					sum += aij
+					work[j] += aij
+				}
+				work[i] = sum + math.Abs(ab[i*ldab+kd])
+			}
+			for _, sum := range work {
+				if sum > value || math.IsNaN(sum) {
+					value = sum
+				}
+			}
+		}
+	case lapack.Frobenius:
+		scale := 0.0
+		ssq := 1.0
+		if uplo == blas.Upper {
+			if kd > 0 {
+				// Sum off-diagonals.
+				for i := 0; i < n-1; i++ {
+					ilen := min(n-i-1, kd)
+					rowscale, rowssq := impl.Dlassq(ilen, ab[i*ldab+1:], 1, 0, 1)
+					scale, ssq = impl.Dcombssq(scale, ssq, rowscale, rowssq)
+				}
+				ssq *= 2
+			}
+			// Sum diagonal.
+			dscale, dssq := impl.Dlassq(n, ab, ldab, 0, 1)
+			scale, ssq = impl.Dcombssq(scale, ssq, dscale, dssq)
+		} else {
+			if kd > 0 {
+				// Sum off-diagonals.
+				for i := 1; i < n; i++ {
+					ilen := min(i, kd)
+					rowscale, rowssq := impl.Dlassq(ilen, ab[i*ldab+kd-ilen:], 1, 0, 1)
+					scale, ssq = impl.Dcombssq(scale, ssq, rowscale, rowssq)
+				}
+				ssq *= 2
+			}
+			// Sum diagonal.
+			dscale, dssq := impl.Dlassq(n, ab[kd:], ldab, 0, 1)
+			scale, ssq = impl.Dcombssq(scale, ssq, dscale, dssq)
+		}
+		value = scale * math.Sqrt(ssq)
+	}
+
+	return value
+}
+
 // Dlansy computes the specified norm of an n×n symmetric matrix. If
 // norm == lapack.MaxColumnSum or norm == lapackMaxRowSum work must have length
 // at least n, otherwise work is unused.
@@ -794,6 +922,54 @@ func (impl Implementation) Dlaswp(n int, a []float64, lda, k1, k2 int, ipiv []in
 	lapacke.Dlaswp(n, a, lda, k1+1, k2+1, ipiv32, incX)
 }
 
+// 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)
+	}
+
+	_ldab := n
+	_ab := make([]float64, (kd+1)*_ldab)
+	convDpbToLapacke(uplo, n, kd, ab, ldab, _ab, _ldab)
+	_rcond := []float64{0}
+	_iwork := make([]int32, n)
+	lapacke.Dpbcon(byte(uplo), n, kd, _ab, _ldab, anorm, _rcond, work, _iwork)
+	return _rcond[0]
+}
+
 // Dpbtrf computes the Cholesky factorization of an n×n symmetric positive
 // definite band matrix
 //  A = U^T * U  if uplo == blas.Upper
@@ -3882,6 +4058,69 @@ func (impl Implementation) Dtgsja(jobU, jobV, jobQ lapack.GSVDJob, m, p, n, k, l
 	return int(ncycle[0]), ok
 }
 
+// Dlassq updates a sum of squares in scaled form. The input parameters scale and
+// sumsq represent the current scale and total sum of squares. These values are
+// updated with the information in the first n elements of the vector specified
+// by x and incX.
+//
+// Dlassq is an internal routine. It is exported for testing purposes.
+func (impl Implementation) Dlassq(n int, x []float64, incx int, scale float64, sumsq float64) (scl, smsq float64) {
+	// This is a pure Go implementation of Dlassq copied verbatim from the gonum
+	// repository. Both implementations should be kept in sync. It must exist
+	// here because it is used by other pure Go functions that LAPACKE doesn't
+	// provide.
+
+	switch {
+	case n < 0:
+		panic(nLT0)
+	case incx <= 0:
+		panic(badIncX)
+	case len(x) < 1+(n-1)*incx:
+		panic(shortX)
+	}
+
+	if n == 0 {
+		return scale, sumsq
+	}
+
+	for ix := 0; ix <= (n-1)*incx; ix += incx {
+		absxi := math.Abs(x[ix])
+		if absxi > 0 || math.IsNaN(absxi) {
+			if scale < absxi {
+				sumsq = 1 + sumsq*(scale/absxi)*(scale/absxi)
+				scale = absxi
+			} else {
+				sumsq += (absxi / scale) * (absxi / scale)
+			}
+		}
+	}
+	return scale, sumsq
+}
+
+// Dcombssq adds two scaled sum-of-squares quantities, V := V1 + V2,
+//  V_scale^2 * V_ssq := V1_scale^2 * V1_ssq + V2_scale^2 * V2_ssq
+// and returns the result V.
+//
+// Dcombssq is an internal routine. It is exported for testing purposes.
+func (Implementation) Dcombssq(scale1, ssq1, scale2, ssq2 float64) (scale, ssq float64) {
+	// This is a pure Go implementation of Dcombssq copied verbatim from the
+	// gonum repository. Both implementations should be kept in sync. It must
+	// exist here because it is used by other pure Go functions that LAPACKE
+	// doesn't provide.
+
+	if scale1 >= scale2 {
+		if scale1 != 0 {
+			return scale1, ssq1 + (scale2/scale1)*(scale2/scale1)*ssq2
+		}
+		// Both scales are zero.
+		if math.IsNaN(ssq1) || math.IsNaN(ssq2) {
+			return 0, math.NaN()
+		}
+		return 0, 0
+	}
+	return scale2, ssq2 + (scale1/scale2)*(scale1/scale2)*ssq1
+}
+
 func min(m, n int) int {
 	if m < n {
 		return m

lapack/netlib/lapack_test.go

@@ -64,6 +64,18 @@ func TestDlange(t *testing.T) {
 	testlapack.DlangeTest(t, impl)
 }
 
+func TestDlansb(t *testing.T) {
+	testlapack.DlansbTest(t, impl)
+}
+
+func TestDlansy(t *testing.T) {
+	testlapack.DlansyTest(t, impl)
+}
+
+func TestDlantr(t *testing.T) {
+	testlapack.DlantrTest(t, impl)
+}
+
 func TestDlarfb(t *testing.T) {
 	testlapack.DlarfbTest(t, impl)
 }
@@ -76,10 +88,6 @@ func TestDlarft(t *testing.T) {
 	testlapack.DlarftTest(t, impl)
 }
 
-func TestDlantr(t *testing.T) {
-	testlapack.DlantrTest(t, impl)
-}
-
 func TestDlapmt(t *testing.T) {
 	testlapack.DlapmtTest(t, impl)
 }
@@ -104,6 +112,10 @@ func TestDlaswp(t *testing.T) {
 	testlapack.DlaswpTest(t, impl)
 }
 
+func TestDpbcon(t *testing.T) {
+	testlapack.DpbconTest(t, impl)
+}
+
 func TestDpbtrf(t *testing.T) {
 	testlapack.DpbtrfTest(t, impl)
 }
@@ -288,3 +300,11 @@ func TestDtrcon(t *testing.T) {
 func TestDtrtri(t *testing.T) {
 	testlapack.DtrtriTest(t, impl)
 }
+
+func TestDlassq(t *testing.T) {
+	testlapack.DlassqTest(t, impl)
+}
+
+func TestDcombssq(t *testing.T) {
+	testlapack.DcombssqTest(t, impl)
+}