// 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 testlapack

import (
	"fmt"
	"math"
	"testing"

	"gonum.org/v1/gonum/blas"
	"gonum.org/v1/gonum/floats"
	"gonum.org/v1/gonum/internal/rand"
	"gonum.org/v1/gonum/lapack"
)

type Dpbconer interface {
	Dpbcon(uplo blas.Uplo, n, kd int, ab []float64, ldab int, anorm float64, work []float64, iwork []int) float64

	Dpbtrser
}

// DpbconTest tests Dpbcon by generating a random symmetric band matrix A and
// checking that the estimated condition number is not too different from the
// condition number computed via the explicit inverse of A.
func DpbconTest(t *testing.T, impl Dpbconer) {
	rnd := rand.New(rand.NewSource(1))
	for _, n := range []int{0, 1, 2, 3, 4, 5, 10, 50} {
		for _, kd := range []int{0, (n + 1) / 4, (3*n - 1) / 4, (5*n + 1) / 4} {
			for _, uplo := range []blas.Uplo{blas.Upper, blas.Lower} {
				for _, ldab := range []int{kd + 1, kd + 1 + 3} {
					dpbconTest(t, impl, uplo, n, kd, ldab, rnd)
				}
			}
		}
	}
}

func dpbconTest(t *testing.T, impl Dpbconer, uplo blas.Uplo, n, kd, ldab int, rnd *rand.Rand) {
	const ratioThresh = 10

	name := fmt.Sprintf("uplo=%v,n=%v,kd=%v,ldab=%v", string(uplo), n, kd, ldab)

	// Generate a random symmetric positive definite band matrix.
	ab := randSymBand(uplo, n, kd, ldab, rnd)

	// Compute the Cholesky decomposition of A.
	abFac := make([]float64, len(ab))
	copy(abFac, ab)
	ok := impl.Dpbtrf(uplo, n, kd, abFac, ldab)
	if !ok {
		t.Fatalf("%v: bad test matrix, Dpbtrf failed", name)
	}

	// Compute the norm of A.
	work := make([]float64, 3*n)
	aNorm := dlansb(lapack.MaxColumnSum, uplo, n, kd, ab, ldab, work)

	// Compute an estimate of rCond.
	iwork := make([]int, n)
	abFacCopy := make([]float64, len(abFac))
	copy(abFacCopy, abFac)
	rCondGot := impl.Dpbcon(uplo, n, kd, abFac, ldab, aNorm, work, iwork)

	if !floats.Equal(abFac, abFacCopy) {
		t.Errorf("%v: unexpected modification of ab", name)
	}

	// Form the inverse of A to compute a good estimate of the condition number
	//  rCondWant := 1/(norm(A) * norm(inv(A)))
	lda := max(1, n)
	aInv := make([]float64, n*lda)
	for i := 0; i < n; i++ {
		aInv[i*lda+i] = 1
	}
	impl.Dpbtrs(uplo, n, kd, n, abFac, ldab, aInv, lda)
	aInvNorm := dlange(lapack.MaxColumnSum, n, n, aInv, lda)
	rCondWant := 1.0
	if aNorm > 0 && aInvNorm > 0 {
		rCondWant = 1 / aNorm / aInvNorm
	}

	ratio := rCondTestRatio(rCondGot, rCondWant)
	if ratio >= ratioThresh {
		t.Errorf("%v: unexpected value of rcond. got=%v, want=%v (ratio=%v)", name, rCondGot, rCondWant, ratio)
	}
}

// rCondTestRatio returns a test ratio to compare two values of the reciprocal
// of the condition number.
//
// This function corresponds to DGET06 in Reference LAPACK.
func rCondTestRatio(rcond, rcondc float64) float64 {
	const eps = dlamchE
	switch {
	case rcond > 0 && rcondc > 0:
		return math.Max(rcond, rcondc)/math.Min(rcond, rcondc) - (1 - eps)
	case rcond > 0:
		return rcond / eps
	case rcondc > 0:
		return rcondc / eps
	default:
		return 0
	}
}