lapack/gonum: use blocking crossover point in Dgebrd
MODIFY
lapack/gonum/dgebrd.go
lapack/gonum/dgebrd.go
@@ -53,46 +53,62 @@ import (
//
// Dgebrd is an internal routine. It is exported for testing purposes.
func (impl Implementation) Dgebrd(m, n int, a []float64, lda int, d, e, tauQ, tauP, work []float64, lwork int) {
- checkMatrix(m, n, a, lda)
- // Calculate optimal work.
- nb := impl.Ilaenv(1, "DGEBRD", " ", m, n, -1, -1)
- var lworkOpt int
- if lwork == -1 {
- if len(work) < 1 {
+ switch {
+ case m < 0:
+ panic(mLT0)
+ case n < 0:
+ panic(nLT0)
+ case lda < max(1, n):
+ panic(badLdA)
+ case lwork < max(1, max(m, n)) && lwork != -1:
panic(badWork)
+ case len(work) < max(1, lwork):
+ panic(shortWork)
}
- lworkOpt = ((m + n) * nb)
- work[0] = float64(max(1, lworkOpt))
+
+ // Quick return if possible.
+ minmn := min(m, n)
+ if minmn == 0 {
+ work[0] = 1
return
}
- minmn := min(m, n)
- if len(d) < minmn {
- panic(badD)
+
+ nb := impl.Ilaenv(1, "DGEBRD", " ", m, n, -1, -1)
+ lwkopt := (m + n) * nb
+ if lwork == -1 {
+ work[0] = float64(lwkopt)
+ return
}
- if len(e) < minmn-1 {
+
+ switch {
+ case len(a) < (m-1)*lda+n:
+ panic("lapack: insufficient length of a")
+ case len(d) < minmn:
+ panic(badD)
+ case len(e) < minmn-1:
panic(badE)
- }
- if len(tauQ) < minmn {
+ case len(tauQ) < minmn:
panic(badTauQ)
- }
- if len(tauP) < minmn {
+ case len(tauP) < minmn:
panic(badTauP)
}
+
+ nx := minmn
ws := max(m, n)
- if lwork < max(1, ws) {
- panic(badWork)
- }
- if len(work) < lwork {
- panic(badWork)
- }
- var nx int
- if nb > 1 && nb < minmn {
+ if 1 < nb && nb < minmn {
+ // At least one blocked operation can be done.
+ // Get the crossover point nx.
nx = max(nb, impl.Ilaenv(3, "DGEBRD", " ", m, n, -1, -1))
+ // Determine when to switch from blocked to unblocked code.
if nx < minmn {
+ // At least one blocked operation will be done.
ws = (m + n) * nb
if lwork < ws {
+ // Not enough work space for the optimal nb,
+ // consider using a smaller block size.
nbmin := impl.Ilaenv(2, "DGEBRD", " ", m, n, -1, -1)
if lwork >= (m+n)*nbmin {
+ // Enough work space for minimum block size.
nb = lwork / (m + n)
} else {
nb = minmn
@@ -100,17 +116,12 @@ func (impl Implementation) Dgebrd(m, n int, a []float64, lda int, d, e, tauQ, ta
}
}
}
- } else {
- nx = minmn
}
bi := blas64.Implementation()
ldworkx := nb
ldworky := nb
var i int
- // Netlib lapack has minmn - nx, but this makes the last nx rows (which by
- // default is large) be unblocked. As written here, the blocking is more
- // consistent.
- for i = 0; i < minmn-nb; i += nb {
+ for i = 0; i < minmn-nx; i += nb {
// Reduce rows and columns i:i+nb to bidiagonal form and return
// the matrices X and Y which are needed to update the unreduced
// part of the matrix.
@@ -146,5 +157,5 @@ func (impl Implementation) Dgebrd(m, n int, a []float64, lda int, d, e, tauQ, ta
}
// Use unblocked code to reduce the remainder of the matrix.
impl.Dgebd2(m-i, n-i, a[i*lda+i:], lda, d[i:], e[i:], tauQ[i:], tauP[i:], work)
- work[0] = float64(lworkOpt)
+ work[0] = float64(ws)
}

lapack/gonum: use blocking crossover point in Dgebrd