diff --git a/native/dsteqr.go b/native/dsteqr.go index d408a7af..d4deefa0 100644 --- a/native/dsteqr.go +++ b/native/dsteqr.go @@ -88,7 +88,6 @@ func (impl Implementation) Dsteqr(compz lapack.EigComp, n int, d, e, z []float64 // Determine where the matrix splits and choose QL or QR iteration for each // block, according to whether top or bottom diagonal element is smaller. - // TODO(btracey): Reformat the gotos to use structural flow techniques. // TODO(btracey): Move these variables closer to actual usage when // gotos are removed. l1 := 0 @@ -104,256 +103,268 @@ func (impl Implementation) Dsteqr(compz lapack.EigComp, n int, d, e, z []float64 ) var iscale scaletype -Ten: - if l1 > n-1 { - goto OneSixty - } - if l1 > 0 { - e[l1-1] = 0 - } - if l1 <= nm1 { - for m = l1; m < nm1; m++ { - test := math.Abs(e[m]) - if test == 0 { - goto Thirty - } - if test <= (math.Sqrt(math.Abs(d[m]))*math.Sqrt(math.Abs(d[m+1])))*eps { - e[m] = 0 - goto Thirty + for { + if l1 > n-1 { + // Order eigenvalues and eigenvectors. + if icompz == 0 { + impl.Dlasrt(lapack.SortIncreasing, n, d) + } else { + // TODO(btracey): Consider replacing this sort with a call to sort.Sort. + for ii := 1; ii < n; ii++ { + i := ii - 1 + k := i + p := d[i] + for j := ii; j < n; j++ { + if d[j] < p { + k = j + p = d[j] + } + } + if k != i { + d[k] = d[i] + d[i] = p + bi.Dswap(n, z[i:], ldz, z[k:], ldz) + } + } } + return true } - } - m = n - 1 -Thirty: - l = l1 - lsv = l - lend = m - lendsv = lend - l1 = m + 1 - if lend == l { - goto Ten - } - - // Scale submatrix in rows and columns L to Lend - anorm = impl.Dlanst(lapack.MaxAbs, lend-l+1, d[l:], e[l:]) - switch { - case anorm == 0: - goto Ten - case anorm > ssfmax: - iscale = down - // TODO(btracey): Why is lda n? - impl.Dlascl(lapack.General, 0, 0, anorm, ssfmax, lend-l+1, 1, d[l:], n) - impl.Dlascl(lapack.General, 0, 0, anorm, ssfmax, lend-l, 1, e[l:], n) - case anorm < ssfmin: - iscale = up - // TODO(btracey): Why is lda n? - impl.Dlascl(lapack.General, 0, 0, anorm, ssfmin, lend-l+1, 1, d[l:], n) - impl.Dlascl(lapack.General, 0, 0, anorm, ssfmin, lend-l, 1, e[l:], n) - } - - // Choose between QL and QR. - if math.Abs(d[lend]) < math.Abs(d[l]) { - lend = lsv - l = lendsv - } - if lend > l { - // QL Iteration. Look for small subdiagonal element. - Fourty: - if l != lend { - lendm1 := lend - 1 - for m = l; m <= lendm1; m++ { - v := math.Abs(e[m]) - if v*v <= (eps2*math.Abs(d[m]))*math.Abs(d[m+1])+safmin { - goto Sixty + if l1 > 0 { + e[l1-1] = 0 + } + if l1 <= nm1 { + for m = l1; m < nm1; m++ { + test := math.Abs(e[m]) + if test == 0 { + break + } + if test <= (math.Sqrt(math.Abs(d[m]))*math.Sqrt(math.Abs(d[m+1])))*eps { + e[m] = 0 + break } } } - m = lend - Sixty: - if m < lend { - e[m] = 0 - } - p = d[l] - if m == l { - goto Eighty + l = l1 + lsv = l + lend = m + lendsv = lend + l1 = m + 1 + if lend == l { + continue } - // If remaining matrix is 2×2, use Dlae2 to compute its eigensystem. - if m == l+1 { - if icompz > 0 { - d[l], d[l+1], work[l], work[n-1+l] = impl.Dlaev2(d[l], e[l], d[l+1]) - impl.Dlasr(blas.Right, lapack.Variable, lapack.Backward, - n, 2, work[l:], work[n-1+l:], z[l:], ldz) - } else { - d[l], d[l+1] = impl.Dlae2(d[l], e[l], d[l+1]) - } - e[l] = 0 - l += 2 - if l <= lend { - goto Fourty - } - goto OneFourty + // Scale submatrix in rows and columns L to Lend + anorm = impl.Dlanst(lapack.MaxAbs, lend-l+1, d[l:], e[l:]) + switch { + case anorm == 0: + continue + case anorm > ssfmax: + iscale = down + // TODO(btracey): Why is lda n? + impl.Dlascl(lapack.General, 0, 0, anorm, ssfmax, lend-l+1, 1, d[l:], n) + impl.Dlascl(lapack.General, 0, 0, anorm, ssfmax, lend-l, 1, e[l:], n) + case anorm < ssfmin: + iscale = up + // TODO(btracey): Why is lda n? + impl.Dlascl(lapack.General, 0, 0, anorm, ssfmin, lend-l+1, 1, d[l:], n) + impl.Dlascl(lapack.General, 0, 0, anorm, ssfmin, lend-l, 1, e[l:], n) } - if jtot == nmaxit { - goto OneFourty + // Choose between QL and QR. + if math.Abs(d[lend]) < math.Abs(d[l]) { + lend = lsv + l = lendsv } - jtot++ - - // Form shift - g = (d[l+1] - p) / (2 * e[l]) - r = impl.Dlapy2(g, 1) - g = d[m] - p + e[l]/(g+math.Copysign(r, g)) - s = 1 - c = 1 - p = 0 - - // Inner loop - for i := m - 1; i >= l; i-- { - f := s * e[i] - b := c * e[i] - c, s, r = impl.Dlartg(g, f) - if i != m-1 { - e[i+1] = r - } - g = d[i+1] - p - r = (d[i]-g)*s + 2*c*b - p = s * r - d[i+1] = g + p - g = c*r - b - - // If eigenvectors are desired, then save rotations. - if icompz > 0 { - work[i] = c - work[n-1+i] = -s - } - } - // If eigenvectors are desired, then apply saved rotations. - if icompz > 0 { - mm := m - l + 1 - impl.Dlasr(blas.Right, lapack.Variable, lapack.Backward, - n, mm, work[l:], work[n-1+l:], z[l:], ldz) - } - d[l] -= p - e[l] = g - goto Fourty - - // Eigenvalue found. - Eighty: - d[l] = p - l++ - if l <= lend { - goto Fourty - } - goto OneFourty - } else { - // QR Iteration. - // Look for small superdiagonal element. - Ninety: - if l != lend { - lendp1 := lend + 1 - for m = l; m >= lendp1; m-- { - v := math.Abs(e[m-1]) - if v*v <= (eps2*math.Abs(d[m])*math.Abs(d[m-1]) + safmin) { - goto OneTen + if lend > l { + // QL Iteration. Look for small subdiagonal element. + for { + if l != lend { + for m = l; m < lend; m++ { + v := math.Abs(e[m]) + if v*v <= (eps2*math.Abs(d[m]))*math.Abs(d[m+1])+safmin { + break + } + } + } else { + m = lend } + if m < lend { + e[m] = 0 + } + p = d[l] + if m == l { + // Eigenvalue found. + d[l] = p + l++ + if l > lend { + break + } + continue + } + + // If remaining matrix is 2×2, use Dlae2 to compute its eigensystem. + if m == l+1 { + if icompz > 0 { + d[l], d[l+1], work[l], work[n-1+l] = impl.Dlaev2(d[l], e[l], d[l+1]) + impl.Dlasr(blas.Right, lapack.Variable, lapack.Backward, + n, 2, work[l:], work[n-1+l:], z[l:], ldz) + } else { + d[l], d[l+1] = impl.Dlae2(d[l], e[l], d[l+1]) + } + e[l] = 0 + l += 2 + if l > lend { + break + } + continue + } + + if jtot == nmaxit { + break + } + jtot++ + + // Form shift + g = (d[l+1] - p) / (2 * e[l]) + r = impl.Dlapy2(g, 1) + g = d[m] - p + e[l]/(g+math.Copysign(r, g)) + s = 1 + c = 1 + p = 0 + + // Inner loop + for i := m - 1; i >= l; i-- { + f := s * e[i] + b := c * e[i] + c, s, r = impl.Dlartg(g, f) + if i != m-1 { + e[i+1] = r + } + g = d[i+1] - p + r = (d[i]-g)*s + 2*c*b + p = s * r + d[i+1] = g + p + g = c*r - b + + // If eigenvectors are desired, then save rotations. + if icompz > 0 { + work[i] = c + work[n-1+i] = -s + } + } + // If eigenvectors are desired, then apply saved rotations. + if icompz > 0 { + mm := m - l + 1 + impl.Dlasr(blas.Right, lapack.Variable, lapack.Backward, + n, mm, work[l:], work[n-1+l:], z[l:], ldz) + } + d[l] -= p + e[l] = g } - } - m = lend - OneTen: - if m > lend { - e[m-1] = 0 - } - p = d[l] - if m == l { - goto OneThirty - } + } else { + // QR Iteration. + // Look for small superdiagonal element. + for { + if l != lend { + for m = l; m > lend; m-- { + v := math.Abs(e[m-1]) + if v*v <= (eps2*math.Abs(d[m])*math.Abs(d[m-1]) + safmin) { + break + } + } + } else { + m = lend + } + if m > lend { + e[m-1] = 0 + } + p = d[l] + if m == l { + // Eigenvalue found + d[l] = p + l-- + if l < lend { + break + } + continue + } - // If remaining matrix is 2×2, use Dlae2 to compute its eigenvalues. - if m == l-1 { - if icompz > 0 { - d[l-1], d[l], work[m], work[n-1+m] = impl.Dlaev2(d[l-1], e[l-1], d[l]) - impl.Dlasr(blas.Right, lapack.Variable, lapack.Forward, - n, 2, work[m:], work[n-1+m:], z[l-1:], ldz) - } else { - d[l-1], d[l] = impl.Dlae2(d[l-1], e[l-1], d[l]) - } - e[l-1] = 0 - l -= 2 - if l >= lend { - goto Ninety - } - goto OneFourty - } - if jtot == nmaxit { - goto OneFourty - } - jtot++ + // If remaining matrix is 2×2, use Dlae2 to compute its eigenvalues. + if m == l-1 { + if icompz > 0 { + d[l-1], d[l], work[m], work[n-1+m] = impl.Dlaev2(d[l-1], e[l-1], d[l]) + impl.Dlasr(blas.Right, lapack.Variable, lapack.Forward, + n, 2, work[m:], work[n-1+m:], z[l-1:], ldz) + } else { + d[l-1], d[l] = impl.Dlae2(d[l-1], e[l-1], d[l]) + } + e[l-1] = 0 + l -= 2 + if l < lend { + break + } + continue + } + if jtot == nmaxit { + break + } + jtot++ - // Form shift. - g = (d[l-1] - p) / (2 * e[l-1]) - r = impl.Dlapy2(g, 1) - g = d[m] - p + (e[l-1])/(g+math.Copysign(r, g)) + // Form shift. + g = (d[l-1] - p) / (2 * e[l-1]) + r = impl.Dlapy2(g, 1) + g = d[m] - p + (e[l-1])/(g+math.Copysign(r, g)) - s = 1 - c = 1 - p = 0 + s = 1 + c = 1 + p = 0 - // Inner loop. - for i := m; i < l; i++ { - f := s * e[i] - b := c * e[i] - c, s, r = impl.Dlartg(g, f) - if i != m { - e[i-1] = r - } - g = d[i] - p - r = (d[i+1]-g)*s + 2*c*b - p = s * r - d[i] = g + p - g = c*r - b + // Inner loop. + for i := m; i < l; i++ { + f := s * e[i] + b := c * e[i] + c, s, r = impl.Dlartg(g, f) + if i != m { + e[i-1] = r + } + g = d[i] - p + r = (d[i+1]-g)*s + 2*c*b + p = s * r + d[i] = g + p + g = c*r - b - // If eigenvectors are desired, then save rotations. - if icompz > 0 { - work[i] = c - work[n-1+i] = s + // If eigenvectors are desired, then save rotations. + if icompz > 0 { + work[i] = c + work[n-1+i] = s + } + } + + // If eigenvectors are desired, then apply saved rotations. + if icompz > 0 { + mm := l - m + 1 + impl.Dlasr(blas.Right, lapack.Variable, lapack.Forward, + n, mm, work[m:], work[n-1+m:], z[m:], ldz) + } + d[l] -= p + e[l-1] = g } } - // If eigenvectors are desired, then apply saved rotations. - if icompz > 0 { - mm := l - m + 1 - impl.Dlasr(blas.Right, lapack.Variable, lapack.Forward, - n, mm, work[m:], work[n-1+m:], z[m:], ldz) + // Undo scaling if necessary. + switch iscale { + case down: + impl.Dlascl(lapack.General, 0, 0, ssfmax, anorm, lendsv-lsv+1, 1, d[lsv:], n) + impl.Dlascl(lapack.General, 0, 0, ssfmax, anorm, lendsv-lsv, 1, e[lsv:], n) + case up: + impl.Dlascl(lapack.General, 0, 0, ssfmin, anorm, lendsv-lsv+1, 1, e[lsv:], n) + impl.Dlascl(lapack.General, 0, 0, ssfmin, anorm, lendsv-lsv, 1, e[lsv:], n) } - d[l] -= p - e[l-1] = g - goto Ninety - // Eigenvalue found - OneThirty: - d[l] = p - l-- - if l >= lend { - goto Ninety + // Check for no convergence to an eigenvalue after a total of n*maxit iterations. + if jtot >= nmaxit { + break } - goto OneFourty - } - - // Undo scaling if necessary. -OneFourty: - switch iscale { - case down: - impl.Dlascl(lapack.General, 0, 0, ssfmax, anorm, lendsv-lsv+1, 1, d[lsv:], n) - impl.Dlascl(lapack.General, 0, 0, ssfmax, anorm, lendsv-lsv, 1, e[lsv:], n) - case up: - impl.Dlascl(lapack.General, 0, 0, ssfmin, anorm, lendsv-lsv+1, 1, e[lsv:], n) - impl.Dlascl(lapack.General, 0, 0, ssfmin, anorm, lendsv-lsv, 1, e[lsv:], n) - } - - // Check for no convergence to an eigenvalue after a total of n*maxit iterations. - if jtot < nmaxit { - goto Ten } for i := 0; i < n-1; i++ { if e[i] != 0 { @@ -361,29 +372,4 @@ OneFourty: } } return true - - // Order eigenvalues and eigenvectors. -OneSixty: - if icompz == 0 { - impl.Dlasrt(lapack.SortIncreasing, n, d) - } else { - // TODO(btracey): Consider replacing this sort with a call to sort.Sort. - for ii := 1; ii < n; ii++ { - i := ii - 1 - k := i - p := d[i] - for j := ii; j < n; j++ { - if d[j] < p { - k = j - p = d[j] - } - } - if k != i { - d[k] = d[i] - d[i] = p - bi.Dswap(n, z[i:], ldz, z[k:], ldz) - } - } - } - return true }