1// Copyright ©2016 The Gonum Authors. All rights reserved.
2// Use of this source code is governed by a BSD-style
3// license that can be found in the LICENSE file.
4
5package gonum
6
7import (
8	"math"
9
10	"gonum.org/v1/gonum/lapack"
11)
12
13// Dsterf computes all eigenvalues of a symmetric tridiagonal matrix using the
14// Pal-Walker-Kahan variant of the QL or QR algorithm.
15//
16// d contains the diagonal elements of the tridiagonal matrix on entry, and
17// contains the eigenvalues in ascending order on exit. d must have length at
18// least n, or Dsterf will panic.
19//
20// e contains the off-diagonal elements of the tridiagonal matrix on entry, and is
21// overwritten during the call to Dsterf. e must have length of at least n-1 or
22// Dsterf will panic.
23//
24// Dsterf is an internal routine. It is exported for testing purposes.
25func (impl Implementation) Dsterf(n int, d, e []float64) (ok bool) {
26	if n < 0 {
27		panic(nLT0)
28	}
29
30	// Quick return if possible.
31	if n == 0 {
32		return true
33	}
34
35	switch {
36	case len(d) < n:
37		panic(shortD)
38	case len(e) < n-1:
39		panic(shortE)
40	}
41
42	if n == 1 {
43		return true
44	}
45
46	const (
47		none = 0 // The values are not scaled.
48		down = 1 // The values are scaled below ssfmax threshold.
49		up   = 2 // The values are scaled below ssfmin threshold.
50	)
51
52	// Determine the unit roundoff for this environment.
53	eps := dlamchE
54	eps2 := eps * eps
55	safmin := dlamchS
56	safmax := 1 / safmin
57	ssfmax := math.Sqrt(safmax) / 3
58	ssfmin := math.Sqrt(safmin) / eps2
59
60	// Compute the eigenvalues of the tridiagonal matrix.
61	maxit := 30
62	nmaxit := n * maxit
63	jtot := 0
64
65	l1 := 0
66
67	for {
68		if l1 > n-1 {
69			impl.Dlasrt(lapack.SortIncreasing, n, d)
70			return true
71		}
72		if l1 > 0 {
73			e[l1-1] = 0
74		}
75		var m int
76		for m = l1; m < n-1; m++ {
77			if math.Abs(e[m]) <= math.Sqrt(math.Abs(d[m]))*math.Sqrt(math.Abs(d[m+1]))*eps {
78				e[m] = 0
79				break
80			}
81		}
82
83		l := l1
84		lsv := l
85		lend := m
86		lendsv := lend
87		l1 = m + 1
88		if lend == 0 {
89			continue
90		}
91
92		// Scale submatrix in rows and columns l to lend.
93		anorm := impl.Dlanst(lapack.MaxAbs, lend-l+1, d[l:], e[l:])
94		iscale := none
95		if anorm == 0 {
96			continue
97		}
98		if anorm > ssfmax {
99			iscale = down
100			impl.Dlascl(lapack.General, 0, 0, anorm, ssfmax, lend-l+1, 1, d[l:], n)
101			impl.Dlascl(lapack.General, 0, 0, anorm, ssfmax, lend-l, 1, e[l:], n)
102		} else if anorm < ssfmin {
103			iscale = up
104			impl.Dlascl(lapack.General, 0, 0, anorm, ssfmin, lend-l+1, 1, d[l:], n)
105			impl.Dlascl(lapack.General, 0, 0, anorm, ssfmin, lend-l, 1, e[l:], n)
106		}
107
108		el := e[l:lend]
109		for i, v := range el {
110			el[i] *= v
111		}
112
113		// Choose between QL and QR iteration.
114		if math.Abs(d[lend]) < math.Abs(d[l]) {
115			lend = lsv
116			l = lendsv
117		}
118		if lend >= l {
119			// QL Iteration.
120			// Look for small sub-diagonal element.
121			for {
122				if l != lend {
123					for m = l; m < lend; m++ {
124						if math.Abs(e[m]) <= eps2*(math.Abs(d[m]*d[m+1])) {
125							break
126						}
127					}
128				} else {
129					m = lend
130				}
131				if m < lend {
132					e[m] = 0
133				}
134				p := d[l]
135				if m == l {
136					// Eigenvalue found.
137					l++
138					if l > lend {
139						break
140					}
141					continue
142				}
143				// If remaining matrix is 2 by 2, use Dlae2 to compute its eigenvalues.
144				if m == l+1 {
145					d[l], d[l+1] = impl.Dlae2(d[l], math.Sqrt(e[l]), d[l+1])
146					e[l] = 0
147					l += 2
148					if l > lend {
149						break
150					}
151					continue
152				}
153				if jtot == nmaxit {
154					break
155				}
156				jtot++
157
158				// Form shift.
159				rte := math.Sqrt(e[l])
160				sigma := (d[l+1] - p) / (2 * rte)
161				r := impl.Dlapy2(sigma, 1)
162				sigma = p - (rte / (sigma + math.Copysign(r, sigma)))
163
164				c := 1.0
165				s := 0.0
166				gamma := d[m] - sigma
167				p = gamma * gamma
168
169				// Inner loop.
170				for i := m - 1; i >= l; i-- {
171					bb := e[i]
172					r := p + bb
173					if i != m-1 {
174						e[i+1] = s * r
175					}
176					oldc := c
177					c = p / r
178					s = bb / r
179					oldgam := gamma
180					alpha := d[i]
181					gamma = c*(alpha-sigma) - s*oldgam
182					d[i+1] = oldgam + (alpha - gamma)
183					if c != 0 {
184						p = (gamma * gamma) / c
185					} else {
186						p = oldc * bb
187					}
188				}
189				e[l] = s * p
190				d[l] = sigma + gamma
191			}
192		} else {
193			for {
194				// QR Iteration.
195				// Look for small super-diagonal element.
196				for m = l; m > lend; m-- {
197					if math.Abs(e[m-1]) <= eps2*math.Abs(d[m]*d[m-1]) {
198						break
199					}
200				}
201				if m > lend {
202					e[m-1] = 0
203				}
204				p := d[l]
205				if m == l {
206					// Eigenvalue found.
207					l--
208					if l < lend {
209						break
210					}
211					continue
212				}
213
214				// If remaining matrix is 2 by 2, use Dlae2 to compute its eigenvalues.
215				if m == l-1 {
216					d[l], d[l-1] = impl.Dlae2(d[l], math.Sqrt(e[l-1]), d[l-1])
217					e[l-1] = 0
218					l -= 2
219					if l < lend {
220						break
221					}
222					continue
223				}
224				if jtot == nmaxit {
225					break
226				}
227				jtot++
228
229				// Form shift.
230				rte := math.Sqrt(e[l-1])
231				sigma := (d[l-1] - p) / (2 * rte)
232				r := impl.Dlapy2(sigma, 1)
233				sigma = p - (rte / (sigma + math.Copysign(r, sigma)))
234
235				c := 1.0
236				s := 0.0
237				gamma := d[m] - sigma
238				p = gamma * gamma
239
240				// Inner loop.
241				for i := m; i < l; i++ {
242					bb := e[i]
243					r := p + bb
244					if i != m {
245						e[i-1] = s * r
246					}
247					oldc := c
248					c = p / r
249					s = bb / r
250					oldgam := gamma
251					alpha := d[i+1]
252					gamma = c*(alpha-sigma) - s*oldgam
253					d[i] = oldgam + alpha - gamma
254					if c != 0 {
255						p = (gamma * gamma) / c
256					} else {
257						p = oldc * bb
258					}
259				}
260				e[l-1] = s * p
261				d[l] = sigma + gamma
262			}
263		}
264
265		// Undo scaling if necessary
266		switch iscale {
267		case down:
268			impl.Dlascl(lapack.General, 0, 0, ssfmax, anorm, lendsv-lsv+1, 1, d[lsv:], n)
269		case up:
270			impl.Dlascl(lapack.General, 0, 0, ssfmin, anorm, lendsv-lsv+1, 1, d[lsv:], n)
271		}
272
273		// Check for no convergence to an eigenvalue after a total of n*maxit iterations.
274		if jtot >= nmaxit {
275			break
276		}
277	}
278	for _, v := range e[:n-1] {
279		if v != 0 {
280			return false
281		}
282	}
283	impl.Dlasrt(lapack.SortIncreasing, n, d)
284	return true
285}
286