vendor/gonum.org/v1/gonum/lapack/gonum/dlantr.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 // Dlantr computes the specified norm of an m×n trapezoidal matrix A. If
  15 // norm == lapack.MaxColumnSum work must have length at least n, otherwise work
  16 // is unused.
  17 func (impl Implementation) Dlantr(norm lapack.MatrixNorm, uplo blas.Uplo, diag blas.Diag, m, n int, a []float64, lda int, work []float64) float64 {
  18         checkMatrix(m, n, a, lda)
  19         switch norm {
  20         case lapack.MaxRowSum, lapack.MaxColumnSum, lapack.NormFrob, lapack.MaxAbs:
  21         default:
  22                 panic(badNorm)
  23         }
  24         if uplo != blas.Upper && uplo != blas.Lower {
  25                 panic(badUplo)
  26         }
  27         if diag != blas.Unit && diag != blas.NonUnit {
  28                 panic(badDiag)
  29         }
  30         if norm == lapack.MaxColumnSum && len(work) < n {
  31                 panic(badWork)
  32         }
  33         if min(m, n) == 0 {
  34                 return 0
  35         }
  36         switch norm {
  37         default:
  38                 panic("unreachable")
  39         case lapack.MaxAbs:
  40                 if diag == blas.Unit {
  41                         value := 1.0
  42                         if uplo == blas.Upper {
  43                                 for i := 0; i < m; i++ {
  44                                         for j := i + 1; j < n; j++ {
  45                                                 tmp := math.Abs(a[i*lda+j])
  46                                                 if math.IsNaN(tmp) {
  47                                                         return tmp
  48                                                 }
  49                                                 if tmp > value {
  50                                                         value = tmp
  51                                                 }
  52                                         }
  53                                 }
  54                                 return value
  55                         }
  56                         for i := 1; i < m; i++ {
  57                                 for j := 0; j < min(i, n); j++ {
  58                                         tmp := math.Abs(a[i*lda+j])
  59                                         if math.IsNaN(tmp) {
  60                                                 return tmp
  61                                         }
  62                                         if tmp > value {
  63                                                 value = tmp
  64                                         }
  65                                 }
  66                         }
  67                         return value
  68                 }
  69                 var value float64
  70                 if uplo == blas.Upper {
  71                         for i := 0; i < m; i++ {
  72                                 for j := i; j < n; j++ {
  73                                         tmp := math.Abs(a[i*lda+j])
  74                                         if math.IsNaN(tmp) {
  75                                                 return tmp
  76                                         }
  77                                         if tmp > value {
  78                                                 value = tmp
  79                                         }
  80                                 }
  81                         }
  82                         return value
  83                 }
  84                 for i := 0; i < m; i++ {
  85                         for j := 0; j <= min(i, n-1); j++ {
  86                                 tmp := math.Abs(a[i*lda+j])
  87                                 if math.IsNaN(tmp) {
  88                                         return tmp
  89                                 }
  90                                 if tmp > value {
  91                                         value = tmp
  92                                 }
  93                         }
  94                 }
  95                 return value
  96         case lapack.MaxColumnSum:
  97                 if diag == blas.Unit {
  98                         for i := 0; i < min(m, n); i++ {
  99                                 work[i] = 1
 100                         }
 101                         for i := min(m, n); i < n; i++ {
 102                                 work[i] = 0
 103                         }
 104                         if uplo == blas.Upper {
 105                                 for i := 0; i < m; i++ {
 106                                         for j := i + 1; j < n; j++ {
 107                                                 work[j] += math.Abs(a[i*lda+j])
 108                                         }
 109                                 }
 110                         } else {
 111                                 for i := 1; i < m; i++ {
 112                                         for j := 0; j < min(i, n); j++ {
 113                                                 work[j] += math.Abs(a[i*lda+j])
 114                                         }
 115                                 }
 116                         }
 117                 } else {
 118                         for i := 0; i < n; i++ {
 119                                 work[i] = 0
 120                         }
 121                         if uplo == blas.Upper {
 122                                 for i := 0; i < m; i++ {
 123                                         for j := i; j < n; j++ {
 124                                                 work[j] += math.Abs(a[i*lda+j])
 125                                         }
 126                                 }
 127                         } else {
 128                                 for i := 0; i < m; i++ {
 129                                         for j := 0; j <= min(i, n-1); j++ {
 130                                                 work[j] += math.Abs(a[i*lda+j])
 131                                         }
 132                                 }
 133                         }
 134                 }
 135                 var max float64
 136                 for _, v := range work[:n] {
 137                         if math.IsNaN(v) {
 138                                 return math.NaN()
 139                         }
 140                         if v > max {
 141                                 max = v
 142                         }
 143                 }
 144                 return max
 145         case lapack.MaxRowSum:
 146                 var maxsum float64
 147                 if diag == blas.Unit {
 148                         if uplo == blas.Upper {
 149                                 for i := 0; i < m; i++ {
 150                                         var sum float64
 151                                         if i < min(m, n) {
 152                                                 sum = 1
 153                                         }
 154                                         for j := i + 1; j < n; j++ {
 155                                                 sum += math.Abs(a[i*lda+j])
 156                                         }
 157                                         if math.IsNaN(sum) {
 158                                                 return math.NaN()
 159                                         }
 160                                         if sum > maxsum {
 161                                                 maxsum = sum
 162                                         }
 163                                 }
 164                                 return maxsum
 165                         } else {
 166                                 for i := 1; i < m; i++ {
 167                                         var sum float64
 168                                         if i < min(m, n) {
 169                                                 sum = 1
 170                                         }
 171                                         for j := 0; j < min(i, n); j++ {
 172                                                 sum += math.Abs(a[i*lda+j])
 173                                         }
 174                                         if math.IsNaN(sum) {
 175                                                 return math.NaN()
 176                                         }
 177                                         if sum > maxsum {
 178                                                 maxsum = sum
 179                                         }
 180                                 }
 181                                 return maxsum
 182                         }
 183                 } else {
 184                         if uplo == blas.Upper {
 185                                 for i := 0; i < m; i++ {
 186                                         var sum float64
 187                                         for j := i; j < n; j++ {
 188                                                 sum += math.Abs(a[i*lda+j])
 189                                         }
 190                                         if math.IsNaN(sum) {
 191                                                 return sum
 192                                         }
 193                                         if sum > maxsum {
 194                                                 maxsum = sum
 195                                         }
 196                                 }
 197                                 return maxsum
 198                         } else {
 199                                 for i := 0; i < m; i++ {
 200                                         var sum float64
 201                                         for j := 0; j <= min(i, n-1); j++ {
 202                                                 sum += math.Abs(a[i*lda+j])
 203                                         }
 204                                         if math.IsNaN(sum) {
 205                                                 return sum
 206                                         }
 207                                         if sum > maxsum {
 208                                                 maxsum = sum
 209                                         }
 210                                 }
 211                                 return maxsum
 212                         }
 213                 }
 214         case lapack.NormFrob:
 215                 var nrm float64
 216                 if diag == blas.Unit {
 217                         if uplo == blas.Upper {
 218                                 for i := 0; i < m; i++ {
 219                                         for j := i + 1; j < n; j++ {
 220                                                 tmp := a[i*lda+j]
 221                                                 nrm += tmp * tmp
 222                                         }
 223                                 }
 224                         } else {
 225                                 for i := 1; i < m; i++ {
 226                                         for j := 0; j < min(i, n); j++ {
 227                                                 tmp := a[i*lda+j]
 228                                                 nrm += tmp * tmp
 229                                         }
 230                                 }
 231                         }
 232                         nrm += float64(min(m, n))
 233                 } else {
 234                         if uplo == blas.Upper {
 235                                 for i := 0; i < m; i++ {
 236                                         for j := i; j < n; j++ {
 237                                                 tmp := math.Abs(a[i*lda+j])
 238                                                 nrm += tmp * tmp
 239                                         }
 240                                 }
 241                         } else {
 242                                 for i := 0; i < m; i++ {
 243                                         for j := 0; j <= min(i, n-1); j++ {
 244                                                 tmp := math.Abs(a[i*lda+j])
 245                                                 nrm += tmp * tmp
 246                                         }
 247                                 }
 248                         }
 249                 }
 250                 return math.Sqrt(nrm)
 251         }
 252 }