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

     1  // Copyright ©2021 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/internal/asm/f64"
    11  	"github.com/gopherd/gonum/lapack"
    12  )
    13  
    14  // Dlangb returns the given norm of an m×n band matrix with kl sub-diagonals and
    15  // ku super-diagonals.
    16  func (impl Implementation) Dlangb(norm lapack.MatrixNorm, m, n, kl, ku int, ab []float64, ldab int) float64 {
    17  	ncol := kl + 1 + ku
    18  	switch {
    19  	case norm != lapack.MaxAbs && norm != lapack.MaxRowSum && norm != lapack.MaxColumnSum && norm != lapack.Frobenius:
    20  		panic(badNorm)
    21  	case m < 0:
    22  		panic(mLT0)
    23  	case n < 0:
    24  		panic(nLT0)
    25  	case kl < 0:
    26  		panic(klLT0)
    27  	case ku < 0:
    28  		panic(kuLT0)
    29  	case ldab < ncol:
    30  		panic(badLdA)
    31  	}
    32  
    33  	// Quick return if possible.
    34  	if m == 0 || n == 0 {
    35  		return 0
    36  	}
    37  
    38  	switch {
    39  	case len(ab) < min(m, n+kl)*ldab:
    40  		panic(shortAB)
    41  	}
    42  
    43  	var value float64
    44  	switch norm {
    45  	case lapack.MaxAbs:
    46  		for i := 0; i < min(m, n+kl); i++ {
    47  			l := max(0, kl-i)
    48  			u := min(n+kl-i, ncol)
    49  			for _, aij := range ab[i*ldab+l : i*ldab+u] {
    50  				aij = math.Abs(aij)
    51  				if aij > value || math.IsNaN(aij) {
    52  					value = aij
    53  				}
    54  			}
    55  		}
    56  	case lapack.MaxRowSum:
    57  		for i := 0; i < min(m, n+kl); i++ {
    58  			l := max(0, kl-i)
    59  			u := min(n+kl-i, ncol)
    60  			sum := f64.L1Norm(ab[i*ldab+l : i*ldab+u])
    61  			if sum > value || math.IsNaN(sum) {
    62  				value = sum
    63  			}
    64  		}
    65  	case lapack.MaxColumnSum:
    66  		for j := 0; j < min(m+ku, n); j++ {
    67  			jb := min(kl+j, ncol-1)
    68  			ib := max(0, j-ku)
    69  			jlen := min(j+kl, m-1) - ib + 1
    70  			sum := f64.L1NormInc(ab[ib*ldab+jb:], jlen, max(1, ldab-1))
    71  			if sum > value || math.IsNaN(sum) {
    72  				value = sum
    73  			}
    74  		}
    75  	case lapack.Frobenius:
    76  		scale := 0.0
    77  		sum := 1.0
    78  		for i := 0; i < min(m, n+kl); i++ {
    79  			l := max(0, kl-i)
    80  			u := min(n+kl-i, ncol)
    81  			ilen := u - l
    82  			scale, sum = impl.Dlassq(ilen, ab[i*ldab+l:], 1, scale, sum)
    83  		}
    84  		value = scale * math.Sqrt(sum)
    85  	}
    86  	return value
    87  }