dgehrd.go

// Copyright ©2016 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"
)

type Dgehrder interface {
	Dgehrd(n, ilo, ihi int, a []float64, lda int, tau, work []float64, lwork int)

	Dorgqr(m, n, k int, a []float64, lda int, tau, work []float64, lwork int)
}

func DgehrdTest(t *testing.T, impl Dgehrder) {
	rnd := rand.New(rand.NewSource(1))

	// Randomized tests for small matrix sizes that will most likely
	// use the unblocked algorithm.
	for _, n := range []int{1, 2, 3, 4, 5, 10, 34} {
		for _, extra := range []int{0, 13} {
			for _, optwork := range []bool{true, false} {
				for cas := 0; cas < 10; cas++ {
					ilo := rnd.Intn(n)
					ihi := rnd.Intn(n)
					if ilo > ihi {
						ilo, ihi = ihi, ilo
					}
					testDgehrd(t, impl, n, ilo, ihi, extra, optwork, rnd)
				}
			}
		}
	}

	// These are selected tests for larger matrix sizes to test the blocked
	// algorithm. Use sizes around several powers of two because that is
	// where the blocked path will most likely start to be taken. For
	// example, at present the blocked algorithm is used for sizes larger
	// than 129.
	for _, test := range []struct {
		n, ilo, ihi int
	}{
		{0, 0, -1},

		{68, 0, 63},
		{68, 0, 64},
		{68, 0, 65},
		{68, 0, 66},
		{68, 0, 67},

		{132, 2, 129},
		{132, 1, 129}, // Size = 129, unblocked.
		{132, 0, 129}, // Size = 130, blocked.
		{132, 1, 130},
		{132, 0, 130},
		{132, 1, 131},
		{132, 0, 131},

		{260, 2, 257},
		{260, 1, 257},
		{260, 0, 257},
		{260, 0, 258},
		{260, 0, 259},
	} {
		for _, extra := range []int{0, 13} {
			for _, optwork := range []bool{true, false} {
				testDgehrd(t, impl, test.n, test.ilo, test.ihi, extra, optwork, rnd)
			}
		}
	}
}

func testDgehrd(t *testing.T, impl Dgehrder, n, ilo, ihi, extra int, optwork bool, rnd *rand.Rand) {
	const tol = 1e-13

	a := randomGeneral(n, n, n+extra, rnd)
	aCopy := a
	aCopy.Data = make([]float64, len(a.Data))
	copy(aCopy.Data, a.Data)

	var tau []float64
	if n > 1 {
		tau = nanSlice(n - 1)
	}

	var work []float64
	if optwork {
		work = nanSlice(1)
		impl.Dgehrd(n, ilo, ihi, a.Data, a.Stride, tau, work, -1)
		work = nanSlice(int(work[0]))
	} else {
		work = nanSlice(max(1, n))
	}

	impl.Dgehrd(n, ilo, ihi, a.Data, a.Stride, tau, work, len(work))

	if n == 0 {
		// Just make sure there is no panic.
		return
	}

	prefix := fmt.Sprintf("Case n=%v, ilo=%v, ihi=%v, extra=%v", n, ilo, ihi, extra)

	// Check any invalid modifications of a.
	if !generalOutsideAllNaN(a) {
		t.Errorf("%v: out-of-range write to A\n%v", prefix, a.Data)
	}
	for i := ilo; i <= ihi; i++ {
		for j := 0; j < min(ilo, i); j++ {
			if a.Data[i*a.Stride+j] != aCopy.Data[i*aCopy.Stride+j] {
				t.Errorf("%v: unexpected modification of A[%v,%v]", prefix, i, j)
			}
		}
	}
	for i := ihi + 1; i < n; i++ {
		for j := 0; j < i; j++ {
			if a.Data[i*a.Stride+j] != aCopy.Data[i*aCopy.Stride+j] {
				t.Errorf("%v: unexpected modification of A[%v,%v]", prefix, i, j)
			}
		}
	}
	for i := 0; i <= ilo; i++ {
		for j := i; j < ilo+1; j++ {
			if a.Data[i*a.Stride+j] != aCopy.Data[i*aCopy.Stride+j] {
				t.Errorf("%v: unexpected modification at A[%v,%v]", prefix, i, j)
			}
		}
		for j := ihi + 1; j < n; j++ {
			if a.Data[i*a.Stride+j] != aCopy.Data[i*aCopy.Stride+j] {
				t.Errorf("%v: unexpected modification at A[%v,%v]", prefix, i, j)
			}
		}
	}
	for i := ihi + 1; i < n; i++ {
		for j := i; j < n; j++ {
			if a.Data[i*a.Stride+j] != aCopy.Data[i*aCopy.Stride+j] {
				t.Errorf("%v: unexpected modification at A[%v,%v]", prefix, i, j)
			}
		}
	}

	// Check that tau has been assigned properly.
	for i, v := range tau {
		if math.IsNaN(v) {
			t.Errorf("%v: unexpected NaN at tau[%v]", prefix, i)
		}
	}

	// Extract Q and check that it is orthogonal.
	q := eye(n, n)
	if ilo != ihi {
		for i := ilo + 2; i <= ihi; i++ {
			for j := ilo + 1; j < ihi; j++ {
				q.Data[i*q.Stride+j] = a.Data[i*a.Stride+j-1]
			}
		}
		nh := ihi - ilo
		impl.Dorgqr(nh, nh, nh, q.Data[(ilo+1)*q.Stride+ilo+1:], q.Stride, tau[ilo:ihi], work, len(work))
	}
	if resid := residualOrthogonal(q, false); resid > tol {
		t.Errorf("%v: Q is not orthogonal; resid=%v, want<=%v", prefix, resid, tol)
	}

	// Construct Qᵀ * AOrig * Q and check that it is upper Hessenberg.
	aq := blas64.General{
		Rows:   n,
		Cols:   n,
		Stride: n,
		Data:   make([]float64, n*n),
	}
	blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, aCopy, q, 0, aq)
	qaq := blas64.General{
		Rows:   n,
		Cols:   n,
		Stride: n,
		Data:   make([]float64, n*n),
	}
	blas64.Gemm(blas.Trans, blas.NoTrans, 1, q, aq, 0, qaq)
	for i := 0; i <= ilo; i++ {
		for j := ilo + 1; j <= ihi; j++ {
			qaqij := qaq.Data[i*qaq.Stride+j]
			diff := qaqij - a.Data[i*a.Stride+j]
			if math.Abs(diff) > tol {
				t.Errorf("%v: Qᵀ*AOrig*Q and A are not equal, diff at [%v,%v]=%v", prefix, i, j, diff)
			}
		}
	}
	for i := ilo + 1; i <= ihi; i++ {
		for j := ilo; j < n; j++ {
			qaqij := qaq.Data[i*qaq.Stride+j]
			if j < i-1 {
				if math.Abs(qaqij) > tol {
					t.Errorf("%v: Qᵀ*AOrig*Q is not upper Hessenberg, [%v,%v]=%v", prefix, i, j, qaqij)
				}
				continue
			}
			diff := qaqij - a.Data[i*a.Stride+j]
			if math.Abs(diff) > tol {
				t.Errorf("%v: Qᵀ*AOrig*Q and A are not equal, diff at [%v,%v]=%v", prefix, i, j, diff)
			}
		}
	}
}