dgesc2.go

// 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/floats"
	"gonum.org/v1/gonum/lapack"
)

type Dgesc2er interface {
	Dgesc2(n int, a []float64, lda int, rhs []float64, ipiv, jpiv []int) (scale float64)

	Dgetc2er
}

func Dgesc2Test(t *testing.T, impl Dgesc2er) {
	rnd := rand.New(rand.NewSource(1))
	for _, n := range []int{0, 1, 2, 3, 4, 5, 10, 20, 50} {
		for _, lda := range []int{n, n + 3} {
			testDgesc2(t, impl, rnd, n, lda, false)
			testDgesc2(t, impl, rnd, n, lda, true)
		}
	}
}

func testDgesc2(t *testing.T, impl Dgesc2er, rnd *rand.Rand, n, lda int, big bool) {
	const tol = 1e-14

	name := fmt.Sprintf("n=%v,lda=%v,big=%v", n, lda, big)

	// Generate random general matrix.
	a := randomGeneral(n, n, max(1, lda), rnd)

	// Generate a random right hand side vector.
	b := randomGeneral(n, 1, 1, rnd)
	if big {
		for i := 0; i < n; i++ {
			b.Data[i] *= bignum
		}
	}

	// Compute the LU factorization of A with full pivoting.
	lu := cloneGeneral(a)
	ipiv := make([]int, n)
	jpiv := make([]int, n)
	impl.Dgetc2(n, lu.Data, lu.Stride, ipiv, jpiv)

	// Make copies of const input to Dgesc2.
	luCopy := cloneGeneral(lu)
	ipivCopy := make([]int, len(ipiv))
	copy(ipivCopy, ipiv)
	jpivCopy := make([]int, len(jpiv))
	copy(jpivCopy, jpiv)

	// Call Dgesc2 to solve A*x = scale*b.
	x := cloneGeneral(b)
	scale := impl.Dgesc2(n, lu.Data, lu.Stride, x.Data, ipiv, jpiv)

	if n == 0 {
		return
	}

	// Check that const input to Dgesc2 hasn't been modified.
	if !floats.Same(lu.Data, luCopy.Data) {
		t.Errorf("%v: unexpected modification in lu", name)
	}
	if !intsEqual(ipiv, ipivCopy) {
		t.Errorf("%v: unexpected modification in ipiv", name)
	}
	if !intsEqual(jpiv, jpivCopy) {
		t.Errorf("%v: unexpected modification in jpiv", name)
	}

	if scale <= 0 || 1 < scale {
		t.Errorf("%v: scale %v out of bounds (0,1]", name, scale)
	}
	if !big && scale != 1 {
		t.Errorf("%v: unexpected scaling, scale=%v", name, scale)
	}

	// Compute the difference rhs := A*x - scale*b.
	diff := b
	for i := 0; i < n; i++ {
		diff.Data[i] *= scale
	}
	blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, a, x, -1, diff)

	// Compute the residual |A*x - scale*b| / |x|.
	xnorm := dlange(lapack.MaxColumnSum, n, 1, x.Data, 1)
	resid := dlange(lapack.MaxColumnSum, n, 1, diff.Data, 1) / xnorm
	if resid > tol || math.IsNaN(resid) {
		t.Errorf("%v: unexpected result; resid=%v, want<=%v", name, resid, tol)
	}
}