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

     1  // Copyright ©2020 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  // Dlangt returns the value of the given norm of an n×n tridiagonal matrix
    14  // represented by the three diagonals.
    15  //
    16  // d must have length at least n and dl and du must have length at least n-1.
    17  func (impl Implementation) Dlangt(norm lapack.MatrixNorm, n int, dl, d, du []float64) float64 {
    18  	switch {
    19  	case norm != lapack.MaxAbs && norm != lapack.MaxRowSum && norm != lapack.MaxColumnSum && norm != lapack.Frobenius:
    20  		panic(badNorm)
    21  	case n < 0:
    22  		panic(nLT0)
    23  	}
    24  
    25  	if n == 0 {
    26  		return 0
    27  	}
    28  
    29  	switch {
    30  	case len(dl) < n-1:
    31  		panic(shortDL)
    32  	case len(d) < n:
    33  		panic(shortD)
    34  	case len(du) < n-1:
    35  		panic(shortDU)
    36  	}
    37  
    38  	dl = dl[:n-1]
    39  	d = d[:n]
    40  	du = du[:n-1]
    41  
    42  	var anorm float64
    43  	switch norm {
    44  	case lapack.MaxAbs:
    45  		for _, diag := range [][]float64{dl, d, du} {
    46  			for _, di := range diag {
    47  				if math.IsNaN(di) {
    48  					return di
    49  				}
    50  				di = math.Abs(di)
    51  				if di > anorm {
    52  					anorm = di
    53  				}
    54  			}
    55  		}
    56  	case lapack.MaxColumnSum:
    57  		if n == 1 {
    58  			return math.Abs(d[0])
    59  		}
    60  		anorm = math.Abs(d[0]) + math.Abs(dl[0])
    61  		if math.IsNaN(anorm) {
    62  			return anorm
    63  		}
    64  		tmp := math.Abs(du[n-2]) + math.Abs(d[n-1])
    65  		if math.IsNaN(tmp) {
    66  			return tmp
    67  		}
    68  		if tmp > anorm {
    69  			anorm = tmp
    70  		}
    71  		for i := 1; i < n-1; i++ {
    72  			tmp = math.Abs(du[i-1]) + math.Abs(d[i]) + math.Abs(dl[i])
    73  			if math.IsNaN(tmp) {
    74  				return tmp
    75  			}
    76  			if tmp > anorm {
    77  				anorm = tmp
    78  			}
    79  		}
    80  	case lapack.MaxRowSum:
    81  		if n == 1 {
    82  			return math.Abs(d[0])
    83  		}
    84  		anorm = math.Abs(d[0]) + math.Abs(du[0])
    85  		if math.IsNaN(anorm) {
    86  			return anorm
    87  		}
    88  		tmp := math.Abs(dl[n-2]) + math.Abs(d[n-1])
    89  		if math.IsNaN(tmp) {
    90  			return tmp
    91  		}
    92  		if tmp > anorm {
    93  			anorm = tmp
    94  		}
    95  		for i := 1; i < n-1; i++ {
    96  			tmp = math.Abs(dl[i-1]) + math.Abs(d[i]) + math.Abs(du[i])
    97  			if math.IsNaN(tmp) {
    98  				return tmp
    99  			}
   100  			if tmp > anorm {
   101  				anorm = tmp
   102  			}
   103  		}
   104  	case lapack.Frobenius:
   105  		scale := 0.0
   106  		ssq := 1.0
   107  		scale, ssq = impl.Dlassq(n, d, 1, scale, ssq)
   108  		if n > 1 {
   109  			scale, ssq = impl.Dlassq(n-1, dl, 1, scale, ssq)
   110  			scale, ssq = impl.Dlassq(n-1, du, 1, scale, ssq)
   111  		}
   112  		anorm = scale * math.Sqrt(ssq)
   113  	}
   114  	return anorm
   115  }