| // 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" |
| "gonum.org/v1/gonum/lapack" |
| ) |
| |
| type Dtrevc3er interface { |
| Dtrevc3(side lapack.EVSide, howmny lapack.EVHowMany, selected []bool, n int, t []float64, ldt int, vl []float64, ldvl int, vr []float64, ldvr int, mm int, work []float64, lwork int) int |
| } |
| |
| func Dtrevc3Test(t *testing.T, impl Dtrevc3er) { |
| rnd := rand.New(rand.NewSource(1)) |
| |
| for _, side := range []lapack.EVSide{lapack.EVRight, lapack.EVLeft, lapack.EVBoth} { |
| var name string |
| switch side { |
| case lapack.EVRight: |
| name = "Rigth" |
| case lapack.EVLeft: |
| name = "Left" |
| case lapack.EVBoth: |
| name = "Both" |
| } |
| t.Run(name, func(t *testing.T) { |
| runDtrevc3Test(t, impl, rnd, side) |
| }) |
| } |
| } |
| |
| func runDtrevc3Test(t *testing.T, impl Dtrevc3er, rnd *rand.Rand, side lapack.EVSide) { |
| for _, n := range []int{0, 1, 2, 3, 4, 5, 6, 7, 10, 34} { |
| for _, extra := range []int{0, 11} { |
| for _, optwork := range []bool{true, false} { |
| for cas := 0; cas < 10; cas++ { |
| dtrevc3Test(t, impl, side, n, extra, optwork, rnd) |
| } |
| } |
| } |
| } |
| } |
| |
| // dtrevc3Test tests Dtrevc3 by generating a random matrix T in Schur canonical |
| // form and performing the following checks: |
| // 1. Compute all eigenvectors of T and check that they are indeed correctly |
| // normalized eigenvectors |
| // 2. Compute selected eigenvectors and check that they are exactly equal to |
| // eigenvectors from check 1. |
| // 3. Compute all eigenvectors multiplied into a matrix Q and check that the |
| // result is equal to eigenvectors from step 1 multiplied by Q and scaled |
| // appropriately. |
| func dtrevc3Test(t *testing.T, impl Dtrevc3er, side lapack.EVSide, n, extra int, optwork bool, rnd *rand.Rand) { |
| const tol = 1e-15 |
| |
| name := fmt.Sprintf("n=%d,extra=%d,optwk=%v", n, extra, optwork) |
| |
| right := side != lapack.EVLeft |
| left := side != lapack.EVRight |
| |
| // Generate a random matrix in Schur canonical form possibly with tiny or zero eigenvalues. |
| // Zero elements of wi signify a real eigenvalue. |
| tmat, wr, wi := randomSchurCanonical(n, n+extra, true, rnd) |
| tmatCopy := cloneGeneral(tmat) |
| |
| // 1. Compute all eigenvectors of T and check that they are indeed correctly |
| // normalized eigenvectors |
| |
| howmny := lapack.EVAll |
| |
| var vr, vl blas64.General |
| if right { |
| // Fill VR and VL with NaN because they should be completely overwritten in Dtrevc3. |
| vr = nanGeneral(n, n, n+extra) |
| } |
| if left { |
| vl = nanGeneral(n, n, n+extra) |
| } |
| |
| var work []float64 |
| if optwork { |
| work = []float64{0} |
| impl.Dtrevc3(side, howmny, nil, n, tmat.Data, tmat.Stride, |
| vl.Data, max(1, vl.Stride), vr.Data, max(1, vr.Stride), n, work, -1) |
| work = make([]float64, int(work[0])) |
| } else { |
| work = make([]float64, max(1, 3*n)) |
| } |
| |
| mGot := impl.Dtrevc3(side, howmny, nil, n, tmat.Data, tmat.Stride, |
| vl.Data, max(1, vl.Stride), vr.Data, max(1, vr.Stride), n, work, len(work)) |
| |
| if !generalOutsideAllNaN(tmat) { |
| t.Errorf("%v: out-of-range write to T", name) |
| } |
| if !equalGeneral(tmat, tmatCopy) { |
| t.Errorf("%v: unexpected modification of T", name) |
| } |
| if !generalOutsideAllNaN(vr) { |
| t.Errorf("%v: out-of-range write to VR", name) |
| } |
| if !generalOutsideAllNaN(vl) { |
| t.Errorf("%v: out-of-range write to VL", name) |
| } |
| |
| mWant := n |
| if mGot != mWant { |
| t.Errorf("%v: unexpected value of m=%d, want %d", name, mGot, mWant) |
| } |
| |
| if right { |
| resid := residualRightEV(tmat, vr, wr, wi) |
| if resid > tol { |
| t.Errorf("%v: unexpected right eigenvectors; residual=%v, want<=%v", name, resid, tol) |
| } |
| resid = residualEVNormalization(vr, wi) |
| if resid > tol { |
| t.Errorf("%v: unexpected normalization of right eigenvectors; residual=%v, want<=%v", name, resid, tol) |
| } |
| } |
| if left { |
| resid := residualLeftEV(tmat, vl, wr, wi) |
| if resid > tol { |
| t.Errorf("%v: unexpected left eigenvectors; residual=%v, want<=%v", name, resid, tol) |
| } |
| resid = residualEVNormalization(vl, wi) |
| if resid > tol { |
| t.Errorf("%v: unexpected normalization of left eigenvectors; residual=%v, want<=%v", name, resid, tol) |
| } |
| } |
| |
| // 2. Compute selected eigenvectors and check that they are exactly equal to |
| // eigenvectors from check 1. |
| |
| howmny = lapack.EVSelected |
| |
| // Follow DCHKHS and select last max(1,n/4) real, max(1,n/4) complex |
| // eigenvectors instead of selecting them randomly. |
| selected := make([]bool, n) |
| selectedWant := make([]bool, n) |
| var nselr, nselc int |
| for j := n - 1; j > 0; { |
| if wi[j] == 0 { |
| if nselr < max(1, n/4) { |
| nselr++ |
| selected[j] = true |
| selectedWant[j] = true |
| } |
| j-- |
| } else { |
| if nselc < max(1, n/4) { |
| nselc++ |
| // Select all columns to check that Dtrevc3 normalizes 'selected' correctly. |
| selected[j] = true |
| selected[j-1] = true |
| selectedWant[j] = false |
| selectedWant[j-1] = true |
| } |
| j -= 2 |
| } |
| } |
| mWant = nselr + 2*nselc |
| |
| var vrSel, vlSel blas64.General |
| if right { |
| vrSel = nanGeneral(n, mWant, n+extra) |
| } |
| if left { |
| vlSel = nanGeneral(n, mWant, n+extra) |
| } |
| |
| if optwork { |
| // Reallocate optimal work in case it depends on howmny and selected. |
| work = []float64{0} |
| impl.Dtrevc3(side, howmny, selected, n, tmat.Data, tmat.Stride, |
| vlSel.Data, max(1, vlSel.Stride), vrSel.Data, max(1, vrSel.Stride), mWant, work, -1) |
| work = make([]float64, int(work[0])) |
| } |
| |
| mGot = impl.Dtrevc3(side, howmny, selected, n, tmat.Data, tmat.Stride, |
| vlSel.Data, max(1, vlSel.Stride), vrSel.Data, max(1, vrSel.Stride), mWant, work, len(work)) |
| |
| if !generalOutsideAllNaN(tmat) { |
| t.Errorf("%v: out-of-range write to T", name) |
| } |
| if !equalGeneral(tmat, tmatCopy) { |
| t.Errorf("%v: unexpected modification of T", name) |
| } |
| if !generalOutsideAllNaN(vrSel) { |
| t.Errorf("%v: out-of-range write to selected VR", name) |
| } |
| if !generalOutsideAllNaN(vlSel) { |
| t.Errorf("%v: out-of-range write to selected VL", name) |
| } |
| |
| if mGot != mWant { |
| t.Errorf("%v: unexpected value of selected m=%d, want %d", name, mGot, mWant) |
| } |
| |
| for i := range selected { |
| if selected[i] != selectedWant[i] { |
| t.Errorf("%v: unexpected selected[%v]", name, i) |
| } |
| } |
| |
| // Check that selected columns of vrSel are equal to the corresponding |
| // columns of vr. |
| var k int |
| match := true |
| if right { |
| loopVR: |
| for j := 0; j < n; j++ { |
| if selected[j] && wi[j] == 0 { |
| for i := 0; i < n; i++ { |
| if vrSel.Data[i*vrSel.Stride+k] != vr.Data[i*vr.Stride+j] { |
| match = false |
| break loopVR |
| } |
| } |
| k++ |
| } else if selected[j] && wi[j] != 0 { |
| for i := 0; i < n; i++ { |
| if vrSel.Data[i*vrSel.Stride+k] != vr.Data[i*vr.Stride+j] || |
| vrSel.Data[i*vrSel.Stride+k+1] != vr.Data[i*vr.Stride+j+1] { |
| match = false |
| break loopVR |
| } |
| } |
| k += 2 |
| } |
| } |
| } |
| if !match { |
| t.Errorf("%v: unexpected selected VR", name) |
| } |
| |
| // Check that selected columns of vlSel are equal to the corresponding |
| // columns of vl. |
| match = true |
| k = 0 |
| if left { |
| loopVL: |
| for j := 0; j < n; j++ { |
| if selected[j] && wi[j] == 0 { |
| for i := 0; i < n; i++ { |
| if vlSel.Data[i*vlSel.Stride+k] != vl.Data[i*vl.Stride+j] { |
| match = false |
| break loopVL |
| } |
| } |
| k++ |
| } else if selected[j] && wi[j] != 0 { |
| for i := 0; i < n; i++ { |
| if vlSel.Data[i*vlSel.Stride+k] != vl.Data[i*vl.Stride+j] || |
| vlSel.Data[i*vlSel.Stride+k+1] != vl.Data[i*vl.Stride+j+1] { |
| match = false |
| break loopVL |
| } |
| } |
| k += 2 |
| } |
| } |
| } |
| if !match { |
| t.Errorf("%v: unexpected selected VL", name) |
| } |
| |
| // 3. Compute all eigenvectors multiplied into a matrix Q and check that the |
| // result is equal to eigenvectors from step 1 multiplied by Q and scaled |
| // appropriately. |
| |
| howmny = lapack.EVAllMulQ |
| |
| var vrMul, qr blas64.General |
| var vlMul, ql blas64.General |
| if right { |
| vrMul = randomGeneral(n, n, n+extra, rnd) |
| qr = cloneGeneral(vrMul) |
| } |
| if left { |
| vlMul = randomGeneral(n, n, n+extra, rnd) |
| ql = cloneGeneral(vlMul) |
| } |
| |
| if optwork { |
| // Reallocate optimal work in case it depends on howmny and selected. |
| work = []float64{0} |
| impl.Dtrevc3(side, howmny, nil, n, tmat.Data, tmat.Stride, |
| vlMul.Data, max(1, vlMul.Stride), vrMul.Data, max(1, vrMul.Stride), n, work, -1) |
| work = make([]float64, int(work[0])) |
| } |
| |
| mGot = impl.Dtrevc3(side, howmny, selected, n, tmat.Data, tmat.Stride, |
| vlMul.Data, max(1, vlMul.Stride), vrMul.Data, max(1, vrMul.Stride), n, work, len(work)) |
| |
| if !generalOutsideAllNaN(tmat) { |
| t.Errorf("%v: out-of-range write to T", name) |
| } |
| if !equalGeneral(tmat, tmatCopy) { |
| t.Errorf("%v: unexpected modification of T", name) |
| } |
| if !generalOutsideAllNaN(vrMul) { |
| t.Errorf("%v: out-of-range write to VRMul", name) |
| } |
| if !generalOutsideAllNaN(vlMul) { |
| t.Errorf("%v: out-of-range write to VLMul", name) |
| } |
| |
| mWant = n |
| if mGot != mWant { |
| t.Errorf("%v: unexpected value of m=%d, want %d", name, mGot, mWant) |
| } |
| |
| if right { |
| // Compute Q * VR explicitly and normalize to match Dtrevc3 output. |
| qvWant := zeros(n, n, n) |
| blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, qr, vr, 0, qvWant) |
| normalizeEV(qvWant, wi) |
| |
| // Compute the difference between Dtrevc3 output and Q * VR. |
| r := zeros(n, n, n) |
| for i := 0; i < n; i++ { |
| for j := 0; j < n; j++ { |
| r.Data[i*r.Stride+j] = vrMul.Data[i*vrMul.Stride+j] - qvWant.Data[i*qvWant.Stride+j] |
| } |
| } |
| qvNorm := dlange(lapack.MaxColumnSum, n, n, qvWant.Data, qvWant.Stride) |
| resid := dlange(lapack.MaxColumnSum, n, n, r.Data, r.Stride) / qvNorm / float64(n) |
| if resid > tol { |
| t.Errorf("%v: unexpected VRMul; resid=%v, want <=%v", name, resid, tol) |
| } |
| } |
| if left { |
| // Compute Q * VL explicitly and normalize to match Dtrevc3 output. |
| qvWant := zeros(n, n, n) |
| blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, ql, vl, 0, qvWant) |
| normalizeEV(qvWant, wi) |
| |
| // Compute the difference between Dtrevc3 output and Q * VL. |
| r := zeros(n, n, n) |
| for i := 0; i < n; i++ { |
| for j := 0; j < n; j++ { |
| r.Data[i*r.Stride+j] = vlMul.Data[i*vlMul.Stride+j] - qvWant.Data[i*qvWant.Stride+j] |
| } |
| } |
| qvNorm := dlange(lapack.MaxColumnSum, n, n, qvWant.Data, qvWant.Stride) |
| resid := dlange(lapack.MaxColumnSum, n, n, r.Data, r.Stride) / qvNorm / float64(n) |
| if resid > tol { |
| t.Errorf("%v: unexpected VLMul; resid=%v, want <=%v", name, resid, tol) |
| } |
| } |
| } |
| |
| // residualEVNormalization returns the maximum normalization error in E: |
| // max |max-norm(E[:,j]) - 1| |
| func residualEVNormalization(emat blas64.General, wi []float64) float64 { |
| n := emat.Rows |
| if n == 0 { |
| return 0 |
| } |
| var ( |
| e = emat.Data |
| lde = emat.Stride |
| enrmin = math.Inf(1) |
| enrmax float64 |
| ipair int |
| ) |
| for j := 0; j < n; j++ { |
| if ipair == 0 && j < n-1 && wi[j] != 0 { |
| ipair = 1 |
| } |
| var nrm float64 |
| switch ipair { |
| case 0: |
| // Real eigenvector |
| for i := 0; i < n; i++ { |
| nrm = math.Max(nrm, math.Abs(e[i*lde+j])) |
| } |
| enrmin = math.Min(enrmin, nrm) |
| enrmax = math.Max(enrmax, nrm) |
| case 1: |
| // Complex eigenvector |
| for i := 0; i < n; i++ { |
| nrm = math.Max(nrm, math.Abs(e[i*lde+j])+math.Abs(e[i*lde+j+1])) |
| } |
| enrmin = math.Min(enrmin, nrm) |
| enrmax = math.Max(enrmax, nrm) |
| ipair = 2 |
| case 2: |
| ipair = 0 |
| } |
| } |
| return math.Max(math.Abs(enrmin-1), math.Abs(enrmin-1)) |
| } |
| |
| // normalizeEV normalizes eigenvectors in the columns of E so that the element |
| // of largest magnitude has magnitude 1. |
| func normalizeEV(emat blas64.General, wi []float64) { |
| n := emat.Rows |
| if n == 0 { |
| return |
| } |
| var ( |
| bi = blas64.Implementation() |
| e = emat.Data |
| lde = emat.Stride |
| ipair int |
| ) |
| for j := 0; j < n; j++ { |
| if ipair == 0 && j < n-1 && wi[j] != 0 { |
| ipair = 1 |
| } |
| switch ipair { |
| case 0: |
| // Real eigenvector |
| ii := bi.Idamax(n, e[j:], lde) |
| remax := 1 / math.Abs(e[ii*lde+j]) |
| bi.Dscal(n, remax, e[j:], lde) |
| case 1: |
| // Complex eigenvector |
| var emax float64 |
| for i := 0; i < n; i++ { |
| emax = math.Max(emax, math.Abs(e[i*lde+j])+math.Abs(e[i*lde+j+1])) |
| } |
| bi.Dscal(n, 1/emax, e[j:], lde) |
| bi.Dscal(n, 1/emax, e[j+1:], lde) |
| ipair = 2 |
| case 2: |
| ipair = 0 |
| } |
| } |
| } |
