blob: f9ebe1ba25080723ce0126dd2833eb06a7acc8df [file] [log] [blame]
 // Copyright ©2020 The Gonum Authors. All rights reserved. // Use of this source code is governed by a BSD-style // license that can be found in the LICENSE file. package gonum import ( "math" "gonum.org/v1/gonum/blas" "gonum.org/v1/gonum/lapack" ) // Dlantb returns the value of the given norm of an n×n triangular band matrix A // with k+1 diagonals. // // When norm is lapack.MaxColumnSum, the length of work must be at least n. func (impl Implementation) Dlantb(norm lapack.MatrixNorm, uplo blas.Uplo, diag blas.Diag, n, k int, a []float64, lda int, work []float64) float64 { switch { case norm != lapack.MaxAbs && norm != lapack.MaxRowSum && norm != lapack.MaxColumnSum && norm != lapack.Frobenius: panic(badNorm) case uplo != blas.Upper && uplo != blas.Lower: panic(badUplo) case n < 0: panic(nLT0) case k < 0: panic(kdLT0) case lda < k+1: panic(badLdA) } // Quick return if possible. if n == 0 { return 0 } switch { case len(a) < (n-1)*lda+k+1: panic(shortAB) case len(work) < n && norm == lapack.MaxColumnSum: panic(shortWork) } var value float64 switch norm { case lapack.MaxAbs: if uplo == blas.Upper { var jfirst int if diag == blas.Unit { value = 1 jfirst = 1 } for i := 0; i < n; i++ { for _, aij := range a[i*lda+jfirst : i*lda+min(n-i, k+1)] { if math.IsNaN(aij) { return aij } aij = math.Abs(aij) if aij > value { value = aij } } } } else { jlast := k + 1 if diag == blas.Unit { value = 1 jlast = k } for i := 0; i < n; i++ { for _, aij := range a[i*lda+max(0, k-i) : i*lda+jlast] { if math.IsNaN(aij) { return math.NaN() } aij = math.Abs(aij) if aij > value { value = aij } } } } case lapack.MaxRowSum: var sum float64 if uplo == blas.Upper { var jfirst int if diag == blas.Unit { jfirst = 1 } for i := 0; i < n; i++ { sum = 0 if diag == blas.Unit { sum = 1 } for _, aij := range a[i*lda+jfirst : i*lda+min(n-i, k+1)] { sum += math.Abs(aij) } if math.IsNaN(sum) { return math.NaN() } if sum > value { value = sum } } } else { jlast := k + 1 if diag == blas.Unit { jlast = k } for i := 0; i < n; i++ { sum = 0 if diag == blas.Unit { sum = 1 } for _, aij := range a[i*lda+max(0, k-i) : i*lda+jlast] { sum += math.Abs(aij) } if math.IsNaN(sum) { return math.NaN() } if sum > value { value = sum } } } case lapack.MaxColumnSum: work = work[:n] if diag == blas.Unit { for i := range work { work[i] = 1 } } else { for i := range work { work[i] = 0 } } if uplo == blas.Upper { var jfirst int if diag == blas.Unit { jfirst = 1 } for i := 0; i < n; i++ { for j, aij := range a[i*lda+jfirst : i*lda+min(n-i, k+1)] { work[i+jfirst+j] += math.Abs(aij) } } } else { jlast := k + 1 if diag == blas.Unit { jlast = k } for i := 0; i < n; i++ { off := max(0, k-i) for j, aij := range a[i*lda+off : i*lda+jlast] { work[i+j+off-k] += math.Abs(aij) } } } for _, wi := range work { if math.IsNaN(wi) { return math.NaN() } if wi > value { value = wi } } case lapack.Frobenius: var scale, ssq float64 switch uplo { case blas.Upper: if diag == blas.Unit { scale = 1 ssq = float64(n) if k > 0 { for i := 0; i < n-1; i++ { ilen := min(n-i-1, k) rowscale, rowssq := impl.Dlassq(ilen, a[i*lda+1:], 1, 0, 1) scale, ssq = impl.Dcombssq(scale, ssq, rowscale, rowssq) } } } else { scale = 0 ssq = 1 for i := 0; i < n; i++ { ilen := min(n-i, k+1) rowscale, rowssq := impl.Dlassq(ilen, a[i*lda:], 1, 0, 1) scale, ssq = impl.Dcombssq(scale, ssq, rowscale, rowssq) } } case blas.Lower: if diag == blas.Unit { scale = 1 ssq = float64(n) if k > 0 { for i := 1; i < n; i++ { ilen := min(i, k) rowscale, rowssq := impl.Dlassq(ilen, a[i*lda+k-ilen:], 1, 0, 1) scale, ssq = impl.Dcombssq(scale, ssq, rowscale, rowssq) } } } else { scale = 0 ssq = 1 for i := 0; i < n; i++ { ilen := min(i, k) + 1 rowscale, rowssq := impl.Dlassq(ilen, a[i*lda+k+1-ilen:], 1, 0, 1) scale, ssq = impl.Dcombssq(scale, ssq, rowscale, rowssq) } } } value = scale * math.Sqrt(ssq) } return value }