lapack/gonum,testlapack: fix bug in Dtrevc3 and rework its test

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)

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

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} {
+		var name string
-			for _, n := range []int{0, 1, 2, 3, 4, 5, 10, 34, 100} {
+		switch side {
-				for _, extra := range []int{0, 11} {
+		case lapack.EVRight:
-					for _, optwork := range []bool{true, false} {
+			name = "Rigth"
-						for cas := 0; cas < 10; cas++ {
+		case lapack.EVLeft:
-							tmat := randomSchurCanonical(n, n+extra, rnd)
+			name = "Left"
-							testDtrevc3(t, impl, side, howmny, tmat, optwork, rnd)
+		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) {
+// dtrevc3Test tests Dtrevc3 by generating a random matrix T in Schur canonical
-	const tol = 1e-14
+// 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
+	// Generate a random matrix in Schur canonical form possibly with tiny or zero eigenvalues.
-	var mWant int // How many columns will the eigenvectors occupy.
+	// Zero elements of wi signify a real eigenvalue.
-	if howmny == lapack.EVSelected {
+	tmat, wr, wi := randomSchurCanonical(n, n+extra, rnd)
-		selected = make([]bool, n)
+	tmatCopy := cloneGeneral(tmat)
-		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
-	}
 
-	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 {
+		// Fill VR and VL with NaN because they should be completely overwritten in Dtrevc3.
-			vr = eye(n, n+extra)
+		vr = nanGeneral(n, n, n+extra)
-		} else {
-			// VR will be overwritten.
-			vr = nanGeneral(n, mWant, n+extra)
-		}
 	}
-
-	var vl blas64.General
 	if left {
-		if howmny == lapack.EVAllMulQ {
+		vl = nanGeneral(n, n, n+extra)
-			vl = eye(n, n+extra)
+	}
-		} else {
+
-			// VL will be overwritten.
+	var work []float64
-			vl = nanGeneral(n, mWant, n+extra)
+	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,
+	mGot = 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))
+		vlSel.Data, max(1, vlSel.Stride), vrSel.Data, max(1, vrSel.Stride), mWant, work, len(work))
-
-	prefix := fmt.Sprintf("Case side=%v, howmny=%v, n=%v, extra=%v, optwk=%v",
-		side, howmny, n, extra, optwork)
 
 	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) {
+	if !equalGeneral(tmat, tmatCopy) {
-		t.Errorf("%v: out-of-range write to VL", prefix)
+		t.Errorf("%v: unexpected modification of T", name)
 	}
-	if !generalOutsideAllNaN(vr) {
+	if !generalOutsideAllNaN(vrSel) {
-		t.Errorf("%v: out-of-range write to VR", prefix)
+		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 {
+	if mGot != mWant {
-		t.Errorf("%v: unexpected value of m. Want %v, got %v", prefix, mWant, m)
+		t.Errorf("%v: unexpected value of selected m=%d, want %d", name, mGot, mWant)
 	}
 
-	if howmny == lapack.EVSelected {
+	for i := range selected {
-		for i := range selected {
+		if selected[i] != selectedWant[i] {
-			if selected[i] != selectedWant[i] {
+			t.Errorf("%v: unexpected selected[%v]", name, i)
-				t.Errorf("%v: unexpected selected[%v]", prefix, i)
-			}
 		}
 	}
 
-	// Check that the columns of VR and VL are actually eigenvectors and
+	// Check that selected columns of vrSel are equal to the corresponding
-	// that the magnitude of their largest element is 1.
+	// columns of vr.
 	var k int
-	for j := 0; j < n; {
+	match := true
-		re := tmat.Data[j*tmat.Stride+j]
+	if right {
-		if j == n-1 || tmat.Data[(j+1)*tmat.Stride+j] == 0 {
+	loopVR:
-			if howmny == lapack.EVSelected && !selected[j] {
+		for j := 0; j < n; j++ {
-				j++
+			if selected[j] && wi[j] == 0 {
-				continue
+				for i := 0; i < n; i++ {
-			}
+					if vrSel.Data[i*vrSel.Stride+k] != vr.Data[i*vr.Stride+j] {
-			if right {
+						match = false
-				ev := columnOf(vr, k)
+						break loopVR
-				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)
 				}
-				if !isRightEigenvectorOf(tmat, ev, nil, complex(re, 0), tol) {
+				k++
-					t.Errorf("%v: VR[:,%v] is not real right eigenvector", prefix, 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 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)
 			}
 		}
-		if left {
+	}
-			evre := columnOf(vl, k)
+	if !match {
-			evim := columnOf(vl, k+1)
+		t.Errorf("%v: unexpected selected VR", name)
-			var evmax float64
+	}
-			for i, v := range evre {
+
-				evmax = math.Max(evmax, math.Abs(v)+math.Abs(evim[i]))
+	// Check that selected columns of vlSel are equal to the corresponding
-			}
+	// columns of vl.
-			if math.Abs(evmax-1) > tol {
+	match = true
-				t.Errorf("%v: magnitude of largest element of VL[:,%v] not 1", prefix, k)
+	k = 0
-			}
+	if left {
-			if !isLeftEigenvectorOf(tmat, evre, evim, complex(re, im), tol) {
+	loopVL:
-				t.Errorf("%v: VL[:,%v:%v] is not complex left eigenvector", prefix, k, k+1)
+		for j := 0; j < n; j++ {
-			}
+			if selected[j] && wi[j] == 0 {
-			floats.Scale(-1, evim)
+				for i := 0; i < n; i++ {
-			if !isLeftEigenvectorOf(tmat, evre, evim, complex(re, -im), tol) {
+					if vlSel.Data[i*vlSel.Stride+k] != vl.Data[i*vl.Stride+j] {
-				t.Errorf("%v: VL[:,%v:%v] is not complex left eigenvector", prefix, k, k+1)
+						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
+		}
 	}
 }