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