linalg_test.go

package tensor

import (
	"bytes"
	"cmp"
	_ "embed"
	"fmt"
	"io"
	"math"
	"math/cmplx"
	"slices"
	"testing"
)

func TestArnoldi(t *testing.T) {
	t.Parallel()
	o := func(opt ArnoldiOptions) *ArnoldiOptions { return &opt }
	type testcase struct {
		a      *Dense
		opt    *ArnoldiOptions
		lambda []complex64
		tol    float32
	}
	tests := []testcase{
		{
			a:      T2([][]complex64{{-2, 0, 0, 0}, {0, -3, -4, 0}, {0, -4, -9, 0}, {0, 0, 0, 5}}),
			opt:    o(NewArnoldiOptions().KrylovSpaceDim(3)),
			lambda: []complex64{-11, -2},
			tol:    1e-4,
		},
		{
			a: T2([][]complex64{
				{-1, 0, 0, 0, -1},
				{0, 1, -1, 0, 0},
				{0, -1, 1, 0, 0},
				{0, 0, 0, -1, 0},
				{-1, 0, 0, 0, -1},
			}),
			lambda: []complex64{-2, -1, 0},
			tol:    5e-5,
		},
		{
			a: T2([][]complex64{
				{-3, 1i, 1i, 0, 1i, 0, 0, 0, 1i, 0, 0, 0, 0, 0, 0, 0},
				{-1i, -1, 0, 1i, 0, 1i, 0, 0, 0, 1i, 0, 0, 0, 0, 0, 0},
				{-1i, 0, 1, 1i, 0, 0, 1i, 0, 0, 0, 1i, 0, 0, 0, 0, 0},
				{0, -1i, -1i, -1, 0, 0, 0, 1i, 0, 0, 0, 1i, 0, 0, 0, 0},
				{-1i, 0, 0, 0, 1, 1i, 1i, 0, 0, 0, 0, 0, 1i, 0, 0, 0},
				{0, -1i, 0, 0, -1i, 3, 0, 1i, 0, 0, 0, 0, 0, 1i, 0, 0},
				{0, 0, -1i, 0, -1i, 0, 1, 1i, 0, 0, 0, 0, 0, 0, 1i, 0},
				{0, 0, 0, -1i, 0, -1i, -1i, -1, 0, 0, 0, 0, 0, 0, 0, 1i},
				{-1i, 0, 0, 0, 0, 0, 0, 0, -1, 1i, 1i, 0, 1i, 0, 0, 0},
				{0, -1i, 0, 0, 0, 0, 0, 0, -1i, 1, 0, 1i, 0, 1i, 0, 0},
				{0, 0, -1i, 0, 0, 0, 0, 0, -1i, 0, 3, 1i, 0, 0, 1i, 0},
				{0, 0, 0, -1i, 0, 0, 0, 0, 0, -1i, -1i, 1, 0, 0, 0, 1i},
				{0, 0, 0, 0, -1i, 0, 0, 0, -1i, 0, 0, 0, -1, 1i, 1i, 0},
				{0, 0, 0, 0, 0, -1i, 0, 0, 0, -1i, 0, 0, -1i, 1, 0, 1i},
				{0, 0, 0, 0, 0, 0, -1i, 0, 0, 0, -1i, 0, -1i, 0, -1, 1i},
				{0, 0, 0, 0, 0, 0, 0, -1i, 0, 0, 0, -1i, 0, -1i, -1i, -3},
			}),
			opt:    o(NewArnoldiOptions().KrylovSpaceDim(5)),
			lambda: []complex64{-4.758770483143627, -4.064177772475917, -2.7587704831436306},
			tol:    5e-5,
		},
		{
			a: T2([][]complex64{
				{-3i, -1, -1, 0, -1, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0},
				{-1, -1i, 0, -1, 0, -1, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0},
				{-1, 0, 1i, -1, 0, 0, -1, 0, 0, 0, -1, 0, 0, 0, 0, 0},
				{0, -1, -1, -1i, 0, 0, 0, -1, 0, 0, 0, -1, 0, 0, 0, 0},
				{-1, 0, 0, 0, 1i, -1, -1, 0, 0, 0, 0, 0, -1, 0, 0, 0},
				{0, -1, 0, 0, -1, 3i, 0, -1, 0, 0, 0, 0, 0, -1, 0, 0},
				{0, 0, -1, 0, -1, 0, 1, -1, 0, 0, 0, 0, 0, 0, -1, 0},
				{0, 0, 0, -1, 0, -1, -1, -1, 0, 0, 0, 0, 0, 0, 0, -1},
				{-1, 0, 0, 0, 0, 0, 0, 0, -1, -1, -1, 0, -1, 0, 0, 0},
				{0, -1, 0, 0, 0, 0, 0, 0, -1, 1, 0, -1, 0, -1, 0, 0},
				{0, 0, -1, 0, 0, 0, 0, 0, -1, 0, 3i, -1, 0, 0, -1, 0},
				{0, 0, 0, -1, 0, 0, 0, 0, 0, -1, -1, 1i, 0, 0, 0, -1},
				{0, 0, 0, 0, -1, 0, 0, 0, -1, 0, 0, 0, -1i, -1, -1, 0},
				{0, 0, 0, 0, 0, -1, 0, 0, 0, -1, 0, 0, -1, 1i, 0, -1},
				{0, 0, 0, 0, 0, 0, -1, 0, 0, 0, -1, 0, -1, 0, -1i, -1},
				{0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, -1, 0, -1, -1, -3i},
			}),
			opt:    o(NewArnoldiOptions().KrylovSpaceDim(9)),
			lambda: []complex64{-3.4558382 - 0.0648431i, -2.1588192 - 0.28336975i, -1.56791219 + 1.64454048i},
			tol:    5e-5,
		},
		{
			a:      readMatrixMarketMust(bytes.NewBuffer(west0479Mtx)),
			lambda: []complex64{-100.88510419 - 66.60624907i, -100.88510419 + 66.60624907i, -74.65352091},
			tol:    0.05,
		},
	}

	for i, test := range tests {
		t.Run(fmt.Sprintf("%d", i), func(t *testing.T) {
			t.Parallel()
			opt := NewArnoldiOptions()
			if test.opt != nil {
				opt = *test.opt
			}
			opt.debug = true

			var bufs [7]*Dense
			for i := range len(bufs) {
				bufs[i] = Zeros(1)
			}
			eigvals, eigvecs := Zeros(1), Zeros(1)
			var err error
			err = Arnoldi(eigvals, eigvecs, test.a, len(test.lambda), bufs, opt)
			if err != nil {
				t.Fatalf("%+v", err)
			}
			sortConj(eigvals, eigvecs, test.lambda)

			for j := range test.lambda {
				// Check eigenvalue.
				lambda := eigvals.At(j)
				if diff := abs(lambda - test.lambda[j]); diff > test.tol {
					t.Fatalf("%d %f lambda %f %f", j, diff, lambda, test.lambda[j])
				}

				// Check eigenvector is normalized.
				vec := eigvecs.Slice([][2]int{{0, eigvecs.Shape()[0]}, {j, j + 1}})
				if vecNorm := vec.FrobeniusNorm(); absf(vecNorm-1) > 50*epsilon {
					t.Fatalf("%f", vecNorm)
				}

				// Check eigenvector relation.
				av := Zeros(1)
				MatMul(av, test.a, vec)
				lambdaVec := Zeros(vec.Shape()...).Set([]int{0, 0}, vec)
				lambdaVec.Mul(lambda)
				if err := equal2(av, lambdaVec.ToSlice2(), test.tol); err != nil {
					t.Fatalf("lambda %v %+v", lambda, err)
				}
			}
		})
	}
}

