gonum.org/v1/gonum@v0.15.1-0.20240517103525-f853624cb1bb/lapack/gonum/dlanhs.go (about)

     1  // Copyright ©2023 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  	"gonum.org/v1/gonum/blas/blas64"
    11  	"gonum.org/v1/gonum/lapack"
    12  )
    13  
    14  // Dlanhs returns the value of the one norm, or the Frobenius norm, or the
    15  // infinity norm, or the element of largest absolute value of a Hessenberg
    16  // matrix A.
    17  //
    18  // If norm is lapack.MaxColumnSum, work must have length at least n.
    19  func (impl Implementation) Dlanhs(norm lapack.MatrixNorm, n int, a []float64, lda int, work []float64) float64 {
    20  	switch {
    21  	case norm != lapack.MaxRowSum && norm != lapack.MaxAbs && norm != lapack.MaxColumnSum && norm != lapack.Frobenius:
    22  		panic(badNorm)
    23  	case n < 0:
    24  		panic(nLT0)
    25  	case lda < max(1, n):
    26  		panic(badLdA)
    27  	}
    28  
    29  	if n == 0 {
    30  		return 0
    31  	}
    32  
    33  	switch {
    34  	case len(a) < (n-1)*lda+n:
    35  		panic(shortA)
    36  	case norm == lapack.MaxColumnSum && len(work) < n:
    37  		panic(shortWork)
    38  	}
    39  
    40  	bi := blas64.Implementation()
    41  	var value float64
    42  	switch norm {
    43  	case lapack.MaxAbs:
    44  		for i := 0; i < n; i++ {
    45  			minj := max(0, i-1)
    46  			for _, v := range a[i*lda+minj : i*lda+n] {
    47  				value = math.Max(value, math.Abs(v))
    48  			}
    49  		}
    50  	case lapack.MaxColumnSum:
    51  		for i := 0; i < n; i++ {
    52  			work[i] = 0
    53  		}
    54  		for i := 0; i < n; i++ {
    55  			for j := max(0, i-1); j < n; j++ {
    56  				work[j] += math.Abs(a[i*lda+j])
    57  			}
    58  		}
    59  		for _, v := range work[:n] {
    60  			value = math.Max(value, v)
    61  		}
    62  	case lapack.MaxRowSum:
    63  		for i := 0; i < n; i++ {
    64  			minj := max(0, i-1)
    65  			sum := bi.Dasum(n-minj, a[i*lda+minj:], 1)
    66  			value = math.Max(value, sum)
    67  		}
    68  	case lapack.Frobenius:
    69  		scale := 0.0
    70  		sum := 1.0
    71  		for i := 0; i < n; i++ {
    72  			minj := max(0, i-1)
    73  			scale, sum = impl.Dlassq(n-minj, a[i*lda+minj:], 1, scale, sum)
    74  		}
    75  		value = scale * math.Sqrt(sum)
    76  	}
    77  	return value
    78  }