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