func TestEig(t *testing.T) {
	t.Parallel()
	type testcase struct {
		a      *Dense
		lambda []complex64
		tol    float32
	}
	tests := []testcase{
		{
			a:      T2([][]complex64{{2i}}),
			lambda: []complex64{2i},
			tol:    epsilon,
		},
		{
			a:      T2([][]complex64{{1, 0}, {1, 3}}),
			lambda: []complex64{1, 3},
			tol:    epsilon,
		},
		{
			a:      T2([][]complex64{{2, -1, -1, 0}, {-1, 3, -1, -1}, {-1, -1, 3, -1}, {0, -1, -1, 2}}),
			lambda: []complex64{0, 2, 4, 4},
			tol:    5e-7,
		},
		{
			a: T2([][]complex64{
				{2, 0, -1, 0},
				{0, 2 + 1i, -1, 0},
				{0, -1, 2 + 1i, 0},
				{0, -1, 0, 2},
			}),
			lambda: []complex64{1 + 1i, 2, 2, 3 + 1i},
			tol:    5e-7,
		},
		{
			a:      T1([]complex64{99, 99, 99, 99, 2, 0, 0, 0, 0, 2 - 1i, -1, -1, -1, -1, 2 - 1i, 0, 0, 0, 0, 2}).Conj().Reshape(5, 4).Slice([][2]int{{1, 5}, {0, 4}}).Transpose(1, 0),
			lambda: []complex64{1 + 1i, 2, 2, 3 + 1i},
			tol:    5e-7,
		},
		{
			a: T2([][]complex64{
				{1 + 2i, 3 + 4i, 5 + 6i, 7 + 8i},
				{9 + 10i, 11 + 12i, 13 + 14i, 15 + 16i},
				{17 + 18i, 19 + 20i, 21 + 22i, 23 + 24i},
				{25 + 26i, 27 + 28i, 29 + 30i, 31 + 32i},
			}),
			lambda: []complex64{-4.65441 - 4.41284i, 0, 0, 68.6544 + 72.4128i},
			tol:    1e-4,
		},
		{
			a: T2([][]complex64{
				{-3, -1, -1, 0, -1, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0},
				{-1, -1, 0, -1, 0, -1, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0},
				{-1, 0, 1, -1, 0, 0, -1, 0, 0, 0, -1, 0, 0, 0, 0, 0},
				{0, -1, -1, -1, 0, 0, 0, -1, 0, 0, 0, -1, 0, 0, 0, 0},
				{-1, 0, 0, 0, 1, -1, -1, 0, 0, 0, 0, 0, -1, 0, 0, 0},
				{0, -1, 0, 0, -1, 3, 0, -1, 0, 0, 0, 0, 0, -1, 0, 0},
				{0, 0, -1, 0, -1, 0, 1, -1, 0, 0, 0, 0, 0, 0, -1, 0},
				{0, 0, 0, -1, 0, -1, -1, -1, 0, 0, 0, 0, 0, 0, 0, -1},
				{-1, 0, 0, 0, 0, 0, 0, 0, -1, -1, -1, 0, -1, 0, 0, 0},
				{0, -1, 0, 0, 0, 0, 0, 0, -1, 1, 0, -1, 0, -1, 0, 0},
				{0, 0, -1, 0, 0, 0, 0, 0, -1, 0, 3, -1, 0, 0, -1, 0},
				{0, 0, 0, -1, 0, 0, 0, 0, 0, -1, -1, 1, 0, 0, 0, -1},
				{0, 0, 0, 0, -1, 0, 0, 0, -1, 0, 0, 0, -1, -1, -1, 0},
				{0, 0, 0, 0, 0, -1, 0, 0, 0, -1, 0, 0, -1, 1, 0, -1},
				{0, 0, 0, 0, 0, 0, -1, 0, 0, 0, -1, 0, -1, 0, -1, -1},
				{0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, -1, 0, -1, -1, -3},
			}),
			lambda: []complex64{-4.758770483143627, -4.064177772475917, -2.7587704831436306, -2.064177772475913, -1.6945927106677185, -1.0, -0.9999999999999996, -0.30540728933227784, 0.3054072893322784, 0.9999999999999996, 1.0000000000000013, 1.6945927106677217, 2.064177772475911, 2.758770483143641, 4.064177772475913, 4.75877048314363},
			tol:    1e-5,
		},
		{
			a: T2([][]complex64{
				{-3i, -1, -1, 0, -1, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0},
				{-1, -1i, 0, -1, 0, -1, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0},
				{-1, 0, 1i, -1, 0, 0, -1, 0, 0, 0, -1, 0, 0, 0, 0, 0},
				{0, -1, -1, -1i, 0, 0, 0, -1, 0, 0, 0, -1, 0, 0, 0, 0},
				{-1, 0, 0, 0, 1i, -1, -1, 0, 0, 0, 0, 0, -1, 0, 0, 0},
				{0, -1, 0, 0, -1, 3i, 0, -1, 0, 0, 0, 0, 0, -1, 0, 0},
				{0, 0, -1, 0, -1, 0, 1, -1, 0, 0, 0, 0, 0, 0, -1, 0},
				{0, 0, 0, -1, 0, -1, -1, -1, 0, 0, 0, 0, 0, 0, 0, -1},
				{-1, 0, 0, 0, 0, 0, 0, 0, -1, -1, -1, 0, -1, 0, 0, 0},
				{0, -1, 0, 0, 0, 0, 0, 0, -1, 1, 0, -1, 0, -1, 0, 0},
				{0, 0, -1, 0, 0, 0, 0, 0, -1, 0, 3i, -1, 0, 0, -1, 0},
				{0, 0, 0, -1, 0, 0, 0, 0, 0, -1, -1, 1i, 0, 0, 0, -1},
				{0, 0, 0, 0, -1, 0, 0, 0, -1, 0, 0, 0, -1i, -1, -1, 0},
				{0, 0, 0, 0, 0, -1, 0, 0, 0, -1, 0, 0, -1, 1i, 0, -1},
				{0, 0, 0, 0, 0, 0, -1, 0, 0, 0, -1, 0, -1, 0, -1i, -1},
				{0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, -1, 0, -1, -1, -3i},
			}),
			lambda: []complex64{-3.455839856357165 - 0.0648422747440924i, -2.1588188751945228 - 0.28336972471472155i, -1.567912187389367 + 1.6445404806590809i, -1.459600689929873 - 1.3983418213792158i, -1.4070590877538072 + 0.21653026290640678i, -0.47270775847220614 + 0.40556302467474487i, -0.273300762797651 - 0.5514005648391669i, -1.665334536937735e-16 + 1.0000000000000002i, 0.05260411888458892 - 2.460493665711876i, 0.11263021588753293 + 2.258005523574051i, 0.5079513527649886 - 0.5835628320875962i, 1.1425522887255422 - 1.5843845132664092i, 1.4325610570119445 + 1.8159456390073405i, 1.6385600581365172 - 0.5627208970384561i, 2.379717436393563 + 0.15180057382597673i, 3.5286626900899227 - 0.0032692108660648464i},
			tol:    1e-5,
		},
	}

	for i, test := range tests {
		t.Run(fmt.Sprintf("%d", i), func(t *testing.T) {
			t.Parallel()
			var bufs [3]*Dense
			for i := range len(bufs) {
				bufs[i] = Zeros(1)
			}
			a := Zeros(test.a.Shape()...).Set([]int{0, 0}, test.a)
			eigvals, eigvecs := Zeros(1), Zeros(1)
			if err := Eig(eigvals, eigvecs, test.a, bufs); err != nil {
				t.Fatalf("%+v", err)
			}
			sortConj(eigvals, eigvecs, test.lambda)

			m := a.Shape()[0]
			for j := range m {
				lambda := eigvals.At(j)
				if diff := abs(lambda - test.lambda[j]); diff > test.tol {
					t.Fatalf("%d %f %v %v", j, diff, lambda, test.lambda[j])
				}

				vec := eigvecs.Slice([][2]int{{0, m}, {j, j + 1}})
				if vecNorm := vec.FrobeniusNorm(); absf(vecNorm-1) > 2*epsilon {
					t.Fatalf("lambda %v |v| %f", lambda, vecNorm)
				}
				av := gemm(a, vec)
				lambdaVec := Zeros(vec.Shape()...).Set([]int{0, 0}, vec)
				lambdaVec.Mul(lambda)
				if err := equal2(av, lambdaVec.ToSlice2(), test.tol); err != nil {
					t.Fatalf("lambda %v %+v", lambda, err)
				}
			}
		})
	}
}

func TestInverseIteration(t *testing.T) {
	t.Parallel()
	type testcase struct {
		a      *Dense
		lambda []complex64
		tol    float32
		long   bool
	}
	tests := []testcase{
		{
			a:      T2([][]complex64{{2, -1, -1, 0}, {-1, 3, -1, -1}, {-1, -1, 3, -1}, {0, -1, -1, 2}}),
			lambda: []complex64{0, 2, 4, 4},
			tol:    5e-7,
		},
	}

	// This test is taken from Section 11.5, Lecture notes of Numerical Methods for Solving Large Scale Eigenvalue Problems, Peter Arbenz.
	tc := testcase{a: readMatrixMarketMust(bytes.NewBuffer(west0479Mtx)), long: true}
	westLambda0 := []complex64{-100.88510419200173 - 66.60624906782236i, -100.88510419200173 + 66.60624906782236i, -74.65352090884986, -35.662104406279056, -35.16048283061636 - 39.397763510664156i, -35.16048283061636 + 39.397763510664156i, -33.73891457387876, -31.67979017809047 - 17.125483696218136i, -31.67979017809047 + 17.125483696218136i, -25.21593825618421 - 25.216770026116333i, -25.21593825618421 + 25.216770026116333i, -23.300845391687574 - 70.68947896043062i, -23.300845391687574 + 70.68947896043062i, -17.82510732753974 - 4.637637141479285i, -17.82510732753974 + 4.637637141479285i, -16.969487311381044 - 30.551622793202235i, -16.969487311381044 + 30.551622793202235i, -16.93249021573771, -14.98768017256926 - 7.282789765247085i, -14.98768017256926 + 7.282789765247085i}
	westLambda1 := []complex64{12.818688736767768, 13.403040912557852 - 12.213944881604004i, 13.403040912557852 + 12.213944881604004i, 17.548546066690633 - 34.23782252989614i, 17.548546066690633 + 34.23782252989614i, 18.109185928572288 - 4.664987087471094i, 18.109185928572288 + 4.664987087471094i, 25.21603743418249 - 25.216850032975813i, 25.21603743418249 + 25.216850032975813i, 33.70695304316405 - 17.55672234252947i, 33.70695304316405 + 17.55672234252947i, 33.87148153603233, 35.6618691257841, 43.06194325775711 - 39.164280664139575i, 43.06194325775711 + 39.164280664139575i, 59.788970139362874 - 43.68881135483681i, 59.788970139362874 + 43.68881135483681i, 74.63543908467804, 108.12525583925499 - 54.06593856030257i, 108.12525583925499 + 54.06593856030257i}
	tc.lambda = make([]complex64, 0)
	tc.lambda = append(tc.lambda, westLambda0...)
	for _ = range tc.a.Shape()[0] - len(westLambda0) - len(westLambda1) {
		tc.lambda = append(tc.lambda, complex64(cmplx.NaN()))
	}
	tc.lambda = append(tc.lambda, westLambda1...)
	tc.tol = 0.002
	tests = append(tests, tc)

	for i, test := range tests {
		t.Run(fmt.Sprintf("%d", i), func(t *testing.T) {
			t.Parallel()
			if testing.Short() && test.long {
				t.Skip()
			}
			a := Zeros(test.a.Shape()...).Set([]int{0, 0}, test.a)
			var bufs [4]*Dense
			for i := range len(bufs) {
				bufs[i] = Zeros(1)
			}
			ebufs := [3]*Dense(bufs[:3])
			eigvals := Zeros(1)
			if err := Eig(eigvals, nil, test.a, ebufs); err != nil {
				t.Fatalf("%+v", err)
			}
			sortConj(eigvals, nil, test.lambda)

			for j := range eigvals.Shape()[0] {
				if cmplx.IsNaN(complex128(test.lambda[j])) {
					continue
				}
				// Copy tensors to run them in parallel.
				aj := Zeros(a.Shape()...).Set([]int{0, 0}, a)
				valsRaw := Zeros(eigvals.Shape()...).Set([]int{0}, eigvals)
				t.Run(fmt.Sprintf("%d", j), func(t *testing.T) {
					t.Parallel()
					bufs := []*Dense{}
					for k := 0; k < 4; k++ {
						bufs = append(bufs, Zeros(1))
					}
					vec := Zeros(1)
					lambda, err := InverseIteration(vec, aj, valsRaw.At(j), bufs)
					if err != nil {
						t.Fatalf("%v %+v", valsRaw.At(j), err)
					}
					if diff := abs(lambda - test.lambda[j]); diff > test.tol {
						t.Fatalf("%d %f %v %v", j, diff, lambda, test.lambda[j])
					}

					if vecNorm := vec.FrobeniusNorm(); absf(vecNorm-1) > epsilon {
						t.Fatalf("%v %f", lambda, vecNorm)
					}
					av := gemm(aj, vec)
					lambdaVec := Zeros(vec.Shape()...).Set([]int{0, 0}, vec)
					lambdaVec.Mul(lambda)
					if err := equal2(av, lambdaVec.ToSlice2(), test.tol); err != nil {
						t.Fatalf("lambda %v %+v", lambda, err)
					}
				})
			}
		})
	}
}

