github.com/gonum/lapack@v0.0.0-20181123203213-e4cdc5a0bff9/native/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 native
     6  
     7  import (
     8  	"math"
     9  
    10  	"github.com/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  	checkMatrix(m, n, a, lda)
    21  	if cfrom == 0 {
    22  		panic(zeroDiv)
    23  	}
    24  	if math.IsNaN(cfrom) || math.IsNaN(cto) {
    25  		panic(nanScale)
    26  	}
    27  	if n == 0 || m == 0 {
    28  		return
    29  	}
    30  	smlnum := dlamchS
    31  	bignum := 1 / smlnum
    32  	cfromc := cfrom
    33  	ctoc := cto
    34  	cfrom1 := cfromc * smlnum
    35  	for {
    36  		var done bool
    37  		var mul, ctol float64
    38  		if cfrom1 == cfromc {
    39  			// cfromc is inf.
    40  			mul = ctoc / cfromc
    41  			done = true
    42  			ctol = ctoc
    43  		} else {
    44  			ctol = ctoc / bignum
    45  			if ctol == ctoc {
    46  				// ctoc is either 0 or inf.
    47  				mul = ctoc
    48  				done = true
    49  				cfromc = 1
    50  			} else if math.Abs(cfrom1) > math.Abs(ctoc) && ctoc != 0 {
    51  				mul = smlnum
    52  				done = false
    53  				cfromc = cfrom1
    54  			} else if math.Abs(ctol) > math.Abs(cfromc) {
    55  				mul = bignum
    56  				done = false
    57  				ctoc = ctol
    58  			} else {
    59  				mul = ctoc / cfromc
    60  				done = true
    61  			}
    62  		}
    63  		switch kind {
    64  		default:
    65  			panic("lapack: not implemented")
    66  		case lapack.General:
    67  			for i := 0; i < m; i++ {
    68  				for j := 0; j < n; j++ {
    69  					a[i*lda+j] = a[i*lda+j] * mul
    70  				}
    71  			}
    72  		case lapack.UpperTri:
    73  			for i := 0; i < m; i++ {
    74  				for j := i; j < n; j++ {
    75  					a[i*lda+j] = a[i*lda+j] * mul
    76  				}
    77  			}
    78  		case lapack.LowerTri:
    79  			for i := 0; i < m; i++ {
    80  				for j := 0; j <= min(i, n-1); j++ {
    81  					a[i*lda+j] = a[i*lda+j] * mul
    82  				}
    83  			}
    84  		}
    85  		if done {
    86  			break
    87  		}
    88  	}
    89  }