vendor/gonum.org/v1/gonum/lapack/gonum/dgebal.go

   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
   5 package gonum
   6
   7 import (
   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.None, 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.None 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.
  57 func (impl Implementation) Dgebal(job lapack.Job, n int, a []float64, lda int, scale []float64) (ilo, ihi int) {
  58         switch job {
  59         default:
  60                 panic(badJob)
  61         case lapack.None, lapack.Permute, lapack.Scale, lapack.PermuteScale:
  62         }
  63         checkMatrix(n, n, a, lda)
  64         if len(scale) != n {
  65                 panic("lapack: bad length of scale")
  66         }
  67
  68         ilo = 0
  69         ihi = n - 1
  70
  71         if n == 0 || job == lapack.None {
  72                 for i := range scale {
  73                         scale[i] = 1
  74                 }
  75                 return ilo, ihi
  76         }
  77
  78         bi := blas64.Implementation()
  79         swapped := true
  80
  81         if job == lapack.Scale {
  82                 goto scaling
  83         }
  84
  85         // Permutation to isolate eigenvalues if possible.
  86         //
  87         // Search for rows isolating an eigenvalue and push them down.
  88         for swapped {
  89                 swapped = false
  90         rows:
  91                 for i := ihi; i >= 0; i-- {
  92                         for j := 0; j <= ihi; j++ {
  93                                 if i == j {
  94                                         continue
  95                                 }
  96                                 if a[i*lda+j] != 0 {
  97                                         continue rows
  98                                 }
  99                         }
 100                         // Row i has only zero off-diagonal elements in the
 101                         // block A[ilo:ihi+1,ilo:ihi+1].
 102                         scale[ihi] = float64(i)
 103                         if i != ihi {
 104                                 bi.Dswap(ihi+1, a[i:], lda, a[ihi:], lda)
 105                                 bi.Dswap(n, a[i*lda:], 1, a[ihi*lda:], 1)
 106                         }
 107                         if ihi == 0 {
 108                                 scale[0] = 1
 109                                 return ilo, ihi
 110                         }
 111                         ihi--
 112                         swapped = true
 113                         break
 114                 }
 115         }
 116         // Search for columns isolating an eigenvalue and push them left.
 117         swapped = true
 118         for swapped {
 119                 swapped = false
 120         columns:
 121                 for j := ilo; j <= ihi; j++ {
 122                         for i := ilo; i <= ihi; i++ {
 123                                 if i == j {
 124                                         continue
 125                                 }
 126                                 if a[i*lda+j] != 0 {
 127                                         continue columns
 128                                 }
 129                         }
 130                         // Column j has only zero off-diagonal elements in the
 131                         // block A[ilo:ihi+1,ilo:ihi+1].
 132                         scale[ilo] = float64(j)
 133                         if j != ilo {
 134                                 bi.Dswap(ihi+1, a[j:], lda, a[ilo:], lda)
 135                                 bi.Dswap(n-ilo, a[j*lda+ilo:], 1, a[ilo*lda+ilo:], 1)
 136                         }
 137                         swapped = true
 138                         ilo++
 139                         break
 140                 }
 141         }
 142
 143 scaling:
 144         for i := ilo; i <= ihi; i++ {
 145                 scale[i] = 1
 146         }
 147
 148         if job == lapack.Permute {
 149                 return ilo, ihi
 150         }
 151
 152         // Balance the submatrix in rows ilo to ihi.
 153
 154         const (
 155                 // sclfac should be a power of 2 to avoid roundoff errors.
 156                 // Elements of scale are restricted to powers of sclfac,
 157                 // therefore the matrix will be only nearly balanced.
 158                 sclfac = 2
 159                 // factor determines the minimum reduction of the row and column
 160                 // norms that is considered non-negligible. It must be less than 1.
 161                 factor = 0.95
 162         )
 163         sfmin1 := dlamchS / dlamchP
 164         sfmax1 := 1 / sfmin1
 165         sfmin2 := sfmin1 * sclfac
 166         sfmax2 := 1 / sfmin2
 167
 168         // Iterative loop for norm reduction.
 169         var conv bool
 170         for !conv {
 171                 conv = true
 172                 for i := ilo; i <= ihi; i++ {
 173                         c := bi.Dnrm2(ihi-ilo+1, a[ilo*lda+i:], lda)
 174                         r := bi.Dnrm2(ihi-ilo+1, a[i*lda+ilo:], 1)
 175                         ica := bi.Idamax(ihi+1, a[i:], lda)
 176                         ca := math.Abs(a[ica*lda+i])
 177                         ira := bi.Idamax(n-ilo, a[i*lda+ilo:], 1)
 178                         ra := math.Abs(a[i*lda+ilo+ira])
 179
 180                         // Guard against zero c or r due to underflow.
 181                         if c == 0 || r == 0 {
 182                                 continue
 183                         }
 184                         g := r / sclfac
 185                         f := 1.0
 186                         s := c + r
 187                         for c < g && math.Max(f, math.Max(c, ca)) < sfmax2 && math.Min(r, math.Min(g, ra)) > sfmin2 {
 188                                 if math.IsNaN(c + f + ca + r + g + ra) {
 189                                         // Panic if NaN to avoid infinite loop.
 190                                         panic("lapack: NaN")
 191                                 }
 192                                 f *= sclfac
 193                                 c *= sclfac
 194                                 ca *= sclfac
 195                                 g /= sclfac
 196                                 r /= sclfac
 197                                 ra /= sclfac
 198                         }
 199                         g = c / sclfac
 200                         for r <= g && math.Max(r, ra) < sfmax2 && math.Min(math.Min(f, c), math.Min(g, ca)) > sfmin2 {
 201                                 f /= sclfac
 202                                 c /= sclfac
 203                                 ca /= sclfac
 204                                 g /= sclfac
 205                                 r *= sclfac
 206                                 ra *= sclfac
 207                         }
 208
 209                         if c+r >= factor*s {
 210                                 // Reduction would be negligible.
 211                                 continue
 212                         }
 213                         if f < 1 && scale[i] < 1 && f*scale[i] <= sfmin1 {
 214                                 continue
 215                         }
 216                         if f > 1 && scale[i] > 1 && scale[i] >= sfmax1/f {
 217                                 continue
 218                         }
 219
 220                         // Now balance.
 221                         scale[i] *= f
 222                         bi.Dscal(n-ilo, 1/f, a[i*lda+ilo:], 1)
 223                         bi.Dscal(ihi+1, f, a[i:], lda)
 224                         conv = false
 225                 }
 226         }
 227         return ilo, ihi
 228 }