github.com/gopherd/gonum@v0.0.4/lapack/gonum/dlansb.go (about)

     1  // Copyright ©2019 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/gopherd/gonum/blas"
    11  	"github.com/gopherd/gonum/lapack"
    12  )
    13  
    14  // Dlansb returns the given norm of an n×n symmetric band matrix with kd
    15  // super-diagonals.
    16  //
    17  // When norm is lapack.MaxColumnSum or lapack.MaxRowSum, the length of work must
    18  // be at least n.
    19  func (impl Implementation) Dlansb(norm lapack.MatrixNorm, uplo blas.Uplo, n, kd int, ab []float64, ldab int, work []float64) float64 {
    20  	switch {
    21  	case norm != lapack.MaxAbs && norm != lapack.MaxRowSum && norm != lapack.MaxColumnSum && norm != lapack.Frobenius:
    22  		panic(badNorm)
    23  	case uplo != blas.Upper && uplo != blas.Lower:
    24  		panic(badUplo)
    25  	case n < 0:
    26  		panic(nLT0)
    27  	case kd < 0:
    28  		panic(kdLT0)
    29  	case ldab < kd+1:
    30  		panic(badLdA)
    31  	}
    32  
    33  	// Quick return if possible.
    34  	if n == 0 {
    35  		return 0
    36  	}
    37  
    38  	switch {
    39  	case len(ab) < (n-1)*ldab+kd+1:
    40  		panic(shortAB)
    41  	case len(work) < n && (norm == lapack.MaxColumnSum || norm == lapack.MaxRowSum):
    42  		panic(shortWork)
    43  	}
    44  
    45  	var value float64
    46  	switch norm {
    47  	case lapack.MaxAbs:
    48  		if uplo == blas.Upper {
    49  			for i := 0; i < n; i++ {
    50  				for j := 0; j < min(n-i, kd+1); j++ {
    51  					aij := math.Abs(ab[i*ldab+j])
    52  					if aij > value || math.IsNaN(aij) {
    53  						value = aij
    54  					}
    55  				}
    56  			}
    57  		} else {
    58  			for i := 0; i < n; i++ {
    59  				for j := max(0, kd-i); j < kd+1; j++ {
    60  					aij := math.Abs(ab[i*ldab+j])
    61  					if aij > value || math.IsNaN(aij) {
    62  						value = aij
    63  					}
    64  				}
    65  			}
    66  		}
    67  	case lapack.MaxColumnSum, lapack.MaxRowSum:
    68  		work = work[:n]
    69  		var sum float64
    70  		if uplo == blas.Upper {
    71  			for i := range work {
    72  				work[i] = 0
    73  			}
    74  			for i := 0; i < n; i++ {
    75  				sum := work[i] + math.Abs(ab[i*ldab])
    76  				for j := i + 1; j < min(i+kd+1, n); j++ {
    77  					aij := math.Abs(ab[i*ldab+j-i])
    78  					sum += aij
    79  					work[j] += aij
    80  				}
    81  				if sum > value || math.IsNaN(sum) {
    82  					value = sum
    83  				}
    84  			}
    85  		} else {
    86  			for i := 0; i < n; i++ {
    87  				sum = 0
    88  				for j := max(0, i-kd); j < i; j++ {
    89  					aij := math.Abs(ab[i*ldab+kd+j-i])
    90  					sum += aij
    91  					work[j] += aij
    92  				}
    93  				work[i] = sum + math.Abs(ab[i*ldab+kd])
    94  			}
    95  			for _, sum := range work {
    96  				if sum > value || math.IsNaN(sum) {
    97  					value = sum
    98  				}
    99  			}
   100  		}
   101  	case lapack.Frobenius:
   102  		scale := 0.0
   103  		sum := 1.0
   104  		if uplo == blas.Upper {
   105  			if kd > 0 {
   106  				// Sum off-diagonals.
   107  				for i := 0; i < n-1; i++ {
   108  					ilen := min(n-i-1, kd)
   109  					scale, sum = impl.Dlassq(ilen, ab[i*ldab+1:], 1, scale, sum)
   110  				}
   111  				sum *= 2
   112  			}
   113  			// Sum diagonal.
   114  			scale, sum = impl.Dlassq(n, ab, ldab, scale, sum)
   115  		} else {
   116  			if kd > 0 {
   117  				// Sum off-diagonals.
   118  				for i := 1; i < n; i++ {
   119  					ilen := min(i, kd)
   120  					scale, sum = impl.Dlassq(ilen, ab[i*ldab+kd-ilen:], 1, scale, sum)
   121  				}
   122  				sum *= 2
   123  			}
   124  			// Sum diagonal.
   125  			scale, sum = impl.Dlassq(n, ab[kd:], ldab, scale, sum)
   126  		}
   127  		value = scale * math.Sqrt(sum)
   128  	}
   129  
   130  	return value
   131  }