func TestFindUnreducedHessenberg(t *testing.T) {
	t.Parallel()
	tests := []struct {
		a *Dense
		p int
		q int
	}{
		{
			a: T2([][]complex64{
				{1, 1, 1, 1, 1, 1},
				{1, 1, 1, 1, 1, 1},
				{0, 1, 1, 1, 1, 1},
				{0, 0, 1, 1, 1, 1},
				{0, 0, 0, 1, 1, 1},
				{0, 0, 0, 0, 1, 1},
			}),
			p: 0, q: 0,
		},
		{
			a: T2([][]complex64{
				{1, 1, 1, 1, 1, 1},
				{1, 1, 1, 1, 1, 1},
				{0, 1, 1, 1, 1, 1},
				{0, 0, 1, 1, 1, 1},
				{0, 0, 0, 1, 1, 1},
				{0, 0, 0, 0, 0.5 * epsilon, 1},
			}),
			p: 0, q: 1,
		},
		{
			a: T2([][]complex64{
				{1, 1, 1, 1, 1, 1},
				{0.5 * epsilon, 1, 1, 1, 1, 1},
				{0, 1, 1, 1, 1, 1},
				{0, 0, 1, 1, 1, 1},
				{0, 0, 0, 1, 1, 1},
				{0, 0, 0, 0, 1, 1},
			}),
			p: 1, q: 0,
		},
		{
			a: T2([][]complex64{
				{1, 1, 1, 1, 1, 1},
				{0.5 * epsilon, 1, 1, 1, 1, 1},
				{0, 1, 1, 1, 1, 1},
				{0, 0, 1, 1, 1, 1},
				{0, 0, 0, 0.5 * epsilon, 1, 1},
				{0, 0, 0, 0, 0.5 * epsilon, 1},
			}),
			p: 1, q: 2,
		},
		{
			a: T2([][]complex64{
				{1, 1, 1, 1, 1, 1},
				{0, 1, 1, 1, 1, 1},
				{0, 1, 1, 1, 1, 1},
				{0, 0, 0, 1, 1, 1},
				{0, 0, 0, 1, 1, 1},
				{0, 0, 0, 0, 0, 1},
			}),
			p: 3, q: 1,
		},
	}
	for i, test := range tests {
		t.Run(fmt.Sprintf("%d", i), func(t *testing.T) {
			t.Parallel()
			p, q := findUnreducedHessenberg(test.a)
			if p != test.p {
				t.Fatalf("%d %d", p, test.p)
			}
			if q != test.q {
				t.Fatalf("%d %d", q, test.q)
			}
		})
	}
}

func TestBalance(t *testing.T) {
	t.Parallel()
	tests := []struct {
		a *Dense
		d [][]complex64
		b [][]complex64
	}{
		{
			a: T2([][]complex64{
				{1, 2, 0},
				{9, 1, 0.01},
				{1, 2, 31.4159},
			}),
			d: [][]complex64{
				{0.5, 0, 0},
				{0, 1, 0},
				{0, 0, 1},
			},
			b: [][]complex64{
				{1, 4, 0},
				{4.5, 1, 0.01},
				{0.5, 2, 31.4159},
			},
		},
	}
	for i, test := range tests {
		t.Run(fmt.Sprintf("%d", i), func(t *testing.T) {
			t.Parallel()
			b := Zeros(test.a.Shape()...).Set([]int{0, 0}, test.a)
			d := Zeros(1)
			balance(b, d)

			// Check that a = d @ b @ dInv
			dInv := Zeros(d.Shape()...)
			for j := range d.Shape()[0] {
				dInv.SetAt([]int{j, j}, 1/d.At(j, j))
			}
			dbd := gemm(d, b, dInv)
			if err := equal2(dbd, test.a.ToSlice2(), epsilon); err != nil {
				t.Fatalf("%+v", err)
			}

			if err := equal2(b, test.b, epsilon); err != nil {
				t.Fatalf("%+v", err)
			}
			if err := equal2(d, test.d, epsilon); err != nil {
				t.Fatalf("%+v", err)
			}
		})
	}
}

func TestHessenberg(t *testing.T) {
	t.Parallel()
	tests := []struct {
		a *Dense
		h [][]complex64
		q [][]complex64
	}{
		{
			a: T2([][]complex64{{2, -1, -1, 0}, {-1, 3, -1, -1}, {-1, -1, 3, -1}, {0, -1, -1, 2}}),
			h: [][]complex64{
				{2, 1.41421356, 0, 0},
				{1.41421356, 2, 1.41421356, 0},
				{0, 1.41421356, 2, 0},
				{0, 0, 0, 4},
			},
			q: [][]complex64{
				{1, 0, 0, 0},
				{0, -0.707106781, 0, -0.707106781},
				{0, -0.707106781, 0, 0.707106781},
				{0, 0, 1, 0},
			},
		},
		{
			a: T2([][]complex64{
				{2, 0, -1, 0},
				{0, 2 + 1i, -1, 0},
				{0, -1, 2 + 1i, 0},
				{0, -1, 0, 2},
			}),
			h: [][]complex64{
				{2, 0, 0.70710678, 0.70710678i},
				{0, 2 + 1i, 0.70710678, 0.70710678i},
				{0, 1.4142136, 2 + 0.5i, -0.5},
				{0, 0, 0.5, 2 + 0.5i},
			},
			q: [][]complex64{
				{1, 0, 0, 0},
				{0, 1, 0, 0},
				{0, 0, -0.70710678, -0.70710678i},
				{0, 0, -0.70710678, 0.70710678i},
			},
		},
		{
			a: readMatrixMarketMust(bytes.NewBuffer(west0479Mtx)),
		},
	}
	for i, test := range tests {
		t.Run(fmt.Sprintf("%d", i), func(t *testing.T) {
			t.Parallel()
			bufs := []*Dense{Zeros(1), Zeros(1)}
			h := Zeros(test.a.Shape()...).Set([]int{0, 0}, test.a)
			q := Zeros(1)
			hessenberg(h, q, bufs)

			// Check q @ h @ q.H == a.
			a := gemm(q, h, q.H())
			if err := equal2(a, test.a.ToSlice2(), 3*epsilon*test.a.FrobeniusNorm()); err != nil {
				t.Fatalf("%+v", err)
			}

			// Check q @ q.H == I.
			delta := 5 * epsilon * float32(test.a.Shape()[0])
			qq := gemm(q, q.H())
			if err := equal2(qq, Zeros(1).Eye(qq.Shape()[0], 0).ToSlice2(), delta); err != nil {
				t.Fatalf("%+v", err)
			}

			// Check h is hessenberg.
			hlow := Zeros(h.Shape()...).Set([]int{0, 0}, h)
			hlow.Tril(-2)
			if err := equal2(hlow, Zeros(hlow.Shape()...).ToSlice2(), delta); err != nil {
				t.Fatalf("%+v", err)
			}

			if len(test.h) > 0 {
				// Match ground truth basis.
				phase := gemm(q.H(), T2(test.q))
				q.Set([]int{0, 0}, gemm(q, phase))
				h.Set([]int{0, 0}, gemm(phase.H(), h, phase))

				// Check results.
				if err := equal2(h, test.h, 2*epsilon*test.a.FrobeniusNorm()); err != nil {
					t.Fatalf("%+v", err)
				}
				if err := equal2(q, test.q, 2*epsilon*test.a.FrobeniusNorm()); err != nil {
					t.Fatalf("%+v", err)
				}
			}
		})
	}
}

