vendor/gonum.org/v1/gonum/lapack/gonum/dlansy.go

   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         "gonum.org/v1/gonum/blas"
  11         "gonum.org/v1/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 }