1 // Copyright ©2015 The Gonum Authors. All rights reserved.
2 // Use of this source code is governed by a BSD-style
3 // license that can be found in the LICENSE file.
8 "gonum.org/v1/gonum/blas"
9 "gonum.org/v1/gonum/blas/blas64"
12 // Dlabrd reduces the first NB rows and columns of a real general m×n matrix
13 // A to upper or lower bidiagonal form by an orthogonal transformation
15 // If m >= n, A is reduced to upper bidiagonal form and upon exit the elements
16 // on and below the diagonal in the first nb columns represent the elementary
17 // reflectors, and the elements above the diagonal in the first nb rows represent
18 // the matrix P. If m < n, A is reduced to lower bidiagonal form and the elements
19 // P is instead stored above the diagonal.
21 // The reduction to bidiagonal form is stored in d and e, where d are the diagonal
22 // elements, and e are the off-diagonal elements.
24 // The matrices Q and P are products of elementary reflectors
25 // Q = H_0 * H_1 * ... * H_{nb-1}
26 // P = G_0 * G_1 * ... * G_{nb-1}
28 // H_i = I - tauQ[i] * v_i * v_i^T
29 // G_i = I - tauP[i] * u_i * u_i^T
31 // As an example, on exit the entries of A when m = 6, n = 5, and nb = 2
38 // and when m = 5, n = 6, and nb = 2
39 // [ 1 u1 u1 u1 u1 u1]
45 // Dlabrd also returns the matrices X and Y which are used with U and V to
46 // apply the transformation to the unreduced part of the matrix
47 // A := A - V*Y^T - X*U^T
48 // and returns the matrices X and Y which are needed to apply the
49 // transformation to the unreduced part of A.
51 // X is an m×nb matrix, Y is an n×nb matrix. d, e, taup, and tauq must all have
52 // length at least nb. Dlabrd will panic if these size constraints are violated.
54 // Dlabrd is an internal routine. It is exported for testing purposes.
55 func (impl Implementation) Dlabrd(m, n, nb int, a []float64, lda int, d, e, tauQ, tauP, x []float64, ldx int, y []float64, ldy int) {
56 checkMatrix(m, n, a, lda)
57 checkMatrix(m, nb, x, ldx)
58 checkMatrix(n, nb, y, ldy)
74 bi := blas64.Implementation()
76 // Reduce to upper bidiagonal form.
77 for i := 0; i < nb; i++ {
78 bi.Dgemv(blas.NoTrans, m-i, i, -1, a[i*lda:], lda, y[i*ldy:], 1, 1, a[i*lda+i:], lda)
79 bi.Dgemv(blas.NoTrans, m-i, i, -1, x[i*ldx:], ldx, a[i:], lda, 1, a[i*lda+i:], lda)
81 a[i*lda+i], tauQ[i] = impl.Dlarfg(m-i, a[i*lda+i], a[min(i+1, m-1)*lda+i:], lda)
84 // Compute Y[i+1:n, i].
86 bi.Dgemv(blas.Trans, m-i, n-i-1, 1, a[i*lda+i+1:], lda, a[i*lda+i:], lda, 0, y[(i+1)*ldy+i:], ldy)
87 bi.Dgemv(blas.Trans, m-i, i, 1, a[i*lda:], lda, a[i*lda+i:], lda, 0, y[i:], ldy)
88 bi.Dgemv(blas.NoTrans, n-i-1, i, -1, y[(i+1)*ldy:], ldy, y[i:], ldy, 1, y[(i+1)*ldy+i:], ldy)
89 bi.Dgemv(blas.Trans, m-i, i, 1, x[i*ldx:], ldx, a[i*lda+i:], lda, 0, y[i:], ldy)
90 bi.Dgemv(blas.Trans, i, n-i-1, -1, a[i+1:], lda, y[i:], ldy, 1, y[(i+1)*ldy+i:], ldy)
91 bi.Dscal(n-i-1, tauQ[i], y[(i+1)*ldy+i:], ldy)
93 // Update A[i, i+1:n].
94 bi.Dgemv(blas.NoTrans, n-i-1, i+1, -1, y[(i+1)*ldy:], ldy, a[i*lda:], 1, 1, a[i*lda+i+1:], 1)
95 bi.Dgemv(blas.Trans, i, n-i-1, -1, a[i+1:], lda, x[i*ldx:], 1, 1, a[i*lda+i+1:], 1)
97 // Generate reflection P[i] to annihilate A[i, i+2:n].
98 a[i*lda+i+1], tauP[i] = impl.Dlarfg(n-i-1, a[i*lda+i+1], a[i*lda+min(i+2, n-1):], 1)
102 // Compute X[i+1:m, i].
103 bi.Dgemv(blas.NoTrans, m-i-1, n-i-1, 1, a[(i+1)*lda+i+1:], lda, a[i*lda+i+1:], 1, 0, x[(i+1)*ldx+i:], ldx)
104 bi.Dgemv(blas.Trans, n-i-1, i+1, 1, y[(i+1)*ldy:], ldy, a[i*lda+i+1:], 1, 0, x[i:], ldx)
105 bi.Dgemv(blas.NoTrans, m-i-1, i+1, -1, a[(i+1)*lda:], lda, x[i:], ldx, 1, x[(i+1)*ldx+i:], ldx)
106 bi.Dgemv(blas.NoTrans, i, n-i-1, 1, a[i+1:], lda, a[i*lda+i+1:], 1, 0, x[i:], ldx)
107 bi.Dgemv(blas.NoTrans, m-i-1, i, -1, x[(i+1)*ldx:], ldx, x[i:], ldx, 1, x[(i+1)*ldx+i:], ldx)
108 bi.Dscal(m-i-1, tauP[i], x[(i+1)*ldx+i:], ldx)
113 // Reduce to lower bidiagonal form.
114 for i := 0; i < nb; i++ {
116 bi.Dgemv(blas.NoTrans, n-i, i, -1, y[i*ldy:], ldy, a[i*lda:], 1, 1, a[i*lda+i:], 1)
117 bi.Dgemv(blas.Trans, i, n-i, -1, a[i:], lda, x[i*ldx:], 1, 1, a[i*lda+i:], 1)
119 // Generate reflection P[i] to annihilate A[i, i+1:n]
120 a[i*lda+i], tauP[i] = impl.Dlarfg(n-i, a[i*lda+i], a[i*lda+min(i+1, n-1):], 1)
124 // Compute X[i+1:m, i].
125 bi.Dgemv(blas.NoTrans, m-i-1, n-i, 1, a[(i+1)*lda+i:], lda, a[i*lda+i:], 1, 0, x[(i+1)*ldx+i:], ldx)
126 bi.Dgemv(blas.Trans, n-i, i, 1, y[i*ldy:], ldy, a[i*lda+i:], 1, 0, x[i:], ldx)
127 bi.Dgemv(blas.NoTrans, m-i-1, i, -1, a[(i+1)*lda:], lda, x[i:], ldx, 1, x[(i+1)*ldx+i:], ldx)
128 bi.Dgemv(blas.NoTrans, i, n-i, 1, a[i:], lda, a[i*lda+i:], 1, 0, x[i:], ldx)
129 bi.Dgemv(blas.NoTrans, m-i-1, i, -1, x[(i+1)*ldx:], ldx, x[i:], ldx, 1, x[(i+1)*ldx+i:], ldx)
130 bi.Dscal(m-i-1, tauP[i], x[(i+1)*ldx+i:], ldx)
132 // Update A[i+1:m, i].
133 bi.Dgemv(blas.NoTrans, m-i-1, i, -1, a[(i+1)*lda:], lda, y[i*ldy:], 1, 1, a[(i+1)*lda+i:], lda)
134 bi.Dgemv(blas.NoTrans, m-i-1, i+1, -1, x[(i+1)*ldx:], ldx, a[i:], lda, 1, a[(i+1)*lda+i:], lda)
136 // Generate reflection Q[i] to annihilate A[i+2:m, i].
137 a[(i+1)*lda+i], tauQ[i] = impl.Dlarfg(m-i-1, a[(i+1)*lda+i], a[min(i+2, m-1)*lda+i:], lda)
138 e[i] = a[(i+1)*lda+i]
141 // Compute Y[i+1:n, i].
142 bi.Dgemv(blas.Trans, m-i-1, n-i-1, 1, a[(i+1)*lda+i+1:], lda, a[(i+1)*lda+i:], lda, 0, y[(i+1)*ldy+i:], ldy)
143 bi.Dgemv(blas.Trans, m-i-1, i, 1, a[(i+1)*lda:], lda, a[(i+1)*lda+i:], lda, 0, y[i:], ldy)
144 bi.Dgemv(blas.NoTrans, n-i-1, i, -1, y[(i+1)*ldy:], ldy, y[i:], ldy, 1, y[(i+1)*ldy+i:], ldy)
145 bi.Dgemv(blas.Trans, m-i-1, i+1, 1, x[(i+1)*ldx:], ldx, a[(i+1)*lda+i:], lda, 0, y[i:], ldy)
146 bi.Dgemv(blas.Trans, i+1, n-i-1, -1, a[i+1:], lda, y[i:], ldy, 1, y[(i+1)*ldy+i:], ldy)
147 bi.Dscal(n-i-1, tauQ[i], y[(i+1)*ldy+i:], ldy)