github.com/gopherd/gonum@v0.0.4/lapack/testlapack/dpstf2.go (about)

     1  // Copyright ©2021 The Gonum Authors. All rights reserved.
     2  // Use of this source code is governed by a BSD-style
     3  // license that can be found in the LICENSE file.
     4  
     5  package testlapack
     6  
     7  import (
     8  	"fmt"
     9  	"math"
    10  	"testing"
    11  
    12  	"math/rand"
    13  
    14  	"github.com/gopherd/gonum/blas"
    15  	"github.com/gopherd/gonum/blas/blas64"
    16  )
    17  
    18  type Dpstf2er interface {
    19  	Dpstf2(uplo blas.Uplo, n int, a []float64, lda int, piv []int, tol float64, work []float64) (rank int, ok bool)
    20  }
    21  
    22  func Dpstf2Test(t *testing.T, impl Dpstf2er) {
    23  	rnd := rand.New(rand.NewSource(1))
    24  	for _, uplo := range []blas.Uplo{blas.Upper, blas.Lower} {
    25  		t.Run(uploToString(uplo), func(t *testing.T) {
    26  			for _, n := range []int{0, 1, 2, 3, 4, 5, 10, 20, 50} {
    27  				for _, lda := range []int{max(1, n), n + 5} {
    28  					for _, rank := range []int{int(0.7 * float64(n)), n} {
    29  						dpstf2Test(t, impl, rnd, uplo, n, lda, rank)
    30  					}
    31  				}
    32  			}
    33  		})
    34  	}
    35  }
    36  
    37  func dpstf2Test(t *testing.T, impl Dpstf2er, rnd *rand.Rand, uplo blas.Uplo, n, lda, rankWant int) {
    38  	const tol = 1e-14
    39  
    40  	name := fmt.Sprintf("n=%v,lda=%v", n, lda)
    41  	bi := blas64.Implementation()
    42  
    43  	// Generate a random, symmetric A with the given rank by applying rankWant
    44  	// rank-1 updates to the zero matrix.
    45  	a := make([]float64, n*lda)
    46  	for i := 0; i < rankWant; i++ {
    47  		x := randomSlice(n, rnd)
    48  		bi.Dsyr(uplo, n, 1, x, 1, a, lda)
    49  	}
    50  
    51  	// Make a copy of A for storing the factorization.
    52  	aFac := make([]float64, len(a))
    53  	copy(aFac, a)
    54  
    55  	// Allocate a slice for pivots and fill it with invalid index values.
    56  	piv := make([]int, n)
    57  	for i := range piv {
    58  		piv[i] = -1
    59  	}
    60  
    61  	// Allocate the work slice.
    62  	work := make([]float64, 2*n)
    63  
    64  	// Call Dpstf2 to Compute the Cholesky factorization with complete pivoting.
    65  	rank, ok := impl.Dpstf2(uplo, n, aFac, lda, piv, -1, work)
    66  
    67  	if ok != (rank == n) {
    68  		t.Errorf("%v: unexpected ok; got %v, want %v", name, ok, rank == n)
    69  	}
    70  	if rank != rankWant {
    71  		t.Errorf("%v: unexpected rank; got %v, want %v", name, rank, rankWant)
    72  	}
    73  
    74  	if n == 0 {
    75  		return
    76  	}
    77  
    78  	// Check that the residual |P*Uᵀ*U*Pᵀ - A| / n or |P*L*Lᵀ*Pᵀ - A| / n is
    79  	// sufficiently small.
    80  	resid := residualDpstrf(uplo, n, a, aFac, lda, rank, piv)
    81  	if resid > tol || math.IsNaN(resid) {
    82  		t.Errorf("%v: residual too large; got %v, want<=%v", name, resid, tol)
    83  	}
    84  }