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

     1  // Copyright ©2018 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  	"testing"
    10  
    11  	"math/rand"
    12  
    13  	"github.com/gopherd/gonum/blas"
    14  	"github.com/gopherd/gonum/blas/blas64"
    15  )
    16  
    17  type Dpotrser interface {
    18  	Dpotrs(uplo blas.Uplo, n, nrhs int, a []float64, lda int, b []float64, ldb int)
    19  
    20  	Dpotrf(uplo blas.Uplo, n int, a []float64, lda int) bool
    21  }
    22  
    23  func DpotrsTest(t *testing.T, impl Dpotrser) {
    24  	const tol = 1e-14
    25  
    26  	rnd := rand.New(rand.NewSource(1))
    27  	bi := blas64.Implementation()
    28  
    29  	for _, uplo := range []blas.Uplo{blas.Upper, blas.Lower} {
    30  		for _, n := range []int{1, 2, 5} {
    31  			for _, nrhs := range []int{1, 2, 5, 10} {
    32  				for _, ld := range []struct{ a, b int }{
    33  					{n, nrhs},
    34  					{n + 7, nrhs},
    35  					{n, nrhs + 3},
    36  					{n + 7, nrhs + 3},
    37  				} {
    38  					// Construct a random SPD matrix A by first making a symmetric matrix
    39  					// and then ensuring that it is diagonally dominant.
    40  					a := nanGeneral(n, n, ld.a)
    41  					for i := 0; i < n; i++ {
    42  						for j := i; j < n; j++ {
    43  							v := rnd.Float64()
    44  							a.Data[i*a.Stride+j] = v
    45  							a.Data[j*a.Stride+i] = v
    46  						}
    47  					}
    48  					for i := 0; i < n; i++ {
    49  						a.Data[i*a.Stride+i] += float64(n)
    50  					}
    51  
    52  					// Generate a random solution X.
    53  					want := nanGeneral(n, nrhs, ld.b)
    54  					for i := 0; i < n; i++ {
    55  						for j := 0; j < nrhs; j++ {
    56  							want.Data[i*want.Stride+j] = rnd.NormFloat64()
    57  						}
    58  					}
    59  
    60  					// Compute the right-hand side matrix as A * X.
    61  					b := nanGeneral(n, nrhs, ld.b)
    62  					bi.Dgemm(blas.NoTrans, blas.NoTrans, n, nrhs, n, 1, a.Data, a.Stride, want.Data, want.Stride, 0, b.Data, b.Stride)
    63  
    64  					// Compute the Cholesky decomposition of A.
    65  					ok := impl.Dpotrf(uplo, n, a.Data, a.Stride)
    66  					if !ok {
    67  						panic("bad test")
    68  					}
    69  
    70  					aCopy := cloneGeneral(a)
    71  
    72  					// Solve A * X = B.
    73  					impl.Dpotrs(uplo, n, nrhs, a.Data, a.Stride, b.Data, b.Stride)
    74  
    75  					name := fmt.Sprintf("uplo=%v,n=%v,nrhs=%v,lda=%v,ldb=%v", uplo, n, nrhs, a.Stride, b.Stride)
    76  
    77  					if !generalOutsideAllNaN(a) {
    78  						t.Errorf("%v: out-of-range modification of A", name)
    79  					}
    80  					if !equalApproxGeneral(a, aCopy, 0) {
    81  						t.Errorf("%v: unexpected modification of A", name)
    82  					}
    83  					if !generalOutsideAllNaN(b) {
    84  						t.Errorf("%v: out-of-range modification of B", name)
    85  					}
    86  					if !equalApproxGeneral(b, want, tol) {
    87  						t.Errorf("%v: unexpected result\ngot  %v\nwant %v", name, b, want)
    88  					}
    89  				}
    90  			}
    91  		}
    92  	}
    93  }