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