func TestQRTall(t *testing.T) {
	t.Parallel()
	type testcase struct {
		a    *Dense
		full bool
		q    [][]complex64
		r    [][]complex64
	}
	tests := []testcase{
		{
			a: T2([][]complex64{{1, 2}, {3, 6}}),
			q: [][]complex64{{0.316228, -0.948683}, {0.948683, 0.316228}},
			r: [][]complex64{{3.16228, 6.32456}, {0, 0}},
		},
		{
			a: T2([][]complex64{{12, -51, 4}, {6, 167, -68}, {-4, 24, -41}}),
			q: [][]complex64{{0.857143, -0.394286, -0.331429}, {0.428571, 0.902857, 0.0342857}, {-0.285714, 0.171429, -0.942857}},
			r: [][]complex64{{14, 21, -14}, {0, 175, -70}, {0, 0, 35}},
		},
		{
			a: T2([][]complex64{{99, 1, 3, 5}, {99, 1 + 2i, 4, 6}}).Conj().Transpose(1, 0).Slice([][2]int{{1, 4}, {0, 2}}),
			q: [][]complex64{{0.169031, -0.113478 - 0.964563i}, {0.507093, 0.156032 + 0.0851085i}, {0.845154, -0.0709238 + 0.141848i}},
			r: [][]complex64{{5.91608, 7.26833 - 0.338062i}, {0, 2.01424}, {0, 0}},
		},
		{
			a: T2([][]complex64{{1, 1 - 2i}, {3, 4}, {5, 6}, {1 + 3i, 4 + 1i}}),
			q: [][]complex64{{0.149071, -0.028991 - 0.446461i}, {0.447214, 0.173946 + 0.22613i}, {0.745356, 0.115964 + 0.376883i}, {0.149071 + 0.447214i, 0.527636 - 0.533434i}},
			r: [][]complex64{{6.7082, 7.45356 - 1.93793i}, {0, 3.83261}, {0, 0}, {0, 0}},
		},
		{
			a: T2([][]complex64{
				{2, 0, -1, 0},
				{0, 2 + 1i, -1, 0},
				{0, -1, 2 + 1i, 0},
				{0, -1, 0, 2},
			}),
			q: [][]complex64{
				{1, 0, 0, 0},
				{0, 0.75592895 + 0.377964473i, 0.07412493 + 0.296499727i, 0.35082321 + 0.263117406i},
				{0, -0.37796447, 0.74124932 + 0.518874522i, 0.0877058 + 0.175411604i},
				{0, -0.37796447, -0.29649973, 0.87705802},
			},
			r: [][]complex64{
				{2, 0, -1, 0},
				{0, 2.64575131, -1.51185789, -0.75592895},
				{0, 0, 1.92724822, -0.59299945},
				{0, 0, 0, 1.75411604},
			},
		},
		{
			a: T2([][]complex64{
				{-3, -1, -1, 0, -1, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0},
				{-1, -1, 0, -1, 0, -1, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0},
				{-1, 0, 1, -1, 0, 0, -1, 0, 0, 0, -1, 0, 0, 0, 0, 0},
				{0, -1, -1, -1, 0, 0, 0, -1, 0, 0, 0, -1, 0, 0, 0, 0},
				{-1, 0, 0, 0, 1, -1, -1, 0, 0, 0, 0, 0, -1, 0, 0, 0},
				{0, -1, 0, 0, -1, 3, 0, -1, 0, 0, 0, 0, 0, -1, 0, 0},
				{0, 0, -1, 0, -1, 0, 1, -1, 0, 0, 0, 0, 0, 0, -1, 0},
				{0, 0, 0, -1, 0, -1, -1, -1, 0, 0, 0, 0, 0, 0, 0, -1},
				{-1, 0, 0, 0, 0, 0, 0, 0, -1, -1, -1, 0, -1, 0, 0, 0},
				{0, -1, 0, 0, 0, 0, 0, 0, -1, 1, 0, -1, 0, -1, 0, 0},
				{0, 0, -1, 0, 0, 0, 0, 0, -1, 0, 3, -1, 0, 0, -1, 0},
				{0, 0, 0, -1, 0, 0, 0, 0, 0, -1, -1, 1, 0, 0, 0, -1},
				{0, 0, 0, 0, -1, 0, 0, 0, -1, 0, 0, 0, -1, -1, -1, 0},
				{0, 0, 0, 0, 0, -1, 0, 0, 0, -1, 0, 0, -1, 1, 0, -1},
				{0, 0, 0, 0, 0, 0, -1, 0, 0, 0, -1, 0, -1, 0, -1, -1},
				{0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, -1, 0, -1, -1, -3},
			}),
			q: [][]complex64{
				{-0.8320502943378436, -0.03962144258751633, -0.24943939134627785, 0.1945297739776206, -0.15932444118273437, -0.015451474519767363, -0.03732659503993789, 0.09952306197527283, 0.1640846421930999, 0.1411390233795794, -0.02918026975551966, 0.023249218246914925, 0.24480383298369912, 0.00642540148334711, 0.10548842571838889, -0.2231713410606682},
				{-0.2773500981126145, -0.35659298328764744, 0.19955151307702232, -0.25532032834562723, 0.11246431142310666, -0.3418638737498538, 0.2775151196447555, 0.0814279597979505, 0.15855369919782702, -0.2846295413933638, 0.3012360233753443, -0.007366236312794194, 0.0018499388453828564, -0.3109610839358521, -0.17766471699939185, 0.3905498468561696},
				{-0.2773500981126145, 0.15848577035006547, 0.508856358346407, -0.37690143708163987, -0.18744051903851117, 0.16996621971744136, -0.05177036442495726, -1.0408340855860843e-16, -0.07743320193382264, 0.31212817396398385, 0.12954211035199034, -0.3701505468166283, -0.2716208349088028, 0.24602889067783687, -0.18876876181185387, -0.05579283526516712},
				{0.0, -0.5150787536377128, -0.3093048452693846, -0.3769014370816399, 0.22492862284621323, 0.11588605889825551, -0.014281479841367434, -0.04976153098763646, 0.0663713159432763, 0.07566916147357444, -0.2613900185865692, -0.35807406883667603, -0.19486970682425736, 0.11058818206598926, 0.4164016804673248, -8.326672684688674e-17},
				{-0.2773500981126145, 0.15848577035006547, 0.01995515130770223, 0.024316221747202566, 0.5342054792597566, 0.07532593828386608, -0.26988751109311615, -0.23523632830519028, 0.011061885990546146, -0.04648329828971118, -0.00764434470434592, 0.36030319922657184, -0.561881134992518, -0.16207548805740937, -0.022208089624924175, -0.05579283526516711},
				{0.0, -0.5150787536377128, 0.1795963617693201, 0.21884599572482327, -0.34676496022124553, 0.5543216483966554, -0.1635878600011191, -0.17642724622889275, 0.06637131594327643, -0.2689586041902713, 0.18421394118140055, 0.20967965171398595, -0.09538480496519804, 0.06801739645587433, -0.027760112031155047, 5.551115123125783e-17},
				{0.0, 0.0, -0.48890120703870465, -0.0972648869888103, -0.3842530640289478, -0.13133753341802298, 0.1871198438306451, -0.348330716913455, 0.10324426924509661, 0.12507142566111232, -0.15110927158137633, 0.09763760781819678, -0.3236050095151754, 0.11629200170562094, -0.4746979157327502, 0.16737850579550115},
				{0.0, 0.0, 0.0, -0.49848254581765283, -0.07497620761540447, -0.17189765403241236, -0.4467831136736914, -0.2940454103814879, -0.03687295330182015, 0.09308806200039813, 0.12509237478300425, 0.40964055092751084, 0.39700355756453415, 0.2011365307759001, 0.13880056015577494, 0.13948208816291763},
				{-0.2773500981126145, 0.15848577035006547, 0.01995515130770223, 0.024316221747202566, 0.0187440519038511, 0.14292613930784848, 0.15612254099313153, -0.14476081741857869, -0.5844363098338504, -0.404432404419647, -0.33559297975642965, -0.05253407083789416, 0.0972405321048404, 0.2077314768653833, 0.07217629128100311, 0.39054984685616945},
				{0.0, -0.5150787536377128, 0.1795963617693201, 0.21884599572482327, 0.16869646713466002, -0.3128923590252899, -0.062319184762331115, 0.04523775544330587, -0.4553809732774797, 0.3367799607304812, -0.19487967621465593, 0.13251143518856925, 0.04360073996037339, 0.1259301039306414, -0.3164652771551668, -0.1673785057955012},
				{0.0, 0.0, -0.48890120703870465, -0.0972648869888103, 0.13120836332695773, 0.2008691687569762, -0.18728213337429708, 0.29856918592581855, -0.41113342931529545, -0.02975143655028218, 0.5712216702754556, -0.13296330404506407, 0.002050048446930741, 0.012723305422879735, -0.2359609522648175, 7.45931094670027e-17},
				{0.0, 0.0, 0.0, -0.49848254581765283, -0.07497620761540447, 0.22790924916656918, 0.23531983829526054, 0.262378981571174, -0.11061885990546058, -0.19625585604459272, -0.2944804899237695, 0.3259503010385879, 0.06763704532314334, -0.24679251958387427, -0.18876876181185415, -0.47423909975392},
				{0.0, 0.0, 0.0, 0.0, -0.5154614273559055, -0.33220670217499915, -0.2560928998827043, 0.18999857286188454, -0.3226383413909269, -0.04373514314013155, -0.011568744548850478, 0.029736721447474135, -0.38769515349584366, -0.3528102877022517, 0.3747615124205923, -1.0408340855860843e-16},
				{0.0, 0.0, 0.0, 0.0, 0.0, -0.39980690319898154, -0.051608074881305435, -0.08142795979795052, 0.06637131594327637, -0.568851034888137, 0.13395777009019907, -0.1335385586993304, -0.12312677631299325, 0.4759522305849843, -0.022208089624923995, -0.4742390997539201},
				{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -0.6304948770876465, 0.18095102177322334, 0.20648853849019322, -0.2336615120135585, -0.39809941201202503, -0.30215559551925764, 0.11289340282161123, -0.16879793169070675, -0.4025216244517471, 0.13948208816291782},
				{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -0.6559474539279349, -0.19911394782982916, -0.02487004494481361, 0.10319297420354073, -0.3588837416230172, 0.2168559637400966, -0.49603411100338457, -0.05274421285919445, -0.3068605939584188},
			},
			r: [][]complex64{
				{3.6055512754639896, 1.1094003924504583, 0.5547001962252291, 0.554700196225229, 0.5547001962252291, 0.554700196225229, 0.554700196225229, 0.0, 1.1094003924504583, 0.554700196225229, 0.554700196225229, 0.0, 0.554700196225229, 0.0, 0.0, 0.0},
				{0.0, 1.9414506867883017, 0.7131859665752946, 0.7131859665752949, 0.7131859665752946, -1.3471290479755564, -0.31697154070013095, 1.0301575072754257, 0.3962144258751636, -0.3169715407001309, -0.31697154070013095, 1.0301575072754257, -0.31697154070013095, 1.0301575072754257, 0.0, 0.0},
				{0.0, 0.0, 2.0454030090394792, -0.39910302615404447, 0.5786993879233647, 0.3192824209232358, -1.017712716692814, 0.6186096905387692, 0.5387890853079601, -0.03991030261540444, -1.9955151307702228, 0.6186096905387691, -0.03991030261540446, -0.3591927235386402, 0.9778024140774093, 0.0},
				{0.0, 0.0, 0.0, 2.0060882941442135, -0.29179466096643103, 1.3860246395905473, 0.7538028741632797, 0.7538028741632798, -0.3404271044608361, 0.9483326481409007, 0.5592731001856591, -0.24316221747202582, -0.048632443494405125, -0.43769199144964654, 0.19452977397762058, 0.9969650916353057},
				{0.0, 0.0, 0.0, 0.0, 1.94000937204859, -1.6119884637311954, -0.6560418166347889, 0.5810656090193845, 0.3561369861731711, 0.11246431142310669, 0.6372977647309377, -0.5998096609232355, -0.03748810380770218, 0.693529920442491, 0.7685061280578955, 0.14995241523080893},
				{0.0, 0.0, 0.0, 0.0, 0.0, 2.501207437887349, -0.20473203738691803, -0.36697251984447565, 0.31675522765523173, 0.0579430294491277, 0.061805898079069665, 0.22404638053662745, 0.513761527782266, -0.309029490395348, 0.2626750668360459, 0.3437953080648247},
				{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.5860557101100565, 0.43753260968553276, 0.38689827206613886, -0.6796686088141732, -0.27102353789867933, 0.49920263627325606, 1.0519608219516408, 0.430391869764849, 0.8867500665140028, 0.8935662273473828},
				{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.5245123584394067, -0.48856775878770314, -0.07238040870928938, 0.597138371851637, 0.6242810251176208, 0.09047551088661168, 0.5157104120536866, 0.33475939028046336, 1.8999857286188457},
				{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.609504411624453, 0.014749181320728, -0.667400454762946, 0.8886381745738674, 0.6231529108007616, 0.9771332624982358, 0.6231529108007617, 0.47197380226329877},
				{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.790948797476222, 0.4329672888629678, -0.5540834967535526, 1.2971633927511852, -0.5680672033434018, 0.20694671098767345, 0.9802904757803308},
				{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 2.6122957821666004, -0.512625439601541, 0.6189477109314518, 0.05299927546876419, -0.1136372163367445, 0.12395083445196897},
				{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0433599803547755, 0.09818830438243636, -0.1465826254263427, 0.6666283119216679, 0.776754527121541},
				{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.8625691298749029, 0.09949647844757825, 0.3795007480023805, -1.104975120616585},
				{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.130849128904105, 0.8886270232678423, 1.2266040229238504},
				{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.7911631928879169, 0.6329305543103333},
				{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.590095805057261},
			},
		},
	}

	var tc testcase
	tc.a = readMatrixMarketMust(bytes.NewBuffer(west0479Mtx))
	tc.a = tc.a.Slice([][2]int{{0, 479}, {80, 112}})
	tests = append(tests, tc)

	tc = testcase{}
	tc.a = readMatrixMarketMust(bytes.NewBuffer(qc324Mtx))
	tc.a = tc.a.Slice([][2]int{{0, 324}, {0, 32}})
	tests = append(tests, tc)

	// Add tests for full QR.
	testsLen := len(tests)
	for i := range testsLen {
		ti := tests[i]

		tc := ti
		tc.a = Zeros(ti.a.Shape()...).Set([]int{0, 0}, ti.a)
		tc.full = true
		tests = append(tests, tc)
	}

	for i, test := range tests {
		t.Run(fmt.Sprintf("%d", i), func(t *testing.T) {
			t.Parallel()
			a := Zeros(test.a.Shape()...).Set([]int{0, 0}, test.a)
			q, bufs := Zeros(1), [2]*Dense{Zeros(1), Zeros(1)}
			r := QR(q, test.a, bufs, QROptions{Full: test.full})

			// Check shape of r.
			m, n := a.Shape()[0], a.Shape()[1]
			k := n
			if test.full {
				k = m
			}
			if !slices.Equal(r.Shape(), []int{k, n}) {
				t.Fatalf("%#v", r.Shape())
			}

			// Check that r is upper triangular.
			rlow := Zeros(r.Shape()...).Set([]int{0, 0}, r)
			rlow.Tril(-1)
			if err := equal2(rlow, Zeros(rlow.Shape()...).ToSlice2(), epsilon); err != nil {
				t.Fatalf("%+v", err)
			}

			// Check that a = q @ r.
			qr := Zeros(1)
			MatMul(qr, q, r)
			if err := equal2(qr, a.ToSlice2(), epsilon*a.FrobeniusNorm()); err != nil {
				t.Fatalf("%+v %#v", err, qr.ToSlice2())
			}

			// Check q is unitary.
			qqt := MatMul(Zeros(1), q.H(), q)
			identity := Zeros(1).Eye(qqt.Shape()[0], 0).ToSlice2()
			if err := equal2(qqt, identity, 5*epsilon); err != nil {
				t.Fatalf("%+v %#v", err, qqt.ToSlice2())
			}

			if test.r != nil {
				rFull := Zeros(a.Shape()...).Set([]int{0, 0}, r)
				qThin := q.Slice([][2]int{{0, q.Shape()[0]}, {0, test.a.Shape()[1]}})
				if err := equal2(rFull, test.r, 5*epsilon*a.FrobeniusNorm()); err != nil {
					t.Fatalf("%+v", err)
				}
				if err := equal2(qThin, test.q, epsilon*a.FrobeniusNorm()); err != nil {
					t.Fatalf("%+v", err)
				}
			}
		})
	}
}

