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