github.com/jingcheng-WU/gonum@v0.9.1-0.20210323123734-f1a2a11a8f7b/lapack/gonum/dlascl.go (about)

     1  // Copyright ©2015 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  	"github.com/jingcheng-WU/gonum/lapack"
    11  )
    12  
    13  // Dlascl multiplies an m×n matrix by the scalar cto/cfrom.
    14  //
    15  // cfrom must not be zero, and cto and cfrom must not be NaN, otherwise Dlascl
    16  // will panic.
    17  //
    18  // Dlascl is an internal routine. It is exported for testing purposes.
    19  func (impl Implementation) Dlascl(kind lapack.MatrixType, kl, ku int, cfrom, cto float64, m, n int, a []float64, lda int) {
    20  	switch kind {
    21  	default:
    22  		panic(badMatrixType)
    23  	case 'H', 'B', 'Q', 'Z': // See dlascl.f.
    24  		panic("not implemented")
    25  	case lapack.General, lapack.UpperTri, lapack.LowerTri:
    26  		if lda < max(1, n) {
    27  			panic(badLdA)
    28  		}
    29  	}
    30  	switch {
    31  	case cfrom == 0:
    32  		panic(zeroCFrom)
    33  	case math.IsNaN(cfrom):
    34  		panic(nanCFrom)
    35  	case math.IsNaN(cto):
    36  		panic(nanCTo)
    37  	case m < 0:
    38  		panic(mLT0)
    39  	case n < 0:
    40  		panic(nLT0)
    41  	}
    42  
    43  	if n == 0 || m == 0 {
    44  		return
    45  	}
    46  
    47  	switch kind {
    48  	case lapack.General, lapack.UpperTri, lapack.LowerTri:
    49  		if len(a) < (m-1)*lda+n {
    50  			panic(shortA)
    51  		}
    52  	}
    53  
    54  	smlnum := dlamchS
    55  	bignum := 1 / smlnum
    56  	cfromc := cfrom
    57  	ctoc := cto
    58  	cfrom1 := cfromc * smlnum
    59  	for {
    60  		var done bool
    61  		var mul, ctol float64
    62  		if cfrom1 == cfromc {
    63  			// cfromc is inf.
    64  			mul = ctoc / cfromc
    65  			done = true
    66  			ctol = ctoc
    67  		} else {
    68  			ctol = ctoc / bignum
    69  			if ctol == ctoc {
    70  				// ctoc is either 0 or inf.
    71  				mul = ctoc
    72  				done = true
    73  				cfromc = 1
    74  			} else if math.Abs(cfrom1) > math.Abs(ctoc) && ctoc != 0 {
    75  				mul = smlnum
    76  				done = false
    77  				cfromc = cfrom1
    78  			} else if math.Abs(ctol) > math.Abs(cfromc) {
    79  				mul = bignum
    80  				done = false
    81  				ctoc = ctol
    82  			} else {
    83  				mul = ctoc / cfromc
    84  				done = true
    85  			}
    86  		}
    87  		switch kind {
    88  		case lapack.General:
    89  			for i := 0; i < m; i++ {
    90  				for j := 0; j < n; j++ {
    91  					a[i*lda+j] = a[i*lda+j] * mul
    92  				}
    93  			}
    94  		case lapack.UpperTri:
    95  			for i := 0; i < m; i++ {
    96  				for j := i; j < n; j++ {
    97  					a[i*lda+j] = a[i*lda+j] * mul
    98  				}
    99  			}
   100  		case lapack.LowerTri:
   101  			for i := 0; i < m; i++ {
   102  				for j := 0; j <= min(i, n-1); j++ {
   103  					a[i*lda+j] = a[i*lda+j] * mul
   104  				}
   105  			}
   106  		}
   107  		if done {
   108  			break
   109  		}
   110  	}
   111  }