github.com/gopherd/gonum@v0.0.4/lapack/gonum/dgebak.go (about) 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 "github.com/gopherd/gonum/blas/blas64" 9 "github.com/gopherd/gonum/lapack" 10 ) 11 12 // Dgebak updates an n×m matrix V as 13 // V = P D V if side == lapack.EVRight, 14 // V = P D^{-1} V if side == lapack.EVLeft, 15 // where P and D are n×n permutation and scaling matrices, respectively, 16 // implicitly represented by job, scale, ilo and ihi as returned by Dgebal. 17 // 18 // Typically, columns of the matrix V contain the right or left (determined by 19 // side) eigenvectors of the balanced matrix output by Dgebal, and Dgebak forms 20 // the eigenvectors of the original matrix. 21 // 22 // Dgebak is an internal routine. It is exported for testing purposes. 23 func (impl Implementation) Dgebak(job lapack.BalanceJob, side lapack.EVSide, n, ilo, ihi int, scale []float64, m int, v []float64, ldv int) { 24 switch { 25 case job != lapack.BalanceNone && job != lapack.Permute && job != lapack.Scale && job != lapack.PermuteScale: 26 panic(badBalanceJob) 27 case side != lapack.EVLeft && side != lapack.EVRight: 28 panic(badEVSide) 29 case n < 0: 30 panic(nLT0) 31 case ilo < 0 || max(0, n-1) < ilo: 32 panic(badIlo) 33 case ihi < min(ilo, n-1) || n <= ihi: 34 panic(badIhi) 35 case m < 0: 36 panic(mLT0) 37 case ldv < max(1, m): 38 panic(badLdV) 39 } 40 41 // Quick return if possible. 42 if n == 0 || m == 0 { 43 return 44 } 45 46 if len(scale) < n { 47 panic(shortScale) 48 } 49 if len(v) < (n-1)*ldv+m { 50 panic(shortV) 51 } 52 53 // Quick return if possible. 54 if job == lapack.BalanceNone { 55 return 56 } 57 58 bi := blas64.Implementation() 59 if ilo != ihi && job != lapack.Permute { 60 // Backward balance. 61 if side == lapack.EVRight { 62 for i := ilo; i <= ihi; i++ { 63 bi.Dscal(m, scale[i], v[i*ldv:], 1) 64 } 65 } else { 66 for i := ilo; i <= ihi; i++ { 67 bi.Dscal(m, 1/scale[i], v[i*ldv:], 1) 68 } 69 } 70 } 71 if job == lapack.Scale { 72 return 73 } 74 // Backward permutation. 75 for i := ilo - 1; i >= 0; i-- { 76 k := int(scale[i]) 77 if k == i { 78 continue 79 } 80 bi.Dswap(m, v[i*ldv:], 1, v[k*ldv:], 1) 81 } 82 for i := ihi + 1; i < n; i++ { 83 k := int(scale[i]) 84 if k == i { 85 continue 86 } 87 bi.Dswap(m, v[i*ldv:], 1, v[k*ldv:], 1) 88 } 89 }