mirror of
https://github.com/gonum/gonum.git
synced 2025-10-21 06:09:26 +08:00
all: switch to subscript+brace notation for elementary reflectors
This commit is contained in:
@@ -254,13 +254,13 @@ func (impl Implementation) Dbdsqr(uplo blas.Uplo, n, ncvt, nru, ncc int, d, e, v
|
|||||||
//
|
//
|
||||||
// The remaining elements of A store the data needed to construct Q and P.
|
// The remaining elements of A store the data needed to construct Q and P.
|
||||||
// The matrices Q and P are products of elementary reflectors
|
// The matrices Q and P are products of elementary reflectors
|
||||||
// if m >= n, Q = H(0) * H(1) * ... * H(n-1),
|
// if m >= n, Q = H_0 * H_1 * ... * H_{n-1},
|
||||||
// P = G(0) * G(1) * ... * G(n-2),
|
// P = G_0 * G_1 * ... * G_{n-2},
|
||||||
// if m < n, Q = H(0) * H(1) * ... * H(m-2),
|
// if m < n, Q = H_0 * H_1 * ... * H_{m-2},
|
||||||
// P = G(0) * G(1) * ... * G(m-1),
|
// P = G_0 * G_1 * ... * G_{m-1},
|
||||||
// where
|
// where
|
||||||
// H(i) = I - tauQ[i] * v_i * v_i^T,
|
// H_i = I - tauQ[i] * v_i * v_i^T,
|
||||||
// G(i) = I - tauP[i] * u_i * u_i^T.
|
// G_i = I - tauP[i] * u_i * u_i^T.
|
||||||
//
|
//
|
||||||
// As an example, on exit the entries of A when m = 6, and n = 5
|
// As an example, on exit the entries of A when m = 6, and n = 5
|
||||||
// [ d e u1 u1 u1]
|
// [ d e u1 u1 u1]
|
||||||
@@ -354,7 +354,7 @@ func (impl Implementation) Dgecon(norm lapack.MatrixNorm, n int, a []float64, ld
|
|||||||
//
|
//
|
||||||
// See Dgeqr2 for a description of the elementary reflectors and orthonormal
|
// See Dgeqr2 for a description of the elementary reflectors and orthonormal
|
||||||
// matrix Q. Q is constructed as a product of these elementary reflectors,
|
// matrix Q. Q is constructed as a product of these elementary reflectors,
|
||||||
// Q = H(k-1) * ... * H(1) * H(0),
|
// Q = H_{k-1} * ... * H_1 * H_0,
|
||||||
// where k = min(m,n).
|
// where k = min(m,n).
|
||||||
//
|
//
|
||||||
// Work is temporary storage of length at least m and this function will panic otherwise.
|
// Work is temporary storage of length at least m and this function will panic otherwise.
|
||||||
@@ -412,10 +412,10 @@ func (impl Implementation) Dgelqf(m, n int, a []float64, lda int, tau, work []fl
|
|||||||
// v[j] = 0 j < i
|
// v[j] = 0 j < i
|
||||||
// v[j] = 1 j == i
|
// v[j] = 1 j == i
|
||||||
// v[j] = a[j*lda+i] j > i
|
// v[j] = a[j*lda+i] j > i
|
||||||
// and computing H(i) = I - tau[i] * v * v^T.
|
// and computing H_i = I - tau[i] * v * v^T.
|
||||||
//
|
//
|
||||||
// The orthonormal matrix Q can be constucted from a product of these elementary
|
// The orthonormal matrix Q can be constucted from a product of these elementary
|
||||||
// reflectors, Q = H(0) * H(1) * ... * H(k-1), where k = min(m,n).
|
// reflectors, Q = H_0 * H_1 * ... * H_{k-1}, where k = min(m,n).
|
||||||
//
|
//
|
||||||
// Work is temporary storage of length at least n and this function will panic otherwise.
|
// Work is temporary storage of length at least n and this function will panic otherwise.
|
||||||
func (impl Implementation) Dgeqr2(m, n int, a []float64, lda int, tau, work []float64) {
|
func (impl Implementation) Dgeqr2(m, n int, a []float64, lda int, tau, work []float64) {
|
||||||
@@ -739,7 +739,7 @@ func (impl Implementation) Dorgbr(vect lapack.DecompUpdate, m, n, k int, a []flo
|
|||||||
|
|
||||||
// Dorglq generates an m×n matrix Q with orthonormal rows defined by the product
|
// Dorglq generates an m×n matrix Q with orthonormal rows defined by the product
|
||||||
// of elementary reflectors
|
// of elementary reflectors
|
||||||
// Q = H(k-1) * ... * H(1) * H(0)
|
// Q = H_{k-1} * ... * H_1 * H_0
|
||||||
// as computed by Dgelqf. Dorglq is the blocked version of Dorgl2 that makes
|
// as computed by Dgelqf. Dorglq is the blocked version of Dorgl2 that makes
|
||||||
// greater use of level-3 BLAS routines.
|
// greater use of level-3 BLAS routines.
|
||||||
//
|
//
|
||||||
@@ -781,7 +781,7 @@ func (impl Implementation) Dorglq(m, n, k int, a []float64, lda int, tau, work [
|
|||||||
|
|
||||||
// Dorgqr generates an m×n matrix Q with orthonormal columns defined by the
|
// Dorgqr generates an m×n matrix Q with orthonormal columns defined by the
|
||||||
// product of elementary reflectors
|
// product of elementary reflectors
|
||||||
// Q = H(0) * H(1) * ... * H(k-1)
|
// Q = H_0 * H_1 * ... * H_{k-1}
|
||||||
// as computed by Dgeqrf. Dorgqr is the blocked version of Dorg2r that makes
|
// as computed by Dgeqrf. Dorgqr is the blocked version of Dorg2r that makes
|
||||||
// greater use of level-3 BLAS routines.
|
// greater use of level-3 BLAS routines.
|
||||||
//
|
//
|
||||||
|
@@ -48,8 +48,8 @@ type Float64 interface {
|
|||||||
type Direct byte
|
type Direct byte
|
||||||
|
|
||||||
const (
|
const (
|
||||||
Forward Direct = 'F' // Reflectors are right-multiplied, H(0) * H(1) * ... * H(k-1).
|
Forward Direct = 'F' // Reflectors are right-multiplied, H_0 * H_1 * ... * H_{k-1}.
|
||||||
Backward Direct = 'B' // Reflectors are left-multiplied, H(k-1) * ... * H(1) * H(0).
|
Backward Direct = 'B' // Reflectors are left-multiplied, H_{k-1} * ... * H_1 * H_0.
|
||||||
)
|
)
|
||||||
|
|
||||||
// Sort is the sorting order.
|
// Sort is the sorting order.
|
||||||
|
@@ -112,10 +112,10 @@ func Gels(trans blas.Transpose, a blas64.General, b blas64.General, work []float
|
|||||||
// v[j] = 0 j < i
|
// v[j] = 0 j < i
|
||||||
// v[j] = 1 j == i
|
// v[j] = 1 j == i
|
||||||
// v[j] = a[j*lda+i] j > i
|
// v[j] = a[j*lda+i] j > i
|
||||||
// and computing H(i) = I - tau[i] * v * v^T.
|
// and computing H_i = I - tau[i] * v * v^T.
|
||||||
//
|
//
|
||||||
// The orthonormal matrix Q can be constucted from a product of these elementary
|
// The orthonormal matrix Q can be constucted from a product of these elementary
|
||||||
// reflectors, Q = H(0) * H(1) * ... * H(k-1), where k = min(m,n).
|
// reflectors, Q = H_0 * H_1 * ... * H_{k-1}, where k = min(m,n).
|
||||||
//
|
//
|
||||||
// Work is temporary storage, and lwork specifies the usable memory length.
|
// Work is temporary storage, and lwork specifies the usable memory length.
|
||||||
// At minimum, lwork >= m and this function will panic otherwise.
|
// At minimum, lwork >= m and this function will panic otherwise.
|
||||||
@@ -135,7 +135,7 @@ func Geqrf(a blas64.General, tau, work []float64, lwork int) {
|
|||||||
//
|
//
|
||||||
// See Geqrf for a description of the elementary reflectors and orthonormal
|
// See Geqrf for a description of the elementary reflectors and orthonormal
|
||||||
// matrix Q. Q is constructed as a product of these elementary reflectors,
|
// matrix Q. Q is constructed as a product of these elementary reflectors,
|
||||||
// Q = H(k-1) * ... * H(1) * H(0).
|
// Q = H_{k-1} * ... * H_1 * H_0.
|
||||||
//
|
//
|
||||||
// Work is temporary storage, and lwork specifies the usable memory length.
|
// Work is temporary storage, and lwork specifies the usable memory length.
|
||||||
// At minimum, lwork >= m and this function will panic otherwise.
|
// At minimum, lwork >= m and this function will panic otherwise.
|
||||||
|
@@ -36,7 +36,7 @@ func (impl Implementation) Dgebd2(m, n int, a []float64, lda int, d, e, tauQ, ta
|
|||||||
a[i*lda+i], tauQ[i] = impl.Dlarfg(m-i, a[i*lda+i], a[min(i+1, m-1)*lda+i:], lda)
|
a[i*lda+i], tauQ[i] = impl.Dlarfg(m-i, a[i*lda+i], a[min(i+1, m-1)*lda+i:], lda)
|
||||||
d[i] = a[i*lda+i]
|
d[i] = a[i*lda+i]
|
||||||
a[i*lda+i] = 1
|
a[i*lda+i] = 1
|
||||||
// Apply H(i) to A[i:m, i+1:n] from the left.
|
// Apply H_i to A[i:m, i+1:n] from the left.
|
||||||
if i < n-1 {
|
if i < n-1 {
|
||||||
impl.Dlarf(blas.Left, m-i, n-i-1, a[i*lda+i:], lda, tauQ[i], a[i*lda+i+1:], lda, work)
|
impl.Dlarf(blas.Left, m-i, n-i-1, a[i*lda+i:], lda, tauQ[i], a[i*lda+i+1:], lda, work)
|
||||||
}
|
}
|
||||||
|
@@ -19,13 +19,13 @@ import (
|
|||||||
//
|
//
|
||||||
// The remaining elements of A store the data needed to construct Q and P.
|
// The remaining elements of A store the data needed to construct Q and P.
|
||||||
// The matrices Q and P are products of elementary reflectors
|
// The matrices Q and P are products of elementary reflectors
|
||||||
// if m >= n, Q = H(0) * H(1) * ... * H(n-1),
|
// if m >= n, Q = H_0 * H_1 * ... * H_{n-1},
|
||||||
// P = G(0) * G(1) * ... * G(n-2),
|
// P = G_0 * G_1 * ... * G_{n-2},
|
||||||
// if m < n, Q = H(0) * H(1) * ... * H(m-2),
|
// if m < n, Q = H_0 * H_1 * ... * H_{m-2},
|
||||||
// P = G(0) * G(1) * ... * G(m-1),
|
// P = G_0 * G_1 * ... * G_{m-1},
|
||||||
// where
|
// where
|
||||||
// H(i) = I - tauQ[i] * v_i * v_i^T,
|
// H_i = I - tauQ[i] * v_i * v_i^T,
|
||||||
// G(i) = I - tauP[i] * u_i * u_i^T.
|
// G_i = I - tauP[i] * u_i * u_i^T.
|
||||||
//
|
//
|
||||||
// As an example, on exit the entries of A when m = 6, and n = 5
|
// As an example, on exit the entries of A when m = 6, and n = 5
|
||||||
// [ d e u1 u1 u1]
|
// [ d e u1 u1 u1]
|
||||||
|
@@ -19,7 +19,7 @@ import "github.com/gonum/blas"
|
|||||||
//
|
//
|
||||||
// See Dgeqr2 for a description of the elementary reflectors and orthonormal
|
// See Dgeqr2 for a description of the elementary reflectors and orthonormal
|
||||||
// matrix Q. Q is constructed as a product of these elementary reflectors,
|
// matrix Q. Q is constructed as a product of these elementary reflectors,
|
||||||
// Q = H(k-1) * ... * H(1) * H(0).
|
// Q = H_{k-1} * ... * H_1 * H_0.
|
||||||
//
|
//
|
||||||
// work is temporary storage of length at least m and this function will panic otherwise.
|
// work is temporary storage of length at least m and this function will panic otherwise.
|
||||||
//
|
//
|
||||||
|
@@ -12,9 +12,9 @@ import "github.com/gonum/blas"
|
|||||||
// where Q is an m×m orthonormal matrix and L is a lower trapezoidal matrix.
|
// where Q is an m×m orthonormal matrix and L is a lower trapezoidal matrix.
|
||||||
//
|
//
|
||||||
// Q is represented as a product of elementary reflectors,
|
// Q is represented as a product of elementary reflectors,
|
||||||
// Q = H(k-1) * ... * H(1) * H(0)
|
// Q = H_{k-1} * ... * H_1 * H_0
|
||||||
// where k = min(m,n) and each H(i) has the form
|
// where k = min(m,n) and each H_i has the form
|
||||||
// H(i) = I - tau[i] * v_i * v_i^T
|
// H_i = I - tau[i] * v_i * v_i^T
|
||||||
// Vector v_i has v[m-k+i+1:m] = 0, v[m-k+i] = 1, and v[:m-k+i+1] is stored on
|
// Vector v_i has v[m-k+i+1:m] = 0, v[m-k+i] = 1, and v[:m-k+i+1] is stored on
|
||||||
// exit in A[0:m-k+i-1, n-k+i].
|
// exit in A[0:m-k+i-1, n-k+i].
|
||||||
//
|
//
|
||||||
@@ -34,10 +34,10 @@ func (impl Implementation) Dgeql2(m, n int, a []float64, lda int, tau, work []fl
|
|||||||
k := min(m, n)
|
k := min(m, n)
|
||||||
var aii float64
|
var aii float64
|
||||||
for i := k - 1; i >= 0; i-- {
|
for i := k - 1; i >= 0; i-- {
|
||||||
// Generate elementary reflector H(i) to annihilate A[0:m-k+i-1, n-k+i].
|
// Generate elementary reflector H_i to annihilate A[0:m-k+i-1, n-k+i].
|
||||||
aii, tau[i] = impl.Dlarfg(m-k+i+1, a[(m-k+i)*lda+n-k+i], a[n-k+i:], lda)
|
aii, tau[i] = impl.Dlarfg(m-k+i+1, a[(m-k+i)*lda+n-k+i], a[n-k+i:], lda)
|
||||||
|
|
||||||
// Apply H(i) to A[0:m-k+i, 0:n-k+i-1] from the left.
|
// Apply H_i to A[0:m-k+i, 0:n-k+i-1] from the left.
|
||||||
a[(m-k+i)*lda+n-k+i] = 1
|
a[(m-k+i)*lda+n-k+i] = 1
|
||||||
impl.Dlarf(blas.Left, m-k+i+1, n-k+i, a[n-k+i:], lda, tau[i], a, lda, work)
|
impl.Dlarf(blas.Left, m-k+i+1, n-k+i, a[n-k+i:], lda, tau[i], a, lda, work)
|
||||||
a[(m-k+i)*lda+n-k+i] = aii
|
a[(m-k+i)*lda+n-k+i] = aii
|
||||||
|
@@ -22,10 +22,10 @@ import "github.com/gonum/blas"
|
|||||||
// v[j] = 0 j < i
|
// v[j] = 0 j < i
|
||||||
// v[j] = 1 j == i
|
// v[j] = 1 j == i
|
||||||
// v[j] = a[j*lda+i] j > i
|
// v[j] = a[j*lda+i] j > i
|
||||||
// and computing H(i) = I - tau[i] * v * v^T.
|
// and computing H_i = I - tau[i] * v * v^T.
|
||||||
//
|
//
|
||||||
// The orthonormal matrix Q can be constructed from a product of these elementary
|
// The orthonormal matrix Q can be constructed from a product of these elementary
|
||||||
// reflectors, Q = H(0) * H(1) * ... * H(k-1), where k = min(m,n).
|
// reflectors, Q = H_0 * H_1 * ... * H_{k-1}, where k = min(m,n).
|
||||||
//
|
//
|
||||||
// work is temporary storage of length at least n and this function will panic otherwise.
|
// work is temporary storage of length at least n and this function will panic otherwise.
|
||||||
//
|
//
|
||||||
@@ -43,7 +43,7 @@ func (impl Implementation) Dgeqr2(m, n int, a []float64, lda int, tau, work []fl
|
|||||||
panic(badTau)
|
panic(badTau)
|
||||||
}
|
}
|
||||||
for i := 0; i < k; i++ {
|
for i := 0; i < k; i++ {
|
||||||
// Generate elementary reflector H(i).
|
// Generate elementary reflector H_i.
|
||||||
a[i*lda+i], tau[i] = impl.Dlarfg(m-i, a[i*lda+i], a[min((i+1), m-1)*lda+i:], lda)
|
a[i*lda+i], tau[i] = impl.Dlarfg(m-i, a[i*lda+i], a[min((i+1), m-1)*lda+i:], lda)
|
||||||
if i < n-1 {
|
if i < n-1 {
|
||||||
aii := a[i*lda+i]
|
aii := a[i*lda+i]
|
||||||
|
@@ -22,11 +22,11 @@ import (
|
|||||||
// elements, and e are the off-diagonal elements.
|
// elements, and e are the off-diagonal elements.
|
||||||
//
|
//
|
||||||
// The matrices Q and P are products of elementary reflectors
|
// The matrices Q and P are products of elementary reflectors
|
||||||
// Q = H(0) * H(1) * ... * H(nb-1)
|
// Q = H_0 * H_1 * ... * H_{nb-1}
|
||||||
// P = G(0) * G(1) * ... * G(nb-1)
|
// P = G_0 * G_1 * ... * G_{nb-1}
|
||||||
// where
|
// where
|
||||||
// H(i) = I - tauQ[i] * v_i * v_i^T
|
// H_i = I - tauQ[i] * v_i * v_i^T
|
||||||
// G(i) = I - tauP[i] * u_i * u_i^T
|
// G_i = I - tauP[i] * u_i * u_i^T
|
||||||
//
|
//
|
||||||
// As an example, on exit the entries of A when m = 6, n = 5, and nb = 2
|
// As an example, on exit the entries of A when m = 6, n = 5, and nb = 2
|
||||||
// [ 1 1 u1 u1 u1]
|
// [ 1 1 u1 u1 u1]
|
||||||
|
@@ -15,8 +15,8 @@ import (
|
|||||||
// H = I - V * T * V^T if store == lapack.ColumnWise
|
// H = I - V * T * V^T if store == lapack.ColumnWise
|
||||||
// H = I - V^T * T * V if store == lapack.RowWise
|
// H = I - V^T * T * V if store == lapack.RowWise
|
||||||
// H is defined by a product of the elementary reflectors where
|
// H is defined by a product of the elementary reflectors where
|
||||||
// H = H(0) * H(1) * ... * H(k-1) if direct == lapack.Forward
|
// H = H_0 * H_1 * ... * H_{k-1} if direct == lapack.Forward
|
||||||
// H = H(k-1) * ... * H(1) * H(0) if direct == lapack.Backward
|
// H = H_{k-1} * ... * H_1 * H_0 if direct == lapack.Backward
|
||||||
//
|
//
|
||||||
// t is a k×k triangular matrix. t is upper triangular if direct = lapack.Forward
|
// t is a k×k triangular matrix. t is upper triangular if direct = lapack.Forward
|
||||||
// and lower triangular otherwise. This function will panic if t is not of
|
// and lower triangular otherwise. This function will panic if t is not of
|
||||||
@@ -25,7 +25,7 @@ import (
|
|||||||
// store describes the storage of the elementary reflectors in v. Please see
|
// store describes the storage of the elementary reflectors in v. Please see
|
||||||
// Dlarfb for a description of layout.
|
// Dlarfb for a description of layout.
|
||||||
//
|
//
|
||||||
// tau contains the scalar factors of the elementary reflectors H(i).
|
// tau contains the scalar factors of the elementary reflectors H_i.
|
||||||
//
|
//
|
||||||
// Dlarft is an internal routine. It is exported for testing purposes.
|
// Dlarft is an internal routine. It is exported for testing purposes.
|
||||||
func (Implementation) Dlarft(direct lapack.Direct, store lapack.StoreV, n, k int,
|
func (Implementation) Dlarft(direct lapack.Direct, store lapack.StoreV, n, k int,
|
||||||
|
@@ -53,11 +53,11 @@ import (
|
|||||||
// H has the form
|
// H has the form
|
||||||
// I - tau * v * v^T
|
// I - tau * v * v^T
|
||||||
// If uplo == blas.Upper,
|
// If uplo == blas.Upper,
|
||||||
// Q = H(n-1) * H(n-2) * ... * H(n-nb)
|
// Q = H_{n-1} * H_{n-2} * ... * H_{n-nb}
|
||||||
// where v[:i-1] is stored in A[:i-1,i], v[i-1] = 1, and v[i:n] = 0.
|
// where v[:i-1] is stored in A[:i-1,i], v[i-1] = 1, and v[i:n] = 0.
|
||||||
//
|
//
|
||||||
// If uplo == blas.Lower,
|
// If uplo == blas.Lower,
|
||||||
// Q = H(0) * H(1) * ... * H(nb-1)
|
// Q = H_0 * H_1 * ... * H_{nb-1}
|
||||||
// where v[:i+1] = 0, v[i+1] = 1, and v[i+2:n] is stored in A[i+2:n,i].
|
// where v[:i+1] = 0, v[i+1] = 1, and v[i+2:n] is stored in A[i+2:n,i].
|
||||||
//
|
//
|
||||||
// The vectors v form the n×nb matrix V which is used with W to apply a
|
// The vectors v form the n×nb matrix V which is used with W to apply a
|
||||||
@@ -89,7 +89,7 @@ func (impl Implementation) Dlatrd(uplo blas.Uplo, n, nb int, a []float64, lda in
|
|||||||
a[i*lda+i+1:], 1, 1, a[i:], lda)
|
a[i*lda+i+1:], 1, 1, a[i:], lda)
|
||||||
}
|
}
|
||||||
if i > 0 {
|
if i > 0 {
|
||||||
// Generate elementary reflector H(i) to annihilate A(0:i-2,i).
|
// Generate elementary reflector H_i to annihilate A(0:i-2,i).
|
||||||
e[i-1], tau[i-1] = impl.Dlarfg(i, a[(i-1)*lda+i], a[i:], lda)
|
e[i-1], tau[i-1] = impl.Dlarfg(i, a[(i-1)*lda+i], a[i:], lda)
|
||||||
a[(i-1)*lda+i] = 1
|
a[(i-1)*lda+i] = 1
|
||||||
|
|
||||||
@@ -119,7 +119,7 @@ func (impl Implementation) Dlatrd(uplo blas.Uplo, n, nb int, a []float64, lda in
|
|||||||
bi.Dgemv(blas.NoTrans, n-i, i, -1, w[i*ldw:], ldw,
|
bi.Dgemv(blas.NoTrans, n-i, i, -1, w[i*ldw:], ldw,
|
||||||
a[i*lda:], 1, 1, a[i*lda+i:], lda)
|
a[i*lda:], 1, 1, a[i*lda+i:], lda)
|
||||||
if i < n-1 {
|
if i < n-1 {
|
||||||
// Generate elementary reflector H(i) to annihilate A(i+2:n,i).
|
// Generate elementary reflector H_i to annihilate A(i+2:n,i).
|
||||||
e[i], tau[i] = impl.Dlarfg(n-i-1, a[(i+1)*lda+i], a[min(i+2, n-1)*lda+i:], lda)
|
e[i], tau[i] = impl.Dlarfg(n-i-1, a[(i+1)*lda+i], a[min(i+2, n-1)*lda+i:], lda)
|
||||||
a[(i+1)*lda+i] = 1
|
a[(i+1)*lda+i] = 1
|
||||||
|
|
||||||
|
@@ -11,7 +11,7 @@ import (
|
|||||||
|
|
||||||
// Dorg2l generates an m×n matrix Q with orthonormal columns which is defined
|
// Dorg2l generates an m×n matrix Q with orthonormal columns which is defined
|
||||||
// as the last n columns of a product of k elementary reflectors of order m.
|
// as the last n columns of a product of k elementary reflectors of order m.
|
||||||
// Q = H(k-1) * ... * H(1) * H(0)
|
// Q = H_{k-1} * ... * H_1 * H_0
|
||||||
// See Dgelqf for more information. It must be that m >= n >= k.
|
// See Dgelqf for more information. It must be that m >= n >= k.
|
||||||
//
|
//
|
||||||
// tau contains the scalar reflectors computed by Dgeqlf. tau must have length
|
// tau contains the scalar reflectors computed by Dgeqlf. tau must have length
|
||||||
@@ -51,7 +51,7 @@ func (impl Implementation) Dorg2l(m, n, k int, a []float64, lda int, tau, work [
|
|||||||
for i := 0; i < k; i++ {
|
for i := 0; i < k; i++ {
|
||||||
ii := n - k + i
|
ii := n - k + i
|
||||||
|
|
||||||
// Apply H(i) to A[0:m-k+i, 0:n-k+i] from the left.
|
// Apply H_i to A[0:m-k+i, 0:n-k+i] from the left.
|
||||||
a[(m-n+ii)*lda+ii] = 1
|
a[(m-n+ii)*lda+ii] = 1
|
||||||
impl.Dlarf(blas.Left, m-n+ii+1, ii, a[ii:], lda, tau[i], a, lda, work)
|
impl.Dlarf(blas.Left, m-n+ii+1, ii, a[ii:], lda, tau[i], a, lda, work)
|
||||||
bi.Dscal(m-n+ii, -tau[i], a[ii:], lda)
|
bi.Dscal(m-n+ii, -tau[i], a[ii:], lda)
|
||||||
|
@@ -11,7 +11,7 @@ import (
|
|||||||
|
|
||||||
// Dorg2r generates an m×n matrix Q with orthonormal columns defined by the
|
// Dorg2r generates an m×n matrix Q with orthonormal columns defined by the
|
||||||
// product of elementary reflectors as computed by Dgeqrf.
|
// product of elementary reflectors as computed by Dgeqrf.
|
||||||
// Q = H(0) * H(2) * ... * H(k-1)
|
// Q = H_0 * H_1 * ... * H_{k-1}
|
||||||
// len(tau) >= k, 0 <= k <= n, 0 <= n <= m, len(work) >= n.
|
// len(tau) >= k, 0 <= k <= n, 0 <= n <= m, len(work) >= n.
|
||||||
// Dorg2r will panic if these conditions are not met.
|
// Dorg2r will panic if these conditions are not met.
|
||||||
//
|
//
|
||||||
|
@@ -11,7 +11,7 @@ import (
|
|||||||
|
|
||||||
// Dorgl2 generates an m×n matrix Q with orthonormal rows defined by the
|
// Dorgl2 generates an m×n matrix Q with orthonormal rows defined by the
|
||||||
// first m rows product of elementary reflectors as computed by Dgelqf.
|
// first m rows product of elementary reflectors as computed by Dgelqf.
|
||||||
// Q = H(0) * H(2) * ... * H(k-1)
|
// Q = H_0 * H_1 * ... * H_{k-1}
|
||||||
// len(tau) >= k, 0 <= k <= m, 0 <= m <= n, len(work) >= m.
|
// len(tau) >= k, 0 <= k <= m, 0 <= m <= n, len(work) >= m.
|
||||||
// Dorgl2 will panic if these conditions are not met.
|
// Dorgl2 will panic if these conditions are not met.
|
||||||
//
|
//
|
||||||
|
@@ -11,7 +11,7 @@ import (
|
|||||||
|
|
||||||
// Dorglq generates an m×n matrix Q with orthonormal columns defined by the
|
// Dorglq generates an m×n matrix Q with orthonormal columns defined by the
|
||||||
// product of elementary reflectors as computed by Dgelqf.
|
// product of elementary reflectors as computed by Dgelqf.
|
||||||
// Q = H(0) * H(2) * ... * H(k-1)
|
// Q = H_0 * H_1 * ... * H_{k-1}
|
||||||
// Dorglq is the blocked version of Dorgl2 that makes greater use of level-3 BLAS
|
// Dorglq is the blocked version of Dorgl2 that makes greater use of level-3 BLAS
|
||||||
// routines.
|
// routines.
|
||||||
//
|
//
|
||||||
|
@@ -10,10 +10,9 @@ import (
|
|||||||
)
|
)
|
||||||
|
|
||||||
// Dorgql generates the m×n matrix Q with orthonormal columns defined as the
|
// Dorgql generates the m×n matrix Q with orthonormal columns defined as the
|
||||||
// last n columns of a product of k elementary reflectors of order m as returned
|
// last n columns of a product of k elementary reflectors of order m
|
||||||
// by Dgelqf. That is,
|
// Q = H_{k-1} * ... * H_1 * H_0
|
||||||
// Q = H(k-1) * ... * H(1) * H(0).
|
// as returned by Dgelqf. See Dgelqf for more information.
|
||||||
// See Dgelqf for more information.
|
|
||||||
//
|
//
|
||||||
// tau must have length at least k, and Dorgql will panic otherwise.
|
// tau must have length at least k, and Dorgql will panic otherwise.
|
||||||
//
|
//
|
||||||
@@ -83,7 +82,7 @@ func (impl Implementation) Dorgql(m, n, k int, a []float64, lda int, tau, work [
|
|||||||
ib := min(nb, k-i)
|
ib := min(nb, k-i)
|
||||||
if n-k+i > 0 {
|
if n-k+i > 0 {
|
||||||
// Form the triangular factor of the block reflector
|
// Form the triangular factor of the block reflector
|
||||||
// H = H(i+ib-1) * ... * H(i+1) * H(i).
|
// H = H_{i+ib-1} * ... * H_{i+1} * H_i.
|
||||||
impl.Dlarft(lapack.Backward, lapack.ColumnWise, m-k+i+ib, ib,
|
impl.Dlarft(lapack.Backward, lapack.ColumnWise, m-k+i+ib, ib,
|
||||||
a[n-k+i:], lda, tau[i:], work, ldwork)
|
a[n-k+i:], lda, tau[i:], work, ldwork)
|
||||||
|
|
||||||
|
@@ -10,8 +10,9 @@ import (
|
|||||||
)
|
)
|
||||||
|
|
||||||
// Dorgqr generates an m×n matrix Q with orthonormal columns defined by the
|
// Dorgqr generates an m×n matrix Q with orthonormal columns defined by the
|
||||||
// product of elementary reflectors as computed by Dgeqrf.
|
// product of elementary reflectors
|
||||||
// Q = H(0) * H(2) * ... * H(k-1)
|
// Q = H_0 * H_1 * ... * H_{k-1}
|
||||||
|
// as computed by Dgeqrf.
|
||||||
// Dorgqr is the blocked version of Dorg2r that makes greater use of level-3 BLAS
|
// Dorgqr is the blocked version of Dorg2r that makes greater use of level-3 BLAS
|
||||||
// routines.
|
// routines.
|
||||||
//
|
//
|
||||||
|
@@ -10,9 +10,9 @@ import "github.com/gonum/blas"
|
|||||||
// of n-1 elementary reflectors of order n as returned by Dsytrd.
|
// of n-1 elementary reflectors of order n as returned by Dsytrd.
|
||||||
//
|
//
|
||||||
// The construction of Q depends on the value of uplo:
|
// The construction of Q depends on the value of uplo:
|
||||||
// Q = H(n-1) * ... * H(1) * H(0) if uplo == blas.Upper
|
// Q = H_{n-1} * ... * H_1 * H_0 if uplo == blas.Upper
|
||||||
// Q = H(0) * H(1) * ... * H(n-1) if uplo == blas.Lower
|
// Q = H_0 * H_1 * ... * H_{n-1} if uplo == blas.Lower
|
||||||
// where H(i) is constructed from the elementary reflectors as computed by Dsytrd.
|
// where H_i is constructed from the elementary reflectors as computed by Dsytrd.
|
||||||
// See the documentation for Dsytrd for more information.
|
// See the documentation for Dsytrd for more information.
|
||||||
//
|
//
|
||||||
// tau must have length at least n-1, and Dorgtr will panic otherwise.
|
// tau must have length at least n-1, and Dorgtr will panic otherwise.
|
||||||
|
@@ -24,11 +24,11 @@ import (
|
|||||||
//
|
//
|
||||||
// Q is represented as a product of elementary reflectors.
|
// Q is represented as a product of elementary reflectors.
|
||||||
// If uplo == blas.Upper
|
// If uplo == blas.Upper
|
||||||
// Q = H(n-2) * ... * H(1) * H(0)
|
// Q = H_{n-2} * ... * H_1 * H_0
|
||||||
// and if uplo == blas.Lower
|
// and if uplo == blas.Lower
|
||||||
// Q = H(0) * H(1) * ... * H(n-2)
|
// Q = H_0 * H_1 * ... * H_{n-2}
|
||||||
// where
|
// where
|
||||||
// H(i) = I - tau * v * v^T
|
// H_i = I - tau * v * v^T
|
||||||
// where tau is stored in tau[i], and v is stored in a.
|
// where tau is stored in tau[i], and v is stored in a.
|
||||||
//
|
//
|
||||||
// If uplo == blas.Upper, v[0:i-1] is stored in A[0:i-1,i+1], v[i] = 1, and
|
// If uplo == blas.Upper, v[0:i-1] is stored in A[0:i-1,i+1], v[i] = 1, and
|
||||||
@@ -66,13 +66,13 @@ func (impl Implementation) Dsytd2(uplo blas.Uplo, n int, a []float64, lda int, d
|
|||||||
if uplo == blas.Upper {
|
if uplo == blas.Upper {
|
||||||
// Reduce the upper triangle of A.
|
// Reduce the upper triangle of A.
|
||||||
for i := n - 2; i >= 0; i-- {
|
for i := n - 2; i >= 0; i-- {
|
||||||
// Generate elementary reflector H(i) = I - tau * v * v^T to
|
// Generate elementary reflector H_i = I - tau * v * v^T to
|
||||||
// annihilate A[i:i-1, i+1].
|
// annihilate A[i:i-1, i+1].
|
||||||
var taui float64
|
var taui float64
|
||||||
a[i*lda+i+1], taui = impl.Dlarfg(i+1, a[i*lda+i+1], a[i+1:], lda)
|
a[i*lda+i+1], taui = impl.Dlarfg(i+1, a[i*lda+i+1], a[i+1:], lda)
|
||||||
e[i] = a[i*lda+i+1]
|
e[i] = a[i*lda+i+1]
|
||||||
if taui != 0 {
|
if taui != 0 {
|
||||||
// Apply H(i) from both sides to A[0:i,0:i].
|
// Apply H_i from both sides to A[0:i,0:i].
|
||||||
a[i*lda+i+1] = 1
|
a[i*lda+i+1] = 1
|
||||||
|
|
||||||
// Compute x := tau * A * v storing x in tau[0:i].
|
// Compute x := tau * A * v storing x in tau[0:i].
|
||||||
@@ -95,13 +95,13 @@ func (impl Implementation) Dsytd2(uplo blas.Uplo, n int, a []float64, lda int, d
|
|||||||
}
|
}
|
||||||
// Reduce the lower triangle of A.
|
// Reduce the lower triangle of A.
|
||||||
for i := 0; i < n-1; i++ {
|
for i := 0; i < n-1; i++ {
|
||||||
// Generate elementary reflector H(i) = I - tau * v * v^T to
|
// Generate elementary reflector H_i = I - tau * v * v^T to
|
||||||
// annihilate A[i+2:n, i].
|
// annihilate A[i+2:n, i].
|
||||||
var taui float64
|
var taui float64
|
||||||
a[(i+1)*lda+i], taui = impl.Dlarfg(n-i-1, a[(i+1)*lda+i], a[min(i+2, n-1)*lda+i:], lda)
|
a[(i+1)*lda+i], taui = impl.Dlarfg(n-i-1, a[(i+1)*lda+i], a[min(i+2, n-1)*lda+i:], lda)
|
||||||
e[i] = a[(i+1)*lda+i]
|
e[i] = a[(i+1)*lda+i]
|
||||||
if taui != 0 {
|
if taui != 0 {
|
||||||
// Apply H(i) from boh sides to A[i+1:n, i+1:n].
|
// Apply H_i from both sides to A[i+1:n, i+1:n].
|
||||||
a[(i+1)*lda+i] = 1
|
a[(i+1)*lda+i] = 1
|
||||||
|
|
||||||
// Compute x := tau * A * v, storing y in tau[i:n-1].
|
// Compute x := tau * A * v, storing y in tau[i:n-1].
|
||||||
|
@@ -21,9 +21,9 @@ import (
|
|||||||
// the product of elementary reflectors.
|
// the product of elementary reflectors.
|
||||||
//
|
//
|
||||||
// If uplo == blas.Upper, Q is constructed with
|
// If uplo == blas.Upper, Q is constructed with
|
||||||
// Q = H(n-2) * ... * H(1) * H(0)
|
// Q = H_{n-2} * ... * H_1 * H_0
|
||||||
// where
|
// where
|
||||||
// H(i) = I - tau * v * v^T
|
// H_i = I - tau * v * v^T
|
||||||
// v is constructed as v[i+1:n] = 0, v[i] = 1, v[0:i-1] is stored in A[0:i-1, i+1],
|
// v is constructed as v[i+1:n] = 0, v[i] = 1, v[0:i-1] is stored in A[0:i-1, i+1],
|
||||||
// and tau is in tau[i]. The elements of A are
|
// and tau is in tau[i]. The elements of A are
|
||||||
// [ d e v2 v3 v4]
|
// [ d e v2 v3 v4]
|
||||||
@@ -33,11 +33,11 @@ import (
|
|||||||
// [ e]
|
// [ e]
|
||||||
//
|
//
|
||||||
// If uplo == blas.Lower, Q is constructed with
|
// If uplo == blas.Lower, Q is constructed with
|
||||||
// Q = H(0) * H(1) * ... * H(n-2)
|
// Q = H_0 * H_1 * ... * H_{n-2}
|
||||||
// where
|
// where
|
||||||
// H(i) = I - tau * v * v^T
|
// H_i = I - tau[i] * v * v^T
|
||||||
// v is constructed as v[0:i+1] = 0, v[i+1] = 1, v[i+2:n] is stored in A[i+2:n, i],
|
// v is constructed as v[0:i+1] = 0, v[i+1] = 1, v[i+2:n] is stored in A[i+2:n, i].
|
||||||
// and tau is in tau[i]. The elements of A are
|
// The elements of A are
|
||||||
// [ d ]
|
// [ d ]
|
||||||
// [ e d ]
|
// [ e d ]
|
||||||
// [v1 e d ]
|
// [v1 e d ]
|
||||||
|
Reference in New Issue
Block a user