github.com/jingcheng-WU/gonum@v0.9.1-0.20210323123734-f1a2a11a8f7b/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/jingcheng-WU/gonum/blas/blas64"
     9  	"github.com/jingcheng-WU/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  }