lapack/gonum: use Dcombssq (and Dlassq) in Dlange, Dlansy, Dlantr
@@ -22,6 +22,8 @@ func (impl Implementation) Dlantr(norm lapack.MatrixNorm, uplo blas.Uplo, diag b
panic(badUplo)
case diag != blas.Unit && diag != blas.NonUnit:
panic(badDiag)
+ case m < 0:
+ panic(mLT0)
case n < 0:
panic(nLT0)
case lda < max(1, n):
@@ -42,8 +44,6 @@ func (impl Implementation) Dlantr(norm lapack.MatrixNorm, uplo blas.Uplo, diag b
}
switch norm {
- default:
- panic(badNorm)
case lapack.MaxAbs:
if diag == blas.Unit {
value := 1.0
@@ -219,42 +219,38 @@ func (impl Implementation) Dlantr(norm lapack.MatrixNorm, uplo blas.Uplo, diag b
return maxsum
}
}
- case lapack.Frobenius:
- var nrm float64
+ default:
+ // lapack.Frobenius:
+ var scale, ssq float64
if diag == blas.Unit {
+ scale = 1
+ ssq = float64(min(m, n))
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
- }
+ for i := 0; i < min(m, n); i++ {
+ rowscale, rowssq := impl.Dlassq(n-i-1, a[i*lda+i+1:], 1, 0, 1)
+ scale, ssq = impl.Dcombssq(scale, ssq, rowscale, rowssq)
}
} else {
for i := 1; i < m; i++ {
- for j := 0; j < min(i, n); j++ {
- tmp := a[i*lda+j]
- nrm += tmp * tmp
- }
+ rowscale, rowssq := impl.Dlassq(min(i, n), a[i*lda:], 1, 0, 1)
+ scale, ssq = impl.Dcombssq(scale, ssq, rowscale, rowssq)
}
}
- nrm += float64(minmn)
} else {
+ scale = 0
+ ssq = 1
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
- }
+ for i := 0; i < min(m, n); i++ {
+ rowscale, rowssq := impl.Dlassq(n-i, a[i*lda+i:], 1, 0, 1)
+ scale, ssq = impl.Dcombssq(scale, ssq, rowscale, rowssq)
}
} 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
- }
+ rowscale, rowssq := impl.Dlassq(min(i+1, n), a[i*lda:], 1, 0, 1)
+ scale, ssq = impl.Dcombssq(scale, ssq, rowscale, rowssq)
}
}
}
- return math.Sqrt(nrm)
+ return scale * math.Sqrt(ssq)
}
}
