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/blas/blas64"
11	"gonum.org/v1/gonum/lapack"
12)
13
14// Dgebal balances an n×n matrix A. Balancing consists of two stages, permuting
15// and scaling. Both steps are optional and depend on the value of job.
16//
17// Permuting consists of applying a permutation matrix P such that the matrix
18// that results from P^T*A*P takes the upper block triangular form
19//            [ T1  X  Y  ]
20//  P^T A P = [  0  B  Z  ],
21//            [  0  0  T2 ]
22// where T1 and T2 are upper triangular matrices and B contains at least one
23// nonzero off-diagonal element in each row and column. The indices ilo and ihi
24// mark the starting and ending columns of the submatrix B. The eigenvalues of A
25// isolated in the first 0 to ilo-1 and last ihi+1 to n-1 elements on the
26// diagonal can be read off without any roundoff error.
27//
28// Scaling consists of applying a diagonal similarity transformation D such that
29// D^{-1}*B*D has the 1-norm of each row and its corresponding column nearly
30// equal. The output matrix is
31//  [ T1     X*D          Y    ]
32//  [  0  inv(D)*B*D  inv(D)*Z ].
33//  [  0      0           T2   ]
34// Scaling may reduce the 1-norm of the matrix, and improve the accuracy of
35// the computed eigenvalues and/or eigenvectors.
36//
37// job specifies the operations that will be performed on A.
38// If job is lapack.BalanceNone, Dgebal sets scale[i] = 1 for all i and returns ilo=0, ihi=n-1.
39// If job is lapack.Permute, only permuting will be done.
40// If job is lapack.Scale, only scaling will be done.
41// If job is lapack.PermuteScale, both permuting and scaling will be done.
42//
43// On return, if job is lapack.Permute or lapack.PermuteScale, it will hold that
44//  A[i,j] == 0,   for i > j and j ∈ {0, ..., ilo-1, ihi+1, ..., n-1}.
45// If job is lapack.BalanceNone or lapack.Scale, or if n == 0, it will hold that
46//  ilo == 0 and ihi == n-1.
47//
48// On return, scale will contain information about the permutations and scaling
49// factors applied to A. If π(j) denotes the index of the column interchanged
50// with column j, and D[j,j] denotes the scaling factor applied to column j,
51// then
52//  scale[j] == π(j),     for j ∈ {0, ..., ilo-1, ihi+1, ..., n-1},
53//           == D[j,j],   for j ∈ {ilo, ..., ihi}.
54// scale must have length equal to n, otherwise Dgebal will panic.
55//
56// Dgebal is an internal routine. It is exported for testing purposes.
57func (impl Implementation) Dgebal(job lapack.BalanceJob, n int, a []float64, lda int, scale []float64) (ilo, ihi int) {
58	switch {
59	case job != lapack.BalanceNone && job != lapack.Permute && job != lapack.Scale && job != lapack.PermuteScale:
60		panic(badBalanceJob)
61	case n < 0:
62		panic(nLT0)
63	case lda < max(1, n):
64		panic(badLdA)
65	}
66
67	ilo = 0
68	ihi = n - 1
69
70	if n == 0 {
71		return ilo, ihi
72	}
73
74	if len(scale) != n {
75		panic(shortScale)
76	}
77
78	if job == lapack.BalanceNone {
79		for i := range scale {
80			scale[i] = 1
81		}
82		return ilo, ihi
83	}
84
85	if len(a) < (n-1)*lda+n {
86		panic(shortA)
87	}
88
89	bi := blas64.Implementation()
90	swapped := true
91
92	if job == lapack.Scale {
93		goto scaling
94	}
95
96	// Permutation to isolate eigenvalues if possible.
97	//
98	// Search for rows isolating an eigenvalue and push them down.
99	for swapped {
100		swapped = false
101	rows:
102		for i := ihi; i >= 0; i-- {
103			for j := 0; j <= ihi; j++ {
104				if i == j {
105					continue
106				}
107				if a[i*lda+j] != 0 {
108					continue rows
109				}
110			}
111			// Row i has only zero off-diagonal elements in the
112			// block A[ilo:ihi+1,ilo:ihi+1].
113			scale[ihi] = float64(i)
114			if i != ihi {
115				bi.Dswap(ihi+1, a[i:], lda, a[ihi:], lda)
116				bi.Dswap(n, a[i*lda:], 1, a[ihi*lda:], 1)
117			}
118			if ihi == 0 {
119				scale[0] = 1
120				return ilo, ihi
121			}
122			ihi--
123			swapped = true
124			break
125		}
126	}
127	// Search for columns isolating an eigenvalue and push them left.
128	swapped = true
129	for swapped {
130		swapped = false
131	columns:
132		for j := ilo; j <= ihi; j++ {
133			for i := ilo; i <= ihi; i++ {
134				if i == j {
135					continue
136				}
137				if a[i*lda+j] != 0 {
138					continue columns
139				}
140			}
141			// Column j has only zero off-diagonal elements in the
142			// block A[ilo:ihi+1,ilo:ihi+1].
143			scale[ilo] = float64(j)
144			if j != ilo {
145				bi.Dswap(ihi+1, a[j:], lda, a[ilo:], lda)
146				bi.Dswap(n-ilo, a[j*lda+ilo:], 1, a[ilo*lda+ilo:], 1)
147			}
148			swapped = true
149			ilo++
150			break
151		}
152	}
153
154scaling:
155	for i := ilo; i <= ihi; i++ {
156		scale[i] = 1
157	}
158
159	if job == lapack.Permute {
160		return ilo, ihi
161	}
162
163	// Balance the submatrix in rows ilo to ihi.
164
165	const (
166		// sclfac should be a power of 2 to avoid roundoff errors.
167		// Elements of scale are restricted to powers of sclfac,
168		// therefore the matrix will be only nearly balanced.
169		sclfac = 2
170		// factor determines the minimum reduction of the row and column
171		// norms that is considered non-negligible. It must be less than 1.
172		factor = 0.95
173	)
174	sfmin1 := dlamchS / dlamchP
175	sfmax1 := 1 / sfmin1
176	sfmin2 := sfmin1 * sclfac
177	sfmax2 := 1 / sfmin2
178
179	// Iterative loop for norm reduction.
180	var conv bool
181	for !conv {
182		conv = true
183		for i := ilo; i <= ihi; i++ {
184			c := bi.Dnrm2(ihi-ilo+1, a[ilo*lda+i:], lda)
185			r := bi.Dnrm2(ihi-ilo+1, a[i*lda+ilo:], 1)
186			ica := bi.Idamax(ihi+1, a[i:], lda)
187			ca := math.Abs(a[ica*lda+i])
188			ira := bi.Idamax(n-ilo, a[i*lda+ilo:], 1)
189			ra := math.Abs(a[i*lda+ilo+ira])
190
191			// Guard against zero c or r due to underflow.
192			if c == 0 || r == 0 {
193				continue
194			}
195			g := r / sclfac
196			f := 1.0
197			s := c + r
198			for c < g && math.Max(f, math.Max(c, ca)) < sfmax2 && math.Min(r, math.Min(g, ra)) > sfmin2 {
199				if math.IsNaN(c + f + ca + r + g + ra) {
200					// Panic if NaN to avoid infinite loop.
201					panic("lapack: NaN")
202				}
203				f *= sclfac
204				c *= sclfac
205				ca *= sclfac
206				g /= sclfac
207				r /= sclfac
208				ra /= sclfac
209			}
210			g = c / sclfac
211			for r <= g && math.Max(r, ra) < sfmax2 && math.Min(math.Min(f, c), math.Min(g, ca)) > sfmin2 {
212				f /= sclfac
213				c /= sclfac
214				ca /= sclfac
215				g /= sclfac
216				r *= sclfac
217				ra *= sclfac
218			}
219
220			if c+r >= factor*s {
221				// Reduction would be negligible.
222				continue
223			}
224			if f < 1 && scale[i] < 1 && f*scale[i] <= sfmin1 {
225				continue
226			}
227			if f > 1 && scale[i] > 1 && scale[i] >= sfmax1/f {
228				continue
229			}
230
231			// Now balance.
232			scale[i] *= f
233			bi.Dscal(n-ilo, 1/f, a[i*lda+ilo:], 1)
234			bi.Dscal(ihi+1, f, a[i:], lda)
235			conv = false
236		}
237	}
238	return ilo, ihi
239}
240