func TestQRFat(t *testing.T) {
	t.Parallel()
	type testcase struct {
		a *Dense
		q [][]complex64
		r [][]complex64
	}
	tests := []testcase{
		{
			a: T2([][]complex64{
				{1 + 1i, 5 + 3i, -1, 4 - 1i, 3 - 2i},
				{-2 + 7i, -4, 2 + 1i, 5, -1 + 3i},
			}),
			q: [][]complex64{{0.13483997 + 0.13483997i, 0.9012942 + 0.3889796i}, {-0.26967994 + 0.94387981i, -0.01897461 - 0.18974615i}},
			r: [][]complex64{
				{7.41619849, 2.15743956 + 3.50583928i, 0.26967994 - 2.02259959i, -0.94387981 - 5.3935989i, 3.23615934 - 0.53935989i},
				{0, 5.74930826, -1.12898958 + 0.74949728i, 3.12132412 - 1.50848187i, 1.37565957 - 3.21619719i},
			},
		},
	}

	var tc testcase
	tc.a = readMatrixMarketMust(bytes.NewBuffer(qc324Mtx))
	tc.a = tc.a.Slice([][2]int{{0, 32}, {0, 324}})
	tests = append(tests, tc)

	for i, test := range tests {
		t.Run(fmt.Sprintf("%d", i), func(t *testing.T) {
			t.Parallel()
			a := Zeros(test.a.Shape()...).Set([]int{0, 0}, test.a)
			q, bufs := Zeros(1), [2]*Dense{Zeros(1), Zeros(1)}
			r := QR(q, test.a, bufs)

			// Check shape of r.
			m, n := a.Shape()[0], a.Shape()[1]
			if !slices.Equal(r.Shape(), []int{m, n}) {
				t.Fatalf("%#v", r.Shape())
			}

			// Check that r is upper triangular.
			rlow := Zeros(r.Shape()...).Set([]int{0, 0}, r)
			rlow.Tril(-1)
			if err := equal2(rlow, Zeros(rlow.Shape()...).ToSlice2(), epsilon); err != nil {
				t.Fatalf("%+v", err)
			}

			// Check that a = q @ r.
			qr := Zeros(1)
			MatMul(qr, q, r)
			if err := equal2(qr, a.ToSlice2(), epsilon*a.FrobeniusNorm()); err != nil {
				t.Fatalf("%+v %#v", err, qr.ToSlice2())
			}

			// Check q is unitary.
			qqt := gemm(q.H(), q)
			identity := Zeros(1).Eye(qqt.Shape()[0], 0).ToSlice2()
			if err := equal2(qqt, identity, 5*epsilon); err != nil {
				t.Fatalf("%+v %#v", err, qqt.ToSlice2())
			}

			if test.r != nil {
				if err := equal2(r, test.r, 5*epsilon*a.FrobeniusNorm()); err != nil {
					t.Fatalf("%+v", err)
				}
				if err := equal2(q, test.q, epsilon*a.FrobeniusNorm()); err != nil {
					t.Fatalf("%+v", err)
				}
			}
		})
	}
}

func TestAreIndependent(t *testing.T) {
	t.Parallel()
	tests := []struct {
		x           *Dense
		independent bool
	}{
		{
			x: T2([][]complex64{
				{1, 7, -2},
				{4, 10, 1},
				{2, -4, 5},
				{-3, -1, -4},
			}),
			independent: false,
		},
		{
			x: T2([][]complex64{
				{0, 1, 1},
				{0, 2, 2},
				{1, 1, 3},
			}),
			independent: false,
		},
		{
			x: T2([][]complex64{
				{1, 2, 3, 4},
				{2, 4, 5, 6},
				{1, 3, 5, 7},
				{1, 2, 6, 7},
				{2, 1, 3, 0},
				{1, 1, 1, 0},
				{9, 1, 4, 0},
			}),
			independent: true,
		},
		{
			x: T2([][]complex64{
				{1, -2, 1},
				{2, 1, 0},
				{3, 0, 1},
			}),
			independent: true,
		},
	}
	for i, _ := range tests {
		t.Run(fmt.Sprintf("%d", i), func(t *testing.T) {
			t.Parallel()
		})
	}
}

//go:embed testdata/impcol_b.mtx
var impcol_bMtx []byte

//go:embed testdata/qc324.mtx
var qc324Mtx []byte

func TestSVD(t *testing.T) {
	t.Parallel()
	type testcase struct {
		a    *Dense
		full bool
		u    [][]complex64
		s    [][]complex64
		v    [][]complex64
	}
	tests := []testcase{
		{
			a: T2([][]complex64{{1, -2}, {0, 1 + 1i}, {-1, 4}, {-1i, 3}}),
			u: [][]complex64{
				{-0.34472358 - 0.15293662i, 0.4764441 - 0.16312692i, 0.58859905 + 0.10224976i, -0.11890409 - 0.48282982i},
				{0.07646831 + 0.22940492i, 0.08156346 + 0.24469038i, -0.52249756 - 0.04600586i, -0.35369817 - 0.6890409i},
				{0.65059682 + 0.30587323i, -0.15019027 + 0.32625384i, 0.57292684 - 0.0797522i, -0.07013681 - 0.12100377i},
				{0.45880985 + 0.26825528i, 0.48938076 - 0.55800756i, -0.18200195 + 0.0156722i, 0.36182605 - 0.04876728i},
			},
			s: [][]complex64{
				{5.70299699, 0},
				{0, 1.2148355},
				{0, 0},
				{0, 0},
			},
			v: [][]complex64{
				{-0.22156345, 0.97514596},
				{0.87219706 + 0.43609853i, 0.19817237 + 0.09908619i},
			},
		},
		{
			a: T2([][]complex64{
				{1, 0, -1, -1i, 3},
				{-2, 1 + 1i, 4, 3, 2},
			}),
			u: [][]complex64{{-0.12724529, 0.99187128}, {-0.99187128i, -0.12724529i}},
			s: [][]complex64{
				{5.94851783, 0, 0, 0, 0},
				{0, 3.40809854, 0, 0, 0},
			},
			v: [][]complex64{
				{-0.0213910916 + 0.33348518i, 0.291033627 + 0.07467231i, 0.461529494 - 0.03769366i, 0.128596775 + 0.23007181i, -0.713406041 - 0.07214495i},
				{-0.166742592 - 0.16674259i, -0.0373361534 - 0.03733615i, -0.591093152 + 0.1416856i, -0.437591477 + 0.21038468i, -0.519152769 + 0.2585968i},
				{0.0213910916 - 0.66697037i, -0.291033627 - 0.14934461i, 0.587645463 + 0.05353303i, -0.262427854 + 0.03581186i, -0.117772423 + 0.12439605i},
				{-0.52161887i, -0 + 0.17902517i, -0.240397136 + 0.08644503i, 0.751316369 + 0.05757735i, -0.218832283 - 0.12215242i},
				{-0.0641732749 - 0.33348518i, 0.873100881 - 0.07467231i, 0.0132236457 - 0.04972348i, -0.149533992 + 0.18568547i, 0.239262014 - 0.00743043i},
			},
		},
		{
			a: T2([][]complex64{
				{1, 0, -1, -1i, 3, -1 + 3i, 2, 4},
				{-2, 1 + 1i, 4, 3, 2, 3 + 1i, 3 + 4i, 1},
				{-2, 1 - 1i, 4, 3i, 2, 1 + 1i, 3 + 2i, 1},
				{-4, 1 + 2i, 2, 3, 2, 2 + 1i, 2 + 4i, -1},
				{-1, 2 + 3i, 4, 3, 2i, 1 - 1i, 1 + 5i, 1i},
			}),
			u: [][]complex64{
				{0.09871383314072041, 0.7973938334803816, -0.44129644303876037, 0.2616993695569026, -0.3019760633046683},
				{0.3060425521966734 + 0.471885727765128i, 0.04885230596061252 - 0.1960745333670441i, -0.11936282028039988 + 0.03447396194948467i, 0.22458683797826312 - 0.05847216402235553i, 0.5981062899812457 - 0.46454788015047627i},
				{0.20307176657797243 + 0.3193410907550637i, 0.25672723215148424 - 0.2306672581482087i, 0.7307977185349942 - 0.34965711310211856i, 0.07819230852857553 + 0.07139527160630535i, -0.2559038977887196 + 0.06814206249745541i},
				{0.2247701791773497 + 0.4538842964030999i, -0.11757931321383443 - 0.13441499018849118i, -0.30411991068813965 + 0.09741441377779265i, -0.5451310365467851 + 0.4867178610730561i, -0.2649961927769231 + 0.07288027153513955i},
				{0.15955893676943228 + 0.49870720711225797i, -0.2805465106493334 + 0.30596800333398777i, -0.11129228421587009 + 0.138317206702829i, 0.2859435511001721 - 0.5006176570840692i, -0.2782049881787783 + 0.3349804664840379i},
			},
			s: [][]complex64{
				{14.551468657151489, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
				{0.0, 7.413728201472017, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
				{0.0, 0.0, 4.60230182303232, 0.0, 0.0, 0.0, 0.0, 0.0},
				{0.0, 0.0, 0.0, 3.2532115388961285, 0.0, 0.0, 0.0, 0.0},
				{0.0, 0.0, 0.0, 0.0, 1.2356482863431857, 0.0, 0.0, 0.0},
			},
			v: [][]complex64{
				{-0.13594191105343975 - 0.26778726749687654i, 0.12640044189572414 + 0.14637379614685872i, -0.07309262310786042 + 0.022247876047339465i, 0.476674096685511 - 0.4525066952380117i, 0.2847160599703515 + 0.13459338269600069i, 0.0937404898020845 - 0.20473595847281478i, 0.09987049446103947 + 0.14673497368654073i, -0.5008004667469982 + 0.04899318112589205i},
				{0.24804624932961664 + 0.08324503893088375i, 0.041891737809336166 + 0.18013235355264667i, 0.2343701007292268 + 0.4022210971922822i, -0.10105382593285267 - 0.12773760642345816i, 0.1123363536911454 + 0.6935893964671251i, 0.09310080546227309 - 0.008687169142269878i, -0.14939075357502138 - 0.004765626736818213i, 0.2220136005270687 + 0.282916119746413i},
				{0.20791850078331212 + 0.41696854374586123i, -0.12576942729999857 - 0.10142334769606866i, 0.3984152166294159 - 0.11138664302244178i, 0.30828900481457233 - 0.3004238931803355i, 0.02263855671147122 - 0.08088131769390598i, -0.2876430450662767 - 0.05901662432447771i, -0.21933000612502163 - 0.4614189839409071i, -0.0761456636227946 - 0.20069161304950367i},
				{0.2081671855992936 + 0.25859453199584265i, -0.234675763781961 - 0.006252242108029146i, -0.5765150759085603 - 0.3961219671056371i, 0.033930739238469856 - 0.05840180813861186i, 0.2987431930662472 + 0.23928670232095034i, -0.015132147700947295 - 0.22257733075558123i, 0.26374639632798275 - 0.12421272777931254i, 0.2390884222510238 + 0.04908834939821921i},
				{0.18976262634449498 + 0.1492017340284553i, 0.4559268784947508 - 0.07570017767813378i, -0.09400186194061252 - 0.04627095625328497i, -0.21543049912307 + 0.1313762815616627i, -0.06600167621814153 - 0.07335510992893612i, 0.4401905528429548 - 0.5302023694182854i, -0.27599769334917473 - 0.20196495763804762i, -0.2257392351553058 - 0.030651788193294935i},
				{0.16341916056656933 + 0.15606720303438723i, -0.23968232409828283 - 0.4913157013980958i, -0.056843805338454775 + 0.21558622239272845i, 0.21092804708512033 + 0.034425728462986485i, 0.3024223423892285 - 0.23812889908246157i, 0.4587241290952476 + 0.2987300115657367i, -0.14068387156065928 + 0.09213853585692233i, -0.046715635549771185 + 0.27643834314791893i},
				{0.6301197801754223 + 0.03112861826199379i, 0.23501921071580187 - 0.010642585145824526i, 0.16339758717033825 + 0.03832723177807859i, -0.006116201521515902 + 0.06383286393986931i, -0.11413142588332176 - 0.11175032351589076i, -0.08066304244605901 + 0.11084708055272419i, 0.6391909038790421 + 0.07338225799916456i, -0.21345183439578086 + 0.10263340217768016i},
				{0.08094761477519581 + 0.01221756989183949i, 0.5285737596643169 - 0.0015889833505214686i, -0.15455608603279622 - 0.06546838784163775i, 0.42852731447958764 - 0.23353486101523852i, -0.21504922128899845 - 0.15464036418880414i, 0.009215596510732844 + 0.10551443946919796i, -0.12741119112680024 + 0.15722823653881218i, 0.5655718054693626 + 0.07316935952434622i},
			},
		},
	}

	s := []complex64{8.639495362414722, 2.897126890781191, 2.8366406167144533, 2.5803873901384407, 2.411006449071305, 2.1899969198186167, 2.136666617057314, 2.113752806954914, 2.105143953798741, 2.0485769077462512, 2.0037722304437473, 1.9835341684022283, 1.9628310979325356, 1.873306846558814, 1.7342803928775503, 1.7230800460647926, 1.692299983821104, 1.652691518442473, 1.643358106894265, 1.63786549718296, 1.6194606089626091, 1.6180339887498996, 1.5656148058917694, 1.5119211165142545, 1.416500577528672, 1.3809530383265074, 1.2817034663543523, 1.266097844245871, 1.2467094479567649, 1.135713250611134, 1.1133839856659247, 1.063757209670854, 1.0278452792156727, 0.9999993918236418, 0.9999325772238793, 0.9998733988809205, 0.9980639022697806, 0.9940142699787611, 0.9516357441535507, 0.9310695392521922, 0.8563664163463203, 0.8346037661066318, 0.721326167746477, 0.6698881938855177, 0.6323502476443208, 0.6180339887499015, 0.5792072095031812, 0.45712964981921844, 0.446196356288565, 0.3989925037636857, 0.3784254215511466, 0.3354549026494168, 0.24450269933082402, 0.2354863863618566, 0.15452065584215918, 0.13468956398368512, 0.0762446393021703, 0.0091023264771938, 5.288146843365497e-05}
	var tc testcase
	tc.a = readMatrixMarketMust(bytes.NewBuffer(impcol_bMtx))
	tc.s = Zeros(tc.a.Shape()...).ToSlice2()
	for i := range len(tc.s[0]) {
		tc.s[i][i] = s[i]
	}
	tests = append(tests, tc)

	s = []complex64{1.5231092179708894, 1.2913163831685748, 1.1814603543690512, 1.1120441349872188, 0.95904024598977, 0.9500090863951773, 0.8283022027003175, 0.771366457825332, 0.7164311954393427, 0.6205892691423285, 0.6192157758420567, 0.5383982798496774, 0.5124244379473237, 0.5108884235406859, 0.5091627768855584, 0.5072618921366342, 0.5051988654116006, 0.5029839842302147, 0.5006254004858618, 0.49812972119326265, 0.49550241798795613, 0.4927481071933595, 0.48987074566671357, 0.48965118298390786, 0.4868737689752382, 0.4837601983751689, 0.48053271398399267, 0.47719371642851416, 0.473745372779706, 0.47018965342470104, 0.4678648818028705, 0.46652836177015994, 0.462763158491201, 0.45889558157180593, 0.4568917638944124, 0.454927063062763, 0.45085894326180354, 0.4466924828706182, 0.4424288735956233, 0.4412401543206885, 0.43806924764501404, 0.433614686829, 0.43098604921381845, 0.4290662312781046, 0.4244249027245899, 0.4212883362216545, 0.41969161775402536, 0.4148678285875984, 0.4121482266902494, 0.40995241108514635, 0.40731882815559667, 0.4049538214666534, 0.4035579518839678, 0.39985399853979997, 0.39553448485530074, 0.3946719247997497, 0.3894109287183109, 0.3880657247748941, 0.3840552219448776, 0.38115968870531863, 0.3794455889591577, 0.3786190967355648, 0.3748148775793116, 0.3730959821060025, 0.36903406866367433, 0.3674901823649941, 0.36381506935149077, 0.3618040662518525, 0.3591621664176325, 0.35603643173495736, 0.35536223623813384, 0.355077320915448, 0.35156418968643693, 0.3501869180125571, 0.34862215000502417, 0.34625746121894585, 0.34447392358365303, 0.3442543487880477, 0.34324254564127143, 0.34226677391666854, 0.34108733909506267, 0.33966522022941337, 0.33825101033407184, 0.33803719088836404, 0.3362359244980474, 0.3342695159848866, 0.3321711984158693, 0.3321438982028766, 0.32989067997352334, 0.32749438099950734, 0.32600164816467203, 0.32496655037717687, 0.32231435525775065, 0.3197657788183362, 0.31953720638825933, 0.3168794297285342, 0.3166467516466896, 0.3136428500931249, 0.3134465352377676, 0.3108281409597114, 0.3105179727251623, 0.30729063138714013, 0.30706104223825803, 0.3055729790406376, 0.30394853374657277, 0.300619065280198, 0.3004907624263908, 0.29695524506001864, 0.2940818611201967, 0.2932989666904353, 0.2921684853760476, 0.28955178560979367, 0.28748556532162056, 0.2860860451823678, 0.28569396085080917, 0.2817511485028488, 0.2808227297348214, 0.27967756332101223, 0.2776963658191623, 0.2740955043979908, 0.27356391143157, 0.2727474884003511, 0.2693178980120895, 0.26894598921471974, 0.2673027552572793, 0.26500370283138963, 0.26057786909257813, 0.2604395968109928, 0.2581981205021026, 0.2560787709047218, 0.253538837035241, 0.2514551899925498, 0.2492395297794737, 0.2468055184533124, 0.24656525710428062, 0.24298916041429183, 0.24199200431669532, 0.2403261385816447, 0.2395527981970572, 0.23891488527770519, 0.23717581725609035, 0.23248783207639437, 0.2321809447062517, 0.22988865054236546, 0.22721564833600774, 0.22537023151764063, 0.22202874011014606, 0.22152927726645988, 0.21822734315750833, 0.21692738086088878, 0.21340204976405078, 0.2129064887471767, 0.2115578243951029, 0.21151354018118723, 0.20787184882912432, 0.20636198514609277, 0.20590799872128598, 0.20295137454643392, 0.20083373191535947, 0.19932325990353802, 0.19558338720223759, 0.19456517771962703, 0.19004162428609953, 0.1899237106235257, 0.18980669105967155, 0.18853320845727684, 0.1843522287657691, 0.1823093965933313, 0.178267193157864, 0.1755680769030667, 0.17291846017265924, 0.17244193303875652, 0.1714044733841509, 0.1698355049006135, 0.16957457871381046, 0.1683021538388738, 0.16657963033174347, 0.16644474073130403, 0.16468271415712035, 0.1633068300410197, 0.16262474154094742, 0.16171421738440392, 0.16041619249365438, 0.15806543878967744, 0.15771643749306527, 0.15702056171514583, 0.15661739068803004, 0.15557930017080987, 0.15296350066452025, 0.15255968295337766, 0.15174081033802894, 0.15041361565062664, 0.15022290187221513, 0.1473618667377935, 0.14726625009687205, 0.14622594342378742, 0.14441388139551659, 0.14438405300424595, 0.14346594889168957, 0.1412929663780092, 0.14043580972532108, 0.1400729182240068, 0.13809172314180854, 0.1362277828211618, 0.13500037452574545, 0.13478326813566213, 0.13262519277046197, 0.1313703778554736, 0.12944200833028657, 0.12942636533374918, 0.1293478915185468, 0.1283517228661296, 0.12785568277009682, 0.12520171846050698, 0.12424197420122889, 0.1241403038628657, 0.12165955979861832, 0.12053175717764182, 0.1187921400752021, 0.11789366255483963, 0.11781108083801124, 0.11672775422597954, 0.11366790050968946, 0.11283290954907285, 0.11273049653187431, 0.11069451280755209, 0.10884986594532084, 0.10852823729030586, 0.10481886460758497, 0.10478220692840218, 0.10423162304348055, 0.10361751161357169, 0.10356802233150646, 0.10063325775833241, 0.09871736420526164, 0.09863589862187803, 0.09849555046570718, 0.09666060867987807, 0.09640703974024337, 0.09619610264353122, 0.0938570458327792, 0.09210870781268, 0.08982929032097695, 0.08912076831965791, 0.0885851772512328, 0.08774304774962127, 0.08748901207848524, 0.08452856442551382, 0.08429266794865399, 0.0833154604113938, 0.08314562335791323, 0.07999321973357928, 0.07922633505741128, 0.0790233710517706, 0.07883675482031612, 0.07662584043478372, 0.07558843618224531, 0.0743161542676456, 0.07161820127819879, 0.0712621421892002, 0.0705924443852334, 0.0702813021552061, 0.07003918286504132, 0.06975620797748555, 0.06705397558996416, 0.06519945205016699, 0.064143115319433, 0.06347903756845852, 0.06293157857054237, 0.0625158295391386, 0.061621651831410794, 0.06062989156798679, 0.05892923346652567, 0.058244218687502604, 0.05611470504811367, 0.055019600746583436, 0.053775968644135336, 0.052636621548710974, 0.05166674952986425, 0.05122771209801271, 0.04754253124541987, 0.04739072540288774, 0.047349292393095245, 0.0465281790651325, 0.043969903289519656, 0.043233668345076484, 0.042593155066873244, 0.04053553410333503, 0.0398769760882596, 0.039439806197284516, 0.03839789183829021, 0.03717814334755745, 0.03614217701955978, 0.03498896789897727, 0.034060672219575025, 0.033802884071942034, 0.0336536806021551, 0.03272840857749781, 0.031257525382860085, 0.030943012604683646, 0.029739479055024665, 0.028809450105783455, 0.028199829722004043, 0.027975805753646193, 0.02556650465885333, 0.02397967852358546, 0.02224262204628253, 0.020308937561522092, 0.018622924616054894, 0.01687955528082615, 0.0151797273222276, 0.013583661713646352, 0.012009447286682316, 0.010516931743708977, 0.009097600220962055, 0.007740586372749171, 0.0064708344143663355, 0.005274718464085628, 0.004164361400977617, 0.003141938299624792, 0.002211101216307765, 0.0013840617422473515, 0.0006357728524245715, 3.287699561753757e-05}
	tc = tc
	tc.a = readMatrixMarketMust(bytes.NewBuffer(qc324Mtx))
	tc.s = Zeros(tc.a.Shape()...).ToSlice2()
	for i := range len(tc.s[0]) {
		tc.s[i][i] = s[i]
	}
	tests = append(tests, tc)

	testsLen := len(tests)
	for i := range testsLen {
		ti := tests[i]

		tst := ti
		tst.a = Zeros(ti.a.Shape()...).Set([]int{0, 0}, ti.a)
		tst.full = true
		tests = append(tests, tst)
	}

	for i, test := range tests {
		t.Run(fmt.Sprintf("%d", i), func(t *testing.T) {
			t.Parallel()
			var bufs [3]*Dense
			for i := range len(bufs) {
				bufs[i] = Zeros(1)
			}
			a := Zeros(test.a.Shape()...).Set([]int{0, 0}, test.a)
			u, v := Zeros(1), Zeros(1)
			s, err := SVD(u, v, test.a, bufs, SVDOptions{Full: test.full})
			if err != nil {
				t.Fatalf("%+v", err)
			}

			// Check shape of s.
			minD := min(a.shape[0], a.shape[1])
			aShape := a.Shape()
			if !test.full {
				aShape = []int{minD, minD}
			}
			if !slices.Equal(s.Shape(), aShape) {
				t.Fatalf("%#v %#v", s.Shape(), aShape)
			}

			// Check a = u @ s @ v.H.
			m := float32(a.Shape()[0])
			usv := gemm(u, s, v.H())
			if err := equal2(usv, a.ToSlice2(), m*epsilon*a.FrobeniusNorm()); err != nil {
				t.Fatalf("%+v", err)
			}

			// Check that singular values are correct.
			// See Corollary 8.6.2, Matrix Computations 4th Ed., G. H. Golub, C. F. Van Loan.
			sFull := Zeros(a.Shape()...).Set([]int{0, 0}, s)
			s0 := real(test.s[0][0])
			if err := equal2(sFull, test.s, m*epsilon*s0); err != nil {
				t.Fatalf("%+v", err)
			}

			// Check unitary.
			if err := equal2(gemm(u.H(), u), Zeros(1).Eye(u.Shape()[1], 0).ToSlice2(), 2*m*epsilon); err != nil {
				t.Fatalf("%+v", err)
			}
			if err := equal2(gemm(v.H(), v), Zeros(1).Eye(v.Shape()[1], 0).ToSlice2(), 2*m*epsilon); err != nil {
				t.Fatalf("%+v", err)
			}

			if test.u != nil {
				testu, testv := test.u, test.v
				if !test.full {
					testu = T2(test.u).Slice([][2]int{{0, len(test.u)}, {0, minD}}).ToSlice2()
					testv = T2(test.v).Slice([][2]int{{0, len(test.v)}, {0, minD}}).ToSlice2()
				}

				// Change to ground truth basis.
				trun := [][2]int{{0, minD}, {0, minD}}
				if u.Shape()[0] < v.Shape()[0] {
					phase := gemm(v.H(), T2(testv))
					u.Set([]int{0, 0}, gemm(u, phase.Slice(trun)))
					v.Set([]int{0, 0}, gemm(v, phase))
				} else {
					phase := gemm(u.H(), T2(testu))
					u.Set([]int{0, 0}, gemm(u, phase))
					v.Set([]int{0, 0}, gemm(v, phase.Slice(trun)))
				}
				// Check u and v in ground truth basis.
				for j, vecsPair := range [][2]*Dense{{u, T2(testu)}, {v, T2(testv)}} {
					vecs, truth := vecsPair[0], vecsPair[1]
					for k := range minD {
						ax := [][2]int{{0, vecs.Shape()[0]}, {k, k + 1}}
						vec, vecTruth := vecs.Slice(ax), truth.Slice(ax).ToSlice2()

						// Compute singular vector bounds.
						// See Theorem 8.6.5, Matrix Computations 4th Ed., G. H. Golub, C. F. Van Loan.
						var gap float32 = float32(math.MaxFloat32)
						if k > 0 {
							gap = min(abs(test.s[k][k]-test.s[k-1][k-1]), gap)
						}
						if k+1 < minD {
							gap = min(abs(test.s[k][k]-test.s[k+1][k+1]), gap)
						}
						tol := 5 * epsilon * a.FrobeniusNorm() / gap

						if err := equal2(vec, vecTruth, tol); err != nil {
							t.Fatalf("%d %d %+v", j, k, err)
						}
					}
				}
			}
		})
	}
}

func TestFindBidiagonal(t *testing.T) {
	t.Parallel()
	tests := []struct {
		a *Dense
		p int
		q int
	}{
		{
			a: T2([][]complex64{
				{1, 1, 0, 0, 0, 0},
				{0, 1, 1, 0, 0, 0},
				{0, 0, 1, 1, 0, 0},
				{0, 0, 0, 1, 1, 0},
				{0, 0, 0, 0, 1, 2.1 * epsilon},
				{0, 0, 0, 0, 0, 1},
			}),
			p: 0, q: 0,
		},
		{
			a: T2([][]complex64{
				{1, 1, 0, 0, 0, 0},
				{0, 1, 1, 0, 0, 0},
				{0, 0, 1, 1, 0, 0},
				{0, 0, 0, 1, 1, 0},
				{0, 0, 0, 0, 1, 1.9 * epsilon},
				{0, 0, 0, 0, 0, 1},
			}),
			p: 0, q: 1,
		},
		{
			a: T2([][]complex64{
				{1, 1.9 * epsilon, 0, 0, 0, 0},
				{0, 1, 1, 0, 0, 0},
				{0, 0, 1, 1, 0, 0},
				{0, 0, 0, 1, 1, 0},
				{0, 0, 0, 0, 1, 1},
				{0, 0, 0, 0, 0, 1},
			}),
			p: 1, q: 0,
		},
		{
			a: T2([][]complex64{
				{1, 1, 0, 0, 0, 0},
				{0, 1, 0, 0, 0, 0},
				{0, 0, 1, 0, 0, 0},
				{0, 0, 0, 1, 0, 0},
				{0, 0, 0, 0, 1, 0},
				{0, 0, 0, 0, 0, 1},
			}),
			p: 0, q: 4,
		},
		{
			a: T2([][]complex64{
				{1, 1, 0, 0, 0, 0, 0, 0, 0},
				{0, 1, 0, 0, 0, 0, 0, 0, 0},
				{0, 0, 1, 1, 0, 0, 0, 0, 0},
				{0, 0, 0, 1, 1.9 * epsilon, 0, 0, 0, 0},
				{0, 0, 0, 0, 1, 1, 0, 0, 0},
				{0, 0, 0, 0, 0, 1, 1, 0, 0},
				{0, 0, 0, 0, 0, 0, 1, 0, 0},
				{0, 0, 0, 0, 0, 0, 0, 1, 0},
				{0, 0, 0, 0, 0, 0, 0, 0, 1},
			}),
			p: 4, q: 2,
		},
	}
	for i, test := range tests {
		t.Run(fmt.Sprintf("%d", i), func(t *testing.T) {
			t.Parallel()
			p, q := findBidiagonal(test.a, epsilon, -1)
			if p != test.p {
				t.Fatalf("%d %d", p, test.p)
			}
			if q != test.q {
				t.Fatalf("%d %d", q, test.q)
			}
		})
	}
}

func TestSortEigen(t *testing.T) {
	t.Parallel()
	tests := []struct {
		val         *Dense
		right       *Dense
		left        *Dense
		fn          func(complex64, complex64) int
		sortedVal   []complex64
		sortedRight [][]complex64
		sortedLeft  [][]complex64
	}{
		{
			val: T1([]complex64{1, 3, 2}),
			right: T2([][]complex64{
				{1, 4, 7},
				{2, 5, 8},
				{3, 6, 9},
			}),
			left: T2([][]complex64{
				{-1i, -4i, -7i},
				{-2i, -5i, -8i},
				{-3i, -6i, -9i},
			}),
			fn:        func(a, b complex64) int { return -cmp.Compare(real(a), real(b)) },
			sortedVal: []complex64{3, 2, 1},
			sortedRight: [][]complex64{
				{4, 7, 1},
				{5, 8, 2},
				{6, 9, 3},
			},
			sortedLeft: [][]complex64{
				{-4i, -7i, -1i},
				{-5i, -8i, -2i},
				{-6i, -9i, -3i},
			},
		},
	}
	for i, test := range tests {
		t.Run(fmt.Sprintf("%d", i), func(t *testing.T) {
			t.Parallel()
			buf := Zeros(1)
			sortEigen(test.val, test.right, test.left, test.fn, buf)
			if err := equal1(test.val, test.sortedVal, 0); err != nil {
				t.Fatalf("%+v", err)
			}
			if err := equal2(test.right, test.sortedRight, 0); err != nil {
				t.Fatalf("%+v", err)
			}
			if err := equal2(test.left, test.sortedLeft, 0); err != nil {
				t.Fatalf("%+v", err)
			}
		})
	}
}

func TestBidiagonalize(t *testing.T) {
	t.Parallel()
	tests := []struct {
		a *Dense
	}{
		{
			a: T2([][]complex64{
				{4, 2, 4, 5},
				{3, 1, 4, 6},
				{0, 2, 0, 1},
				{2, 1, 3, 3},
				{5, 6, 0, 7},
			}),
		},
		{
			a: T2([][]complex64{{1, -2}, {0, 1 + 1i}, {-1, 4}, {-1i, 3}}),
		},
		{
			a: T2([][]complex64{
				{1, -2, 1 + 2i},
				{0, 1 - 1i, 3i},
				{-1, 4, 2},
				{1i, 3, 0},
				{3, 2, 2 - 3i},
			}),
		},
	}
	for i, test := range tests {
		t.Run(fmt.Sprintf("%d", i), func(t *testing.T) {
			t.Parallel()
			bufs := []*Dense{Zeros(1), Zeros(1)}
			b := Zeros(test.a.Shape()...).Set([]int{0, 0}, test.a)
			u, v := Zeros(1), Zeros(1)
			bidiagonalize(u, v, b, bufs)

			// Check u @ b @ v.H == a.
			ubvh := gemm(u, b, v.H())
			delta := epsilon * test.a.FrobeniusNorm()
			if err := equal2(ubvh, test.a.ToSlice2(), delta); err != nil {
				t.Fatalf("%+v", err)
			}

			// Check u and v are unitary.
			uu := gemm(u, u.H())
			if err := equal2(uu, Zeros(1).Eye(test.a.Shape()[0], 0).ToSlice2(), 2*epsilon); err != nil {
				t.Fatalf("%+v", err)
			}
			vv := gemm(v, v.H())
			if err := equal2(vv, Zeros(1).Eye(test.a.Shape()[1], 0).ToSlice2(), 2*epsilon); err != nil {
				t.Fatalf("%+v", err)
			}

			// Check b is bidiagonal.
			low := Zeros(b.Shape()...).Set([]int{0, 0}, b)
			low.Tril(-1)
			delta = epsilon * test.a.FrobeniusNorm()
			if err := equal2(low, Zeros(low.Shape()...).ToSlice2(), delta); err != nil {
				t.Fatalf("%+v", err)
			}
			up := Zeros(b.Shape()...).Set([]int{0, 0}, b)
			up.Triu(2)
			delta = epsilon * test.a.FrobeniusNorm()
			if err := equal2(up, Zeros(up.Shape()...).ToSlice2(), delta); err != nil {
				t.Fatalf("%+v", err)
			}

			// Check that b is all real.
			for j := range b.Shape()[1] {
				if imag(b.At(j, j)) != 0 {
					t.Fatalf("%d %d %v", i, j, b.At(i, j))
				}
				if j+1 < b.Shape()[1] {
					if imag(b.At(j, j+1)) != 0 {
						t.Fatalf("%d %d %v", i, j, b.At(i, j+1))
					}
				}
			}
		})
	}
}

func TestConstant(t *testing.T) {
	t.Parallel()

	// Constants for float32
	const (
		minExponent = -126
		maxExponent = 127
	)

	epsilonE := math.Nextafter32(1, 2) - 1
	if epsilon != epsilonE {
		t.Fatalf("%f %f", epsilon, epsilonE)
	}

	safminE := float32(math.Pow(radix, max(minExponent-1, 1-maxExponent)))
	if safmin != safminE {
		t.Fatalf("%f %f", safmin, safminE)
	}
}

// sortConj sorts eigenpairs according to whether the next one is a conjugate.
func sortConj(eigvals, eigvecs *Dense, expected []complex64) {
	var m int = -1
	var veci *Dense
	if eigvecs != nil {
		m = eigvecs.Shape()[0]
		veci = Zeros(m, 1)
	}

	for i := range len(expected) - 1 {
		if expected[i] == conj(expected[i+1]) {
			vali := eigvals.At(i)
			if imag(vali) > imag(eigvals.At(i+1)) {
				eigvals.SetAt([]int{i}, eigvals.At(i+1))
				eigvals.SetAt([]int{i + 1}, vali)

				if eigvecs != nil {
					veci.Set([]int{0, 0}, eigvecs.Slice([][2]int{{0, m}, {i, i + 1}}))
					eigvecs.Set([]int{0, i}, eigvecs.Slice([][2]int{{0, m}, {i + 1, i + 2}}))
					eigvecs.Set([]int{0, i + 1}, veci)
				}
			}
		}
	}
}

func readMatrixMarketMust(r io.Reader) *Dense {
	t, err := ReadMatrixMarket(r)
	if err != nil {
		panic(fmt.Sprintf("%+v", err))
	}
	return t
}