diff --git a/lapack/gonum/dtrevc3.go b/lapack/gonum/dtrevc3.go index 4ed43ee4..c6568665 100644 --- a/lapack/gonum/dtrevc3.go +++ b/lapack/gonum/dtrevc3.go @@ -647,7 +647,7 @@ leftev: // when forming the right-hand side. beta := math.Max(norms[j], norms[j+1]) if beta > vcrit { - bi.Dscal(n-ki+1, 1/vmax, b[ki*ldb+iv:], 1) + bi.Dscal(n-ki, 1/vmax, b[ki*ldb+iv:], ldb) vmax = 1 } b[j*ldb+iv] -= bi.Ddot(j-ki-1, t[(ki+1)*ldt+j:], ldt, b[(ki+1)*ldb+iv:], ldb) diff --git a/lapack/testlapack/dgeev.go b/lapack/testlapack/dgeev.go index df68fdde..ad8894b9 100644 --- a/lapack/testlapack/dgeev.go +++ b/lapack/testlapack/dgeev.go @@ -783,9 +783,9 @@ func residualRightEV(a, e blas64.General, wr, wi []float64) float64 { return math.Min(errnorm/anorm, 1) } -// residualRightEV returns the residual +// residualLeftEV returns the residual // | Aᵀ E - E Wᵀ|_1 / ( |Aᵀ|_1 |E|_1 ) -// where the columns of E contain the right eigenvectors of A and W is a block diagonal matrix with +// where the columns of E contain the left eigenvectors of A and W is a block diagonal matrix with // a 1×1 block for each real eigenvalue and a 2×2 block for each complex conjugate pair. func residualLeftEV(a, e blas64.General, wr, wi []float64) float64 { // The implementation follows DGET22 routine from the Reference LAPACK's diff --git a/lapack/testlapack/dtrevc3.go b/lapack/testlapack/dtrevc3.go index 0715e6df..89dfe58b 100644 --- a/lapack/testlapack/dtrevc3.go +++ b/lapack/testlapack/dtrevc3.go @@ -11,8 +11,8 @@ import ( "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" ) @@ -22,202 +22,421 @@ type Dtrevc3er interface { func Dtrevc3Test(t *testing.T, impl Dtrevc3er) { rnd := rand.New(rand.NewSource(1)) + for _, side := range []lapack.EVSide{lapack.EVRight, lapack.EVLeft, lapack.EVBoth} { - for _, howmny := range []lapack.EVHowMany{lapack.EVAll, lapack.EVAllMulQ, lapack.EVSelected} { - for _, n := range []int{0, 1, 2, 3, 4, 5, 10, 34, 100} { - for _, extra := range []int{0, 11} { - for _, optwork := range []bool{true, false} { - for cas := 0; cas < 10; cas++ { - tmat := randomSchurCanonical(n, n+extra, rnd) - testDtrevc3(t, impl, side, howmny, tmat, optwork, rnd) - } - } + 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) } } } } } -func testDtrevc3(t *testing.T, impl Dtrevc3er, side lapack.EVSide, howmny lapack.EVHowMany, tmat blas64.General, optwork bool, rnd *rand.Rand) { - const tol = 1e-14 +// 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) - n := tmat.Rows - extra := tmat.Stride - tmat.Cols right := side != lapack.EVLeft left := side != lapack.EVRight - var selected, selectedWant []bool - var mWant int // How many columns will the eigenvectors occupy. - if howmny == lapack.EVSelected { - selected = make([]bool, n) - selectedWant = make([]bool, n) - // Dtrevc3 will compute only selected eigenvectors. Pick them - // randomly disregarding whether they are real or complex. - for i := range selected { - if rnd.Float64() < 0.5 { - selected[i] = true - } - } - // Dtrevc3 will modify (standardize) the slice selected based on - // whether the corresponding eigenvalues are real or complex. Do - // the same process here to fill selectedWant. - for i := 0; i < n; { - if i == n-1 || tmat.Data[(i+1)*tmat.Stride+i] == 0 { - // Real eigenvalue. - if selected[i] { - selectedWant[i] = true - mWant++ // Real eigenvectors occupy one column. - } - i++ - } else { - // Complex eigenvalue. - if selected[i] || selected[i+1] { - // Dtrevc3 will modify selected so that - // only the first element of the pair is - // true. - selectedWant[i] = true - mWant += 2 // Complex eigenvectors occupy two columns. - } - i += 2 - } - } - } else { - // All eigenvectors occupy n columns. - mWant = n - } + // 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, rnd) + tmatCopy := cloneGeneral(tmat) - var vr blas64.General + // 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 { - if howmny == lapack.EVAllMulQ { - vr = eye(n, n+extra) - } else { - // VR will be overwritten. - vr = nanGeneral(n, mWant, n+extra) - } + // Fill VR and VL with NaN because they should be completely overwritten in Dtrevc3. + vr = nanGeneral(n, n, n+extra) } - - var vl blas64.General if left { - if howmny == lapack.EVAllMulQ { - vl = eye(n, n+extra) - } else { - // VL will be overwritten. - vl = nanGeneral(n, mWant, n+extra) + 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) } } - work := make([]float64, max(1, 3*n)) + // 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, - vl.Data, max(1, vl.Stride), vr.Data, max(1, vr.Stride), mWant, work, -1) + vlSel.Data, max(1, vlSel.Stride), vrSel.Data, max(1, vrSel.Stride), mWant, work, -1) work = make([]float64, int(work[0])) } - m := impl.Dtrevc3(side, howmny, selected, n, tmat.Data, tmat.Stride, - vl.Data, max(1, vl.Stride), vr.Data, max(1, vr.Stride), mWant, work, len(work)) - - prefix := fmt.Sprintf("Case side=%v, howmny=%v, n=%v, extra=%v, optwk=%v", - side, howmny, n, extra, optwork) + 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", prefix) + t.Errorf("%v: out-of-range write to T", name) } - if !generalOutsideAllNaN(vl) { - t.Errorf("%v: out-of-range write to VL", prefix) + if !equalGeneral(tmat, tmatCopy) { + t.Errorf("%v: unexpected modification of T", name) } - if !generalOutsideAllNaN(vr) { - t.Errorf("%v: out-of-range write to VR", prefix) + 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 m != mWant { - t.Errorf("%v: unexpected value of m. Want %v, got %v", prefix, mWant, m) + if mGot != mWant { + t.Errorf("%v: unexpected value of selected m=%d, want %d", name, mGot, mWant) } - if howmny == lapack.EVSelected { - for i := range selected { - if selected[i] != selectedWant[i] { - t.Errorf("%v: unexpected selected[%v]", prefix, i) - } + for i := range selected { + if selected[i] != selectedWant[i] { + t.Errorf("%v: unexpected selected[%v]", name, i) } } - // Check that the columns of VR and VL are actually eigenvectors and - // that the magnitude of their largest element is 1. + // Check that selected columns of vrSel are equal to the corresponding + // columns of vr. var k int - for j := 0; j < n; { - re := tmat.Data[j*tmat.Stride+j] - if j == n-1 || tmat.Data[(j+1)*tmat.Stride+j] == 0 { - if howmny == lapack.EVSelected && !selected[j] { - j++ - continue - } - if right { - ev := columnOf(vr, k) - norm := floats.Norm(ev, math.Inf(1)) - if math.Abs(norm-1) > tol { - t.Errorf("%v: magnitude of largest element of VR[:,%v] not 1", prefix, k) + 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 + } } - if !isRightEigenvectorOf(tmat, ev, nil, complex(re, 0), tol) { - t.Errorf("%v: VR[:,%v] is not real right eigenvector", prefix, k) + 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 + } } - } - if left { - ev := columnOf(vl, k) - norm := floats.Norm(ev, math.Inf(1)) - if math.Abs(norm-1) > tol { - t.Errorf("%v: magnitude of largest element of VL[:,%v] not 1", prefix, k) - } - if !isLeftEigenvectorOf(tmat, ev, nil, complex(re, 0), tol) { - t.Errorf("%v: VL[:,%v] is not real left eigenvector", prefix, k) - } - } - k++ - j++ - continue - } - if howmny == lapack.EVSelected && !selected[j] { - j += 2 - continue - } - im := math.Sqrt(math.Abs(tmat.Data[(j+1)*tmat.Stride+j])) * - math.Sqrt(math.Abs(tmat.Data[j*tmat.Stride+j+1])) - if right { - evre := columnOf(vr, k) - evim := columnOf(vr, k+1) - var evmax float64 - for i, v := range evre { - evmax = math.Max(evmax, math.Abs(v)+math.Abs(evim[i])) - } - if math.Abs(evmax-1) > tol { - t.Errorf("%v: magnitude of largest element of VR[:,%v] not 1", prefix, k) - } - if !isRightEigenvectorOf(tmat, evre, evim, complex(re, im), tol) { - t.Errorf("%v: VR[:,%v:%v] is not complex right eigenvector", prefix, k, k+1) - } - floats.Scale(-1, evim) - if !isRightEigenvectorOf(tmat, evre, evim, complex(re, -im), tol) { - t.Errorf("%v: VR[:,%v:%v] is not complex right eigenvector", prefix, k, k+1) + k += 2 } } - if left { - evre := columnOf(vl, k) - evim := columnOf(vl, k+1) - var evmax float64 - for i, v := range evre { - evmax = math.Max(evmax, math.Abs(v)+math.Abs(evim[i])) - } - if math.Abs(evmax-1) > tol { - t.Errorf("%v: magnitude of largest element of VL[:,%v] not 1", prefix, k) - } - if !isLeftEigenvectorOf(tmat, evre, evim, complex(re, im), tol) { - t.Errorf("%v: VL[:,%v:%v] is not complex left eigenvector", prefix, k, k+1) - } - floats.Scale(-1, evim) - if !isLeftEigenvectorOf(tmat, evre, evim, complex(re, -im), tol) { - t.Errorf("%v: VL[:,%v:%v] is not complex left eigenvector", prefix, k, k+1) + } + 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 } } - k += 2 - j += 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 + } } }