This repository has been archived by the owner on Nov 24, 2018. It is now read-only.
-
Notifications
You must be signed in to change notification settings - Fork 11
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Add dorgqr and tests, update dorg2r tests, and fix cgo dorgqr
- Loading branch information
Showing
12 changed files
with
520 additions
and
36 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,111 @@ | ||
package native | ||
|
||
import ( | ||
"github.com/gonum/blas" | ||
"github.com/gonum/lapack" | ||
) | ||
|
||
// Dorglq generates an m×n matrix Q with orthonormal columns defined by the | ||
// product of elementary reflectors as computed by Dgelqf. | ||
// Q = H(0) * H(2) * ... * H(k-1) | ||
// Dorglq is the blocked version of dorgl2 that makes greater use of level-3 BLAS | ||
// routines. | ||
// | ||
// len(tau) >= k, 0 <= k <= n, and 0 <= n <= m. | ||
// | ||
// Work is temporary storage, and lwork specifies the usable memory length. At minimum, | ||
// lwork >= m, and the amount of blocking is limited by the usable length. | ||
// If lwork == -1, instead of computing Dorglq the optimal work length is stored | ||
// into work[0]. | ||
// | ||
// Dorglq will panic if the conditions on input values are not met. | ||
func (impl Implementation) Dorglq(m, n, k int, a []float64, lda int, tau, work []float64, lwork int) { | ||
nb := impl.Ilaenv(1, "DORGLQ", " ", m, n, k, -1) | ||
// work is treated as an n×nb matrix | ||
if lwork == -1 { | ||
work[0] = float64(max(1, m) * nb) | ||
return | ||
} | ||
checkMatrix(m, n, a, lda) | ||
if k < 0 { | ||
panic(kLT0) | ||
} | ||
if k > m { | ||
panic(kGTM) | ||
} | ||
if m > n { | ||
panic(nLTM) | ||
} | ||
if len(tau) < k { | ||
panic(badTau) | ||
} | ||
if len(work) < lwork { | ||
panic(shortWork) | ||
} | ||
if lwork < m { | ||
panic(badWork) | ||
} | ||
if m == 0 { | ||
return | ||
} | ||
nbmin := 2 // Minimum number of blocks | ||
var nx int // Minimum number of rows | ||
iws := m // Length of work needed | ||
var ldwork int | ||
if nb > 1 && nb < k { | ||
nx = max(0, impl.Ilaenv(3, "DORGLQ", " ", m, n, k, -1)) | ||
if nx < k { | ||
ldwork = nb | ||
iws = m * ldwork | ||
if lwork < iws { | ||
nb = lwork / m | ||
ldwork = nb | ||
nbmin = max(2, impl.Ilaenv(2, "DORGLQ", " ", m, n, k, -1)) | ||
} | ||
} | ||
} | ||
var ki, kk int | ||
if nb >= nbmin && nb < k && nx < k { | ||
// The first kk rows are handled by the blocked method. | ||
// Note: lapack has nx here, but this means the last nx rows are handled | ||
// serially which could be quite different than nb. | ||
ki = ((k - nb - 1) / nb) * nb | ||
kk = min(k, ki+nb) | ||
for i := kk; i < m; i++ { | ||
for j := 0; j < kk; j++ { | ||
a[i*lda+j] = 0 | ||
} | ||
} | ||
} | ||
if kk < m { | ||
// Perform the operation on colums kk to the end. | ||
impl.Dorgl2(m-kk, n-kk, k-kk, a[kk*lda+kk:], lda, tau[kk:], work) | ||
} | ||
if kk == 0 { | ||
return | ||
} | ||
// Perform the operation on column-blocks | ||
for i := ki; i >= 0; i -= nb { | ||
ib := min(nb, k-i) | ||
if i+ib < m { | ||
impl.Dlarft(lapack.Forward, lapack.RowWise, | ||
n-i, ib, | ||
a[i*lda+i:], lda, | ||
tau[i:], | ||
work, ldwork) | ||
|
||
impl.Dlarfb(blas.Right, blas.Trans, lapack.Forward, lapack.RowWise, | ||
m-i-ib, n-i, ib, | ||
a[i*lda+i:], lda, | ||
work, ldwork, | ||
a[(i+ib)*lda+i:], lda, | ||
work[ib*ldwork:], ldwork) | ||
} | ||
impl.Dorgl2(ib, n-i, ib, a[i*lda+i:], lda, tau[i:], work) | ||
for l := i; l < i+ib; l++ { | ||
for j := 0; j < i; j++ { | ||
a[l*lda+j] = 0 | ||
} | ||
} | ||
} | ||
} |
Oops, something went wrong.