1 // Copyright ©2017 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.
10 "gonum.org/v1/gonum/blas"
11 "gonum.org/v1/gonum/blas/blas64"
14 // Dlaqps computes a step of QR factorization with column pivoting
15 // of an m×n matrix A by using Blas-3. It tries to factorize nb
16 // columns from A starting from the row offset, and updates all
17 // of the matrix with Dgemm.
19 // In some cases, due to catastrophic cancellations, it cannot
20 // factorize nb columns. Hence, the actual number of factorized
21 // columns is returned in kb.
23 // Dlaqps computes a QR factorization with column pivoting of the
24 // block A[offset:m, 0:nb] of the m×n matrix A. The block
25 // A[0:offset, 0:n] is accordingly pivoted, but not factorized.
27 // On exit, the upper triangle of block A[offset:m, 0:kb] is the
28 // triangular factor obtained. The elements in block A[offset:m, 0:n]
29 // below the diagonal, together with tau, represent the orthogonal
30 // matrix Q as a product of elementary reflectors.
32 // offset is number of rows of the matrix A that must be pivoted but
33 // not factorized. offset must not be negative otherwise Dlaqps will panic.
35 // On exit, jpvt holds the permutation that was applied; the jth column
36 // of A*P was the jpvt[j] column of A. jpvt must have length n,
37 // otherwise Dlapqs will panic.
39 // On exit tau holds the scalar factors of the elementary reflectors.
40 // It must have length nb, otherwise Dlapqs will panic.
42 // vn1 and vn2 hold the partial and complete column norms respectively.
43 // They must have length n, otherwise Dlapqs will panic.
45 // auxv must have length nb, otherwise Dlaqps will panic.
47 // f and ldf represent an n×nb matrix F that is overwritten during the
50 // Dlaqps is an internal routine. It is exported for testing purposes.
51 func (impl Implementation) Dlaqps(m, n, offset, nb int, a []float64, lda int, jpvt []int, tau, vn1, vn2, auxv, f []float64, ldf int) (kb int) {
52 checkMatrix(m, n, a, lda)
53 checkMatrix(n, nb, f, ldf)
76 lastrk := min(m, n+offset)
78 tol3z := math.Sqrt(dlamchE)
80 bi := blas64.Implementation()
83 for ; k < nb && lsticc == -1; k++ {
86 // Determine kth pivot column and swap if necessary.
87 p := k + bi.Idamax(n-k, vn1[k:], 1)
89 bi.Dswap(m, a[p:], lda, a[k:], lda)
90 bi.Dswap(k, f[p*ldf:], 1, f[k*ldf:], 1)
91 jpvt[p], jpvt[k] = jpvt[k], jpvt[p]
96 // Apply previous Householder reflectors to column K:
98 // A[rk:m, k] = A[rk:m, k] - A[rk:m, 0:k-1]*F[k, 0:k-1]^T.
100 bi.Dgemv(blas.NoTrans, m-rk, k, -1,
107 // Generate elementary reflector H_k.
109 a[rk*lda+k], tau[k] = impl.Dlarfg(m-rk, a[rk*lda+k], a[(rk+1)*lda+k:], lda)
117 // Compute kth column of F:
119 // Compute F[k+1:n, k] = tau[k]*A[rk:m, k+1:n]^T*A[rk:m, k].
121 bi.Dgemv(blas.Trans, m-rk, n-k-1, tau[k],
125 f[(k+1)*ldf+k:], ldf)
128 // Padding F[0:k, k] with zeros.
129 for j := 0; j < k; j++ {
133 // Incremental updating of F:
135 // F[0:n, k] := F[0:n, k] - tau[k]*F[0:n, 0:k-1]*A[rk:m, 0:k-1]^T*A[rk:m,k].
137 bi.Dgemv(blas.Trans, m-rk, k, -tau[k],
142 bi.Dgemv(blas.NoTrans, n, k, 1,
149 // Update the current row of A:
151 // A[rk, k+1:n] = A[rk, k+1:n] - A[rk, 0:k]*F[k+1:n, 0:k]^T.
153 bi.Dgemv(blas.NoTrans, n-k-1, k+1, -1,
160 // Update partial column norms.
162 for j := k + 1; j < n; j++ {
167 // The following marked lines follow from the
168 // analysis in Lapack Working Note 176.
169 r := math.Abs(a[rk*lda+j]) / vn1[j] // *
170 temp := math.Max(0, 1-r*r) // *
171 r = vn1[j] / vn2[j] // *
172 temp2 := temp * r * r // *
174 // vn2 is used here as a collection of
175 // indices into vn2 and also a collection
177 vn2[j] = float64(lsticc)
180 vn1[j] *= math.Sqrt(temp) // *
190 // Apply the block reflector to the rest of the matrix:
192 // A[offset+kb+1:m, kb+1:n] := A[offset+kb+1:m, kb+1:n] - A[offset+kb+1:m, 1:kb]*F[kb+1:n, 1:kb]^T.
193 if kb < min(n, m-offset) {
194 bi.Dgemm(blas.NoTrans, blas.Trans,
202 // Recomputation of difficult columns.
204 itemp := int(vn2[lsticc])
206 // NOTE: The computation of vn1[lsticc] relies on the fact that
207 // Dnrm2 does not fail on vectors with norm below the value of
209 v := bi.Dnrm2(m-rk, a[rk*lda+lsticc:], lda)