Repository: gonum/gonum
Pull request #1726 (2021-09-27 17:35:05)

lapack/gonum: add Dpstf2 and its test
Please take a look.

Updates #1717

<!--
Checklist:

- API changes have been discussed
- code is goformated correctly (goimports)
- packages with generated code have had code generation run
- tests pass locally
- linked to relevant issues

Please make sure your commit message summary line and pull request
title match the Go convention; a one-line summary of the change,
prefixed by the primary affected package that should complete the
sentence, "This change modifies Gonum to _____."
-->



Patch:
diff --git a/README.md b/README.md
index f7ab18ea..2808820d 100644
--- a/README.md
+++ b/README.md
@@ -18,7 +18,7 @@ go get -u gonum.org/v1/gonum/...
 
 ## Supported Go versions
 
-Gonum supports and tests using the gc compiler on the [two most recent Go releases](https://github.com/gonum/gonum/blob/master/.github/workflows/ci.yml#L14-15) on Linux (386, amd64 and arm64), macOS and Windows (both on amd64).
+Gonum supports and tests using the gc compiler on the [two most recent Go releases](https://github.com/gonum/gonum/blob/master/.github/workflows/ci.yml#L14-L15) on Linux (386, amd64 and arm64), macOS and Windows (both on amd64).
 
 ## Release schedule
 
diff --git a/lapack/gonum/dpstf2.go b/lapack/gonum/dpstf2.go
new file mode 100644
index 00000000..565bb454
--- /dev/null
+++ b/lapack/gonum/dpstf2.go
@@ -0,0 +1,190 @@
+// Copyright ©2021 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"
+	"gonum.org/v1/gonum/blas/blas64"
+)
+
+// Dpstf2 computes the Cholesky factorization with complete pivoting of an n×n
+// symmetric positive semidefinite matrix A.
+//
+// The factorization has the form
+//  Pᵀ * A * P = Uᵀ * U ,  if uplo = blas.Upper,
+//  Pᵀ * A * P = L  * Lᵀ,  if uplo = blas.Lower,
+// where U is an upper triangular matrix and L is lower triangular, and P is
+// stored as vector piv.
+//
+// Dpstf2 does not attempt to check that A is positive semidefinite.
+//
+// The length of piv must be n and the length of work must be at least 2*n,
+// otherwise Dpstf2 will panic.
+//
+// Dpstf2 is an internal routine. It is exported for testing purposes.
+func (Implementation) Dpstf2(uplo blas.Uplo, n int, a []float64, lda int, piv []int, tol float64, work []float64) (rank int, ok bool) {
+	switch {
+	case uplo != blas.Upper && uplo != blas.Lower:
+		panic(badUplo)
+	case n < 0:
+		panic(nLT0)
+	case lda < max(1, n):
+		panic(badLdA)
+	}
+
+	// Quick return if possible.
+	if n == 0 {
+		return 0, true
+	}
+
+	switch {
+	case len(a) < (n-1)*lda+n:
+		panic(shortA)
+	case len(piv) != n:
+		panic(badLenPiv)
+	case len(work) < 2*n:
+		panic(shortWork)
+	}
+
+	// Initialize piv.
+	for i := range piv[:n] {
+		piv[i] = i
+	}
+
+	// Compute the first pivot.
+	pvt := 0
+	ajj := a[0]
+	for i := 1; i < n; i++ {
+		aii := a[i*lda+i]
+		if aii > ajj {
+			pvt = i
+			ajj = aii
+		}
+	}
+	if ajj <= 0 || math.IsNaN(ajj) {
+		return 0, false
+	}
+
+	// Compute stopping value if not supplied.
+	dstop := tol
+	if dstop < 0 {
+		dstop = float64(n) * dlamchE * ajj
+	}
+
+	// Set first half of work to zero, holds dot products.
+	dots := work[:n]
+	for i := range dots {
+		dots[i] = 0
+	}
+	work2 := work[n : 2*n]
+
+	bi := blas64.Implementation()
+	if uplo == blas.Upper {
+		// Compute the Cholesky factorization  Pᵀ * A * P = Uᵀ * U.
+		for j := 0; j < n; j++ {
+			// Update dot products and compute possible pivots which are stored
+			// in the second half of work.
+			for i := j; i < n; i++ {
+				if j > 0 {
+					tmp := a[(j-1)*lda+i]
+					dots[i] += tmp * tmp
+				}
+				work2[i] = a[i*lda+i] - dots[i]
+			}
+			if j > 0 {
+				// Find the pivot.
+				pvt = j
+				ajj = work2[pvt]
+				for k := j + 1; k < n; k++ {
+					wk := work2[k]
+					if wk > ajj {
+						pvt = k
+						ajj = wk
+					}
+				}
+				// Test for exit.
+				if ajj <= dstop || math.IsNaN(ajj) {
+					a[j*lda+j] = ajj
+					return j, false
+				}
+			}
+			if j != pvt {
+				// Swap pivot rows and columns.
+				a[pvt*lda+pvt] = a[j*lda+j]
+				bi.Dswap(j, a[j:], lda, a[pvt:], lda)
+				if pvt < n-1 {
+					bi.Dswap(n-pvt-1, a[j*lda+(pvt+1):], 1, a[pvt*lda+(pvt+1):], 1)
+				}
+				bi.Dswap(pvt-j-1, a[j*lda+(j+1):], 1, a[(j+1)*lda+pvt:], lda)
+				// Swap dot products and piv.
+				dots[j], dots[pvt] = dots[pvt], dots[j]
+				piv[j], piv[pvt] = piv[pvt], piv[j]
+			}
+			ajj = math.Sqrt(ajj)
+			a[j*lda+j] = ajj
+			// Compute elements j+1:n of row j.
+			if j < n-1 {
+				bi.Dgemv(blas.Trans, j, n-j-1,
+					-1, a[j+1:], lda, a[j:], lda,
+					1, a[j*lda+j+1:], 1)
+				bi.Dscal(n-j-1, 1/ajj, a[j*lda+j+1:], 1)
+			}
+		}
+	} else {
+		// Compute the Cholesky factorization  Pᵀ * A * P = L * Lᵀ.
+		for j := 0; j < n; j++ {
+			// Update dot products and compute possible pivots which are stored
+			// in the second half of work.
+			for i := j; i < n; i++ {
+				if j > 0 {
+					tmp := a[i*lda+(j-1)]
+					dots[i] += tmp * tmp
+				}
+				work2[i] = a[i*lda+i] - dots[i]
+			}
+			if j > 0 {
+				// Find the pivot.
+				pvt = j
+				ajj = work2[pvt]
+				for k := j + 1; k < n; k++ {
+					wk := work2[k]
+					if wk > ajj {
+						pvt = k
+						ajj = wk
+					}
+				}
+				// Test for exit.
+				if ajj <= dstop || math.IsNaN(ajj) {
+					a[j*lda+j] = ajj
+					return j, false
+				}
+			}
+			if j != pvt {
+				// Swap pivot rows and columns.
+				a[pvt*lda+pvt] = a[j*lda+j]
+				bi.Dswap(j, a[j*lda:], 1, a[pvt*lda:], 1)
+				if pvt < n-1 {
+					bi.Dswap(n-pvt-1, a[(pvt+1)*lda+j:], lda, a[(pvt+1)*lda+pvt:], lda)
+				}
+				bi.Dswap(pvt-j-1, a[(j+1)*lda+j:], lda, a[pvt*lda+(j+1):], 1)
+				// Swap dot products and piv.
+				dots[j], dots[pvt] = dots[pvt], dots[j]
+				piv[j], piv[pvt] = piv[pvt], piv[j]
+			}
+			ajj = math.Sqrt(ajj)
+			a[j*lda+j] = ajj
+			// Compute elements j+1:n of column j.
+			if j < n-1 {
+				bi.Dgemv(blas.NoTrans, n-j-1, j,
+					-1, a[(j+1)*lda:], lda, a[j*lda:], 1,
+					1, a[(j+1)*lda+j:], lda)
+				bi.Dscal(n-j-1, 1/ajj, a[(j+1)*lda+j:], lda)
+			}
+		}
+	}
+	return n, true
+}
diff --git a/lapack/gonum/errors.go b/lapack/gonum/errors.go
index 1429e65d..406a8c10 100644
--- a/lapack/gonum/errors.go
+++ b/lapack/gonum/errors.go
@@ -104,6 +104,7 @@ const (
 	badLenJpiv     = "lapack: bad length of jpiv"
 	badLenJpvt     = "lapack: bad length of jpvt"
 	badLenK        = "lapack: bad length of k"
+	badLenPiv      = "lapack: bad length of piv"
 	badLenSelected = "lapack: bad length of selected"
 	badLenSi       = "lapack: bad length of si"
 	badLenSr       = "lapack: bad length of sr"


Test patch:
diff --git a/lapack/gonum/lapack_test.go b/lapack/gonum/lapack_test.go
index 8d30502e..ab00b6fe 100644
--- a/lapack/gonum/lapack_test.go
+++ b/lapack/gonum/lapack_test.go
@@ -553,6 +553,11 @@ func TestDpotrs(t *testing.T) {
 	testlapack.DpotrsTest(t, impl)
 }
 
+func TestDpstf2(t *testing.T) {
+	t.Parallel()
+	testlapack.Dpstf2Test(t, impl)
+}
+
 func TestDrscl(t *testing.T) {
 	t.Parallel()
 	testlapack.DrsclTest(t, impl)
diff --git a/lapack/testlapack/dpstf2.go b/lapack/testlapack/dpstf2.go
new file mode 100644
index 00000000..f7223688
--- /dev/null
+++ b/lapack/testlapack/dpstf2.go
@@ -0,0 +1,166 @@
+// Copyright ©2021 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/blas"
+	"gonum.org/v1/gonum/blas/blas64"
+	"gonum.org/v1/gonum/lapack"
+)
+
+type Dpstf2er interface {
+	Dpstf2(uplo blas.Uplo, n int, a []float64, lda int, piv []int, tol float64, work []float64) (rank int, ok bool)
+}
+
+func Dpstf2Test(t *testing.T, impl Dpstf2er) {
+	rnd := rand.New(rand.NewSource(1))
+	for _, uplo := range []blas.Uplo{blas.Upper, blas.Lower} {
+		t.Run(uploToString(uplo), func(t *testing.T) {
+			for _, n := range []int{0, 1, 2, 3, 4, 5, 10, 20, 50} {
+				for _, lda := range []int{max(1, n), n + 5} {
+					for _, rank := range []int{int(0.7 * float64(n)), n} {
+						dpstf2Test(t, impl, rnd, uplo, n, lda, rank)
+					}
+				}
+			}
+		})
+	}
+}
+
+func dpstf2Test(t *testing.T, impl Dpstf2er, rnd *rand.Rand, uplo blas.Uplo, n, lda, rankWant int) {
+	const tol = 1e-14
+
+	name := fmt.Sprintf("n=%v,lda=%v", n, lda)
+	bi := blas64.Implementation()
+
+	// Generate a random, symmetric A with the given rank by applying rankWant
+	// rank-1 updates to the zero matrix.
+	a := make([]float64, n*lda)
+	for i := 0; i < rankWant; i++ {
+		x := randomSlice(n, rnd)
+		bi.Dsyr(uplo, n, 1, x, 1, a, lda)
+	}
+
+	// Make a copy of A for storing the factorization.
+	aFac := make([]float64, len(a))
+	copy(aFac, a)
+
+	// Allocate a slice for pivots and fill it with invalid index values.
+	piv := make([]int, n)
+	for i := range piv {
+		piv[i] = -1
+	}
+
+	// Allocate the work slice.
+	work := make([]float64, 2*n)
+
+	// Call Dpstf2 to Compute the Cholesky factorization with complete pivoting.
+	rank, ok := impl.Dpstf2(uplo, n, aFac, lda, piv, -1, work)
+
+	if ok != (rank == n) {
+		t.Errorf("%v: unexpected ok; got %v, want %v", name, ok, rank == n)
+	}
+	if rank != rankWant {
+		t.Errorf("%v: unexpected rank; got %v, want %v", name, rank, rankWant)
+	}
+
+	if n == 0 {
+		return
+	}
+
+	// Reconstruct the symmetric positive semi-definite matrix A from its L or U
+	// factors and the permutation matrix P.
+	perm := zeros(n, n, n)
+	if uplo == blas.Upper {
+		// Change notation.
+		u, ldu := aFac, lda
+		// Zero out last n-rank rows of the factor U.
+		for i := rank; i < n; i++ {
+			for j := i; j < n; j++ {
+				u[i*ldu+j] = 0
+			}
+		}
+		// Extract U to aRec.
+		aRec := zeros(n, n, n)
+		for i := 0; i < n; i++ {
+			for j := i; j < n; j++ {
+				aRec.Data[i*aRec.Stride+j] = u[i*ldu+j]
+			}
+		}
+		// Multiply U by Uᵀ from the left.
+		bi.Dtrmm(blas.Left, blas.Upper, blas.Trans, blas.NonUnit, n, n,
+			1, u, ldu, aRec.Data, aRec.Stride)
+		// Form P * Uᵀ * U * Pᵀ.
+		for i := 0; i < n; i++ {
+			for j := 0; j < n; j++ {
+				if piv[i] > piv[j] {
+					// Don't set the lower triangle.
+					continue
+				}
+				if i <= j {
+					perm.Data[piv[i]*perm.Stride+piv[j]] = aRec.Data[i*aRec.Stride+j]
+				} else {
+					perm.Data[piv[i]*perm.Stride+piv[j]] = aRec.Data[j*aRec.Stride+i]
+				}
+			}
+		}
+		// Compute the difference P*Uᵀ*U*Pᵀ - A.
+		for i := 0; i < n; i++ {
+			for j := i; j < n; j++ {
+				perm.Data[i*perm.Stride+j] -= a[i*lda+j]
+			}
+		}
+	} else {
+		// Change notation.
+		l, ldl := aFac, lda
+		// Zero out last n-rank columns of the factor L.
+		for i := rank; i < n; i++ {
+			for j := rank; j <= i; j++ {
+				l[i*ldl+j] = 0
+			}
+		}
+		// Extract L to aRec.
+		aRec := zeros(n, n, n)
+		for i := 0; i < n; i++ {
+			for j := 0; j <= i; j++ {
+				aRec.Data[i*aRec.Stride+j] = l[i*ldl+j]
+			}
+		}
+		// Multiply L by Lᵀ from the right.
+		bi.Dtrmm(blas.Right, blas.Lower, blas.Trans, blas.NonUnit, n, n,
+			1, l, ldl, aRec.Data, aRec.Stride)
+		// Form P * L * Lᵀ * Pᵀ.
+		for i := 0; i < n; i++ {
+			for j := 0; j < n; j++ {
+				if piv[i] < piv[j] {
+					// Don't set the upper triangle.
+					continue
+				}
+				if i >= j {
+					perm.Data[piv[i]*perm.Stride+piv[j]] = aRec.Data[i*aRec.Stride+j]
+				} else {
+					perm.Data[piv[i]*perm.Stride+piv[j]] = aRec.Data[j*aRec.Stride+i]
+				}
+			}
+		}
+		// Compute the difference P*L*Lᵀ*Pᵀ - A.
+		for i := 0; i < n; i++ {
+			for j := 0; j <= i; j++ {
+				perm.Data[i*perm.Stride+j] -= a[i*lda+j]
+			}
+		}
+	}
+	// Compute |P*Uᵀ*U*Pᵀ - A| / n or |P*L*Lᵀ*Pᵀ - A| / n.
+	resid := dlansy(lapack.MaxColumnSum, uplo, n, perm.Data, perm.Stride) / float64(n)
+	if resid > tol || math.IsNaN(resid) {
+		t.Errorf("%v: residual too large; got %v, want<=%v", name, resid, tol)
+	}
+}