+++ /dev/null
-// Copyright ©2015 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"
-)
-
-// Dlantr computes the specified norm of an m×n trapezoidal matrix A. If
-// norm == lapack.MaxColumnSum work must have length at least n, otherwise work
-// is unused.
-func (impl Implementation) Dlantr(norm lapack.MatrixNorm, uplo blas.Uplo, diag blas.Diag, m, n int, a []float64, lda int, work []float64) float64 {
- checkMatrix(m, n, a, lda)
- switch norm {
- case lapack.MaxRowSum, lapack.MaxColumnSum, lapack.NormFrob, lapack.MaxAbs:
- default:
- panic(badNorm)
- }
- if uplo != blas.Upper && uplo != blas.Lower {
- panic(badUplo)
- }
- if diag != blas.Unit && diag != blas.NonUnit {
- panic(badDiag)
- }
- if norm == lapack.MaxColumnSum && len(work) < n {
- panic(badWork)
- }
- if min(m, n) == 0 {
- return 0
- }
- switch norm {
- default:
- panic("unreachable")
- case lapack.MaxAbs:
- if diag == blas.Unit {
- value := 1.0
- if uplo == blas.Upper {
- for i := 0; i < m; i++ {
- for j := i + 1; j < n; j++ {
- tmp := math.Abs(a[i*lda+j])
- if math.IsNaN(tmp) {
- return tmp
- }
- if tmp > value {
- value = tmp
- }
- }
- }
- return value
- }
- for i := 1; i < m; i++ {
- for j := 0; j < min(i, n); j++ {
- tmp := math.Abs(a[i*lda+j])
- if math.IsNaN(tmp) {
- return tmp
- }
- if tmp > value {
- value = tmp
- }
- }
- }
- return value
- }
- var value float64
- if uplo == blas.Upper {
- for i := 0; i < m; i++ {
- for j := i; j < n; j++ {
- tmp := math.Abs(a[i*lda+j])
- if math.IsNaN(tmp) {
- return tmp
- }
- if tmp > value {
- value = tmp
- }
- }
- }
- return value
- }
- for i := 0; i < m; i++ {
- for j := 0; j <= min(i, n-1); j++ {
- tmp := math.Abs(a[i*lda+j])
- if math.IsNaN(tmp) {
- return tmp
- }
- if tmp > value {
- value = tmp
- }
- }
- }
- return value
- case lapack.MaxColumnSum:
- if diag == blas.Unit {
- for i := 0; i < min(m, n); i++ {
- work[i] = 1
- }
- for i := min(m, n); i < n; i++ {
- work[i] = 0
- }
- if uplo == blas.Upper {
- for i := 0; i < m; i++ {
- for j := i + 1; j < n; j++ {
- work[j] += math.Abs(a[i*lda+j])
- }
- }
- } else {
- for i := 1; i < m; i++ {
- for j := 0; j < min(i, n); j++ {
- work[j] += math.Abs(a[i*lda+j])
- }
- }
- }
- } else {
- for i := 0; i < n; i++ {
- work[i] = 0
- }
- if uplo == blas.Upper {
- for i := 0; i < m; i++ {
- for j := i; j < n; j++ {
- work[j] += math.Abs(a[i*lda+j])
- }
- }
- } else {
- for i := 0; i < m; i++ {
- for j := 0; j <= min(i, n-1); j++ {
- work[j] += math.Abs(a[i*lda+j])
- }
- }
- }
- }
- var max float64
- for _, v := range work[:n] {
- if math.IsNaN(v) {
- return math.NaN()
- }
- if v > max {
- max = v
- }
- }
- return max
- case lapack.MaxRowSum:
- var maxsum float64
- if diag == blas.Unit {
- if uplo == blas.Upper {
- for i := 0; i < m; i++ {
- var sum float64
- if i < min(m, n) {
- sum = 1
- }
- for j := i + 1; j < n; j++ {
- sum += math.Abs(a[i*lda+j])
- }
- if math.IsNaN(sum) {
- return math.NaN()
- }
- if sum > maxsum {
- maxsum = sum
- }
- }
- return maxsum
- } else {
- for i := 1; i < m; i++ {
- var sum float64
- if i < min(m, n) {
- sum = 1
- }
- for j := 0; j < min(i, n); j++ {
- sum += math.Abs(a[i*lda+j])
- }
- if math.IsNaN(sum) {
- return math.NaN()
- }
- if sum > maxsum {
- maxsum = sum
- }
- }
- return maxsum
- }
- } else {
- if uplo == blas.Upper {
- for i := 0; i < m; i++ {
- var sum float64
- for j := i; j < n; j++ {
- sum += math.Abs(a[i*lda+j])
- }
- if math.IsNaN(sum) {
- return sum
- }
- if sum > maxsum {
- maxsum = sum
- }
- }
- return maxsum
- } else {
- for i := 0; i < m; i++ {
- var sum float64
- for j := 0; j <= min(i, n-1); j++ {
- sum += math.Abs(a[i*lda+j])
- }
- if math.IsNaN(sum) {
- return sum
- }
- if sum > maxsum {
- maxsum = sum
- }
- }
- return maxsum
- }
- }
- case lapack.NormFrob:
- var nrm float64
- if diag == blas.Unit {
- if uplo == blas.Upper {
- for i := 0; i < m; i++ {
- for j := i + 1; j < n; j++ {
- tmp := a[i*lda+j]
- nrm += tmp * tmp
- }
- }
- } else {
- for i := 1; i < m; i++ {
- for j := 0; j < min(i, n); j++ {
- tmp := a[i*lda+j]
- nrm += tmp * tmp
- }
- }
- }
- nrm += float64(min(m, n))
- } else {
- if uplo == blas.Upper {
- for i := 0; i < m; i++ {
- for j := i; j < n; j++ {
- tmp := math.Abs(a[i*lda+j])
- nrm += tmp * tmp
- }
- }
- } else {
- for i := 0; i < m; i++ {
- for j := 0; j <= min(i, n-1); j++ {
- tmp := math.Abs(a[i*lda+j])
- nrm += tmp * tmp
- }
- }
- }
- }
- return math.Sqrt(nrm)
- }
-}
OSDN Git Service
