ztrsm.go

// Copyright ©2019 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 testblas

import (
	"fmt"
	"testing"

	"golang.org/x/exp/rand"

	"github.com/ArkaGPL/gonum/blas"
)

type Ztrsmer interface {
	Ztrsm(side blas.Side, uplo blas.Uplo, transA blas.Transpose, diag blas.Diag, m, n int, alpha complex128, a []complex128, lda int, b []complex128, ldb int)
}

func ZtrsmTest(t *testing.T, impl Ztrsmer) {
	for _, side := range []blas.Side{blas.Left, blas.Right} {
		for _, uplo := range []blas.Uplo{blas.Lower, blas.Upper} {
			for _, trans := range []blas.Transpose{blas.NoTrans, blas.Trans, blas.ConjTrans} {
				for _, diag := range []blas.Diag{blas.Unit, blas.NonUnit} {
					name := sideString(side) + "-" + uploString(uplo) + "-" + transString(trans) + "-" + diagString(diag)
					t.Run(name, func(t *testing.T) {
						for _, m := range []int{0, 1, 2, 3, 4, 5} {
							for _, n := range []int{0, 1, 2, 3, 4, 5} {
								ztrsmTest(t, impl, side, uplo, trans, diag, m, n)
							}
						}
					})
				}
			}
		}
	}
}

func ztrsmTest(t *testing.T, impl Ztrsmer, side blas.Side, uplo blas.Uplo, trans blas.Transpose, diag blas.Diag, m, n int) {
	const tol = 1e-13

	rnd := rand.New(rand.NewSource(1))

	nA := m
	if side == blas.Right {
		nA = n
	}
	for _, lda := range []int{max(1, nA), nA + 2} {
		for _, ldb := range []int{max(1, n), n + 3} {
			for _, alpha := range []complex128{0, 1, complex(0.7, -0.9)} {
				// Allocate the matrix A and fill it with random numbers.
				a := make([]complex128, nA*lda)
				for i := range a {
					a[i] = rndComplex128(rnd)
				}
				// Set some elements of A to 0 and 1 to cover special cases in Ztrsm.
				if nA > 2 {
					if uplo == blas.Upper {
						a[nA-2] = 1
						a[nA-1] = 0
					} else {
						a[(nA-2)*lda] = 1
						a[(nA-1)*lda] = 0
					}
				}
				// Create a copy of A for checking that Ztrsm
				// does not modify its triangle opposite to uplo.
				aCopy := make([]complex128, len(a))
				copy(aCopy, a)
				// Create a dense representation of A for
				// computing the right-hand side matrix using zmm.
				aTri := make([]complex128, len(a))
				copy(aTri, a)
				if uplo == blas.Upper {
					for i := 0; i < nA; i++ {
						// Zero out the lower triangle.
						for j := 0; j < i; j++ {
							aTri[i*lda+j] = 0
						}
						if diag == blas.Unit {
							aTri[i*lda+i] = 1
						}
					}
				} else {
					for i := 0; i < nA; i++ {
						if diag == blas.Unit {
							aTri[i*lda+i] = 1
						}
						// Zero out the upper triangle.
						for j := i + 1; j < nA; j++ {
							aTri[i*lda+j] = 0
						}
					}
				}

				// Allocate the right-hand side matrix B and fill it with random numbers.
				b := make([]complex128, m*ldb)
				for i := range b {
					b[i] = rndComplex128(rnd)
				}
				// Set some elements of B to 0 to cover special cases in Ztrsm.
				if m > 1 && n > 1 {
					b[0] = 0
					b[(m-1)*ldb+n-1] = 0
				}
				bCopy := make([]complex128, len(b))
				copy(bCopy, b)

				// Compute the solution matrix X using Ztrsm.
				// X is overwritten on B.
				impl.Ztrsm(side, uplo, trans, diag, m, n, alpha, a, lda, b, ldb)
				x := b

				prefix := fmt.Sprintf("m=%v,n=%v,lda=%v,ldb=%v,alpha=%v", m, n, lda, ldb, alpha)

				if !zsame(a, aCopy) {
					t.Errorf("%v: unexpected modification of A", prefix)
					continue
				}

				// Compute the left-hand side matrix of op(A)*X=alpha*B or X*op(A)=alpha*B
				// using an internal Zgemm implementation.
				var lhs []complex128
				if side == blas.Left {
					lhs = zmm(trans, blas.NoTrans, m, n, m, 1, aTri, lda, x, ldb, 0, b, ldb)
				} else {
					lhs = zmm(blas.NoTrans, trans, m, n, n, 1, x, ldb, aTri, lda, 0, b, ldb)
				}

				// Compute the right-hand side matrix alpha*B.
				rhs := bCopy
				for i := 0; i < m; i++ {
					for j := 0; j < n; j++ {
						rhs[i*ldb+j] *= alpha
					}
				}

				if !zEqualApprox(lhs, rhs, tol) {
					t.Errorf("%v: unexpected result", prefix)
				}
			}
		}
	}
}