+++ /dev/null
-// Copyright ©2016 The Gonum Authors. All rights reserved.
-// Use of this source code is governed by a BSD-style
-// license that can be found in the LICENSE file.
-
-package gonum
-
-import (
- "math"
-
- "gonum.org/v1/gonum/blas/blas64"
- "gonum.org/v1/gonum/lapack"
-)
-
-// Dgebal balances an n×n matrix A. Balancing consists of two stages, permuting
-// and scaling. Both steps are optional and depend on the value of job.
-//
-// Permuting consists of applying a permutation matrix P such that the matrix
-// that results from P^T*A*P takes the upper block triangular form
-// [ T1 X Y ]
-// P^T A P = [ 0 B Z ],
-// [ 0 0 T2 ]
-// where T1 and T2 are upper triangular matrices and B contains at least one
-// nonzero off-diagonal element in each row and column. The indices ilo and ihi
-// mark the starting and ending columns of the submatrix B. The eigenvalues of A
-// isolated in the first 0 to ilo-1 and last ihi+1 to n-1 elements on the
-// diagonal can be read off without any roundoff error.
-//
-// Scaling consists of applying a diagonal similarity transformation D such that
-// D^{-1}*B*D has the 1-norm of each row and its corresponding column nearly
-// equal. The output matrix is
-// [ T1 X*D Y ]
-// [ 0 inv(D)*B*D inv(D)*Z ].
-// [ 0 0 T2 ]
-// Scaling may reduce the 1-norm of the matrix, and improve the accuracy of
-// the computed eigenvalues and/or eigenvectors.
-//
-// job specifies the operations that will be performed on A.
-// If job is lapack.None, Dgebal sets scale[i] = 1 for all i and returns ilo=0, ihi=n-1.
-// If job is lapack.Permute, only permuting will be done.
-// If job is lapack.Scale, only scaling will be done.
-// If job is lapack.PermuteScale, both permuting and scaling will be done.
-//
-// On return, if job is lapack.Permute or lapack.PermuteScale, it will hold that
-// A[i,j] == 0, for i > j and j ∈ {0, ..., ilo-1, ihi+1, ..., n-1}.
-// If job is lapack.None or lapack.Scale, or if n == 0, it will hold that
-// ilo == 0 and ihi == n-1.
-//
-// On return, scale will contain information about the permutations and scaling
-// factors applied to A. If π(j) denotes the index of the column interchanged
-// with column j, and D[j,j] denotes the scaling factor applied to column j,
-// then
-// scale[j] == π(j), for j ∈ {0, ..., ilo-1, ihi+1, ..., n-1},
-// == D[j,j], for j ∈ {ilo, ..., ihi}.
-// scale must have length equal to n, otherwise Dgebal will panic.
-//
-// Dgebal is an internal routine. It is exported for testing purposes.
-func (impl Implementation) Dgebal(job lapack.Job, n int, a []float64, lda int, scale []float64) (ilo, ihi int) {
- switch job {
- default:
- panic(badJob)
- case lapack.None, lapack.Permute, lapack.Scale, lapack.PermuteScale:
- }
- checkMatrix(n, n, a, lda)
- if len(scale) != n {
- panic("lapack: bad length of scale")
- }
-
- ilo = 0
- ihi = n - 1
-
- if n == 0 || job == lapack.None {
- for i := range scale {
- scale[i] = 1
- }
- return ilo, ihi
- }
-
- bi := blas64.Implementation()
- swapped := true
-
- if job == lapack.Scale {
- goto scaling
- }
-
- // Permutation to isolate eigenvalues if possible.
- //
- // Search for rows isolating an eigenvalue and push them down.
- for swapped {
- swapped = false
- rows:
- for i := ihi; i >= 0; i-- {
- for j := 0; j <= ihi; j++ {
- if i == j {
- continue
- }
- if a[i*lda+j] != 0 {
- continue rows
- }
- }
- // Row i has only zero off-diagonal elements in the
- // block A[ilo:ihi+1,ilo:ihi+1].
- scale[ihi] = float64(i)
- if i != ihi {
- bi.Dswap(ihi+1, a[i:], lda, a[ihi:], lda)
- bi.Dswap(n, a[i*lda:], 1, a[ihi*lda:], 1)
- }
- if ihi == 0 {
- scale[0] = 1
- return ilo, ihi
- }
- ihi--
- swapped = true
- break
- }
- }
- // Search for columns isolating an eigenvalue and push them left.
- swapped = true
- for swapped {
- swapped = false
- columns:
- for j := ilo; j <= ihi; j++ {
- for i := ilo; i <= ihi; i++ {
- if i == j {
- continue
- }
- if a[i*lda+j] != 0 {
- continue columns
- }
- }
- // Column j has only zero off-diagonal elements in the
- // block A[ilo:ihi+1,ilo:ihi+1].
- scale[ilo] = float64(j)
- if j != ilo {
- bi.Dswap(ihi+1, a[j:], lda, a[ilo:], lda)
- bi.Dswap(n-ilo, a[j*lda+ilo:], 1, a[ilo*lda+ilo:], 1)
- }
- swapped = true
- ilo++
- break
- }
- }
-
-scaling:
- for i := ilo; i <= ihi; i++ {
- scale[i] = 1
- }
-
- if job == lapack.Permute {
- return ilo, ihi
- }
-
- // Balance the submatrix in rows ilo to ihi.
-
- const (
- // sclfac should be a power of 2 to avoid roundoff errors.
- // Elements of scale are restricted to powers of sclfac,
- // therefore the matrix will be only nearly balanced.
- sclfac = 2
- // factor determines the minimum reduction of the row and column
- // norms that is considered non-negligible. It must be less than 1.
- factor = 0.95
- )
- sfmin1 := dlamchS / dlamchP
- sfmax1 := 1 / sfmin1
- sfmin2 := sfmin1 * sclfac
- sfmax2 := 1 / sfmin2
-
- // Iterative loop for norm reduction.
- var conv bool
- for !conv {
- conv = true
- for i := ilo; i <= ihi; i++ {
- c := bi.Dnrm2(ihi-ilo+1, a[ilo*lda+i:], lda)
- r := bi.Dnrm2(ihi-ilo+1, a[i*lda+ilo:], 1)
- ica := bi.Idamax(ihi+1, a[i:], lda)
- ca := math.Abs(a[ica*lda+i])
- ira := bi.Idamax(n-ilo, a[i*lda+ilo:], 1)
- ra := math.Abs(a[i*lda+ilo+ira])
-
- // Guard against zero c or r due to underflow.
- if c == 0 || r == 0 {
- continue
- }
- g := r / sclfac
- f := 1.0
- s := c + r
- for c < g && math.Max(f, math.Max(c, ca)) < sfmax2 && math.Min(r, math.Min(g, ra)) > sfmin2 {
- if math.IsNaN(c + f + ca + r + g + ra) {
- // Panic if NaN to avoid infinite loop.
- panic("lapack: NaN")
- }
- f *= sclfac
- c *= sclfac
- ca *= sclfac
- g /= sclfac
- r /= sclfac
- ra /= sclfac
- }
- g = c / sclfac
- for r <= g && math.Max(r, ra) < sfmax2 && math.Min(math.Min(f, c), math.Min(g, ca)) > sfmin2 {
- f /= sclfac
- c /= sclfac
- ca /= sclfac
- g /= sclfac
- r *= sclfac
- ra *= sclfac
- }
-
- if c+r >= factor*s {
- // Reduction would be negligible.
- continue
- }
- if f < 1 && scale[i] < 1 && f*scale[i] <= sfmin1 {
- continue
- }
- if f > 1 && scale[i] > 1 && scale[i] >= sfmax1/f {
- continue
- }
-
- // Now balance.
- scale[i] *= f
- bi.Dscal(n-ilo, 1/f, a[i*lda+ilo:], 1)
- bi.Dscal(ihi+1, f, a[i:], lda)
- conv = false
- }
- }
- return ilo, ihi
-}
OSDN Git Service
