blob: cc3bc06db66d135002d82ff7a0146e031bcee15b [file] [log] [blame]
 // Copyright ©2017 The Gonum Authors. All rights reserved. // Use of this source code is governed by a BSD-style // license that can be found in the LICENSE file. package gonum import ( "math" "gonum.org/v1/gonum/blas" "gonum.org/v1/gonum/blas/blas64" ) // Dlaqp2 computes a QR factorization with column pivoting of the block A[offset:m, 0:n] // of the m×n matrix A. The block A[0:offset, 0:n] is accordingly pivoted, but not factorized. // // On exit, the upper triangle of block A[offset:m, 0:n] is the triangular factor obtained. // The elements in block A[offset:m, 0:n] below the diagonal, together with tau, represent // the orthogonal matrix Q as a product of elementary reflectors. // // offset is number of rows of the matrix A that must be pivoted but not factorized. // offset must not be negative otherwise Dlaqp2 will panic. // // On exit, jpvt holds the permutation that was applied; the jth column of A*P was the // jpvt[j] column of A. jpvt must have length n, otherwise Dlaqp2 will panic. // // On exit tau holds the scalar factors of the elementary reflectors. It must have length // at least min(m-offset, n) otherwise Dlaqp2 will panic. // // vn1 and vn2 hold the partial and complete column norms respectively. They must have length n, // otherwise Dlaqp2 will panic. // // work must have length n, otherwise Dlaqp2 will panic. // // Dlaqp2 is an internal routine. It is exported for testing purposes. func (impl Implementation) Dlaqp2(m, n, offset int, a []float64, lda int, jpvt []int, tau, vn1, vn2, work []float64) { switch { case m < 0: panic(mLT0) case n < 0: panic(nLT0) case offset < 0: panic(offsetLT0) case offset > m: panic(offsetGTM) case lda < max(1, n): panic(badLdA) } // Quick return if possible. if m == 0 || n == 0 { return } mn := min(m-offset, n) switch { case len(a) < (m-1)*lda+n: panic(shortA) case len(jpvt) != n: panic(badLenJpvt) case len(tau) < mn: panic(shortTau) case len(vn1) < n: panic(shortVn1) case len(vn2) < n: panic(shortVn2) case len(work) < n: panic(shortWork) } tol3z := math.Sqrt(dlamchE) bi := blas64.Implementation() // Compute factorization. for i := 0; i < mn; i++ { offpi := offset + i // Determine ith pivot column and swap if necessary. p := i + bi.Idamax(n-i, vn1[i:], 1) if p != i { bi.Dswap(m, a[p:], lda, a[i:], lda) jpvt[p], jpvt[i] = jpvt[i], jpvt[p] vn1[p] = vn1[i] vn2[p] = vn2[i] } // Generate elementary reflector H_i. if offpi < m-1 { a[offpi*lda+i], tau[i] = impl.Dlarfg(m-offpi, a[offpi*lda+i], a[(offpi+1)*lda+i:], lda) } else { tau[i] = 0 } if i < n-1 { // Apply H_iᵀ to A[offset+i:m, i:n] from the left. aii := a[offpi*lda+i] a[offpi*lda+i] = 1 impl.Dlarf(blas.Left, m-offpi, n-i-1, a[offpi*lda+i:], lda, tau[i], a[offpi*lda+i+1:], lda, work) a[offpi*lda+i] = aii } // Update partial column norms. for j := i + 1; j < n; j++ { if vn1[j] == 0 { continue } // The following marked lines follow from the // analysis in Lapack Working Note 176. r := math.Abs(a[offpi*lda+j]) / vn1[j] // * temp := math.Max(0, 1-r*r) // * r = vn1[j] / vn2[j] // * temp2 := temp * r * r // * if temp2 < tol3z { var v float64 if offpi < m-1 { v = bi.Dnrm2(m-offpi-1, a[(offpi+1)*lda+j:], lda) } vn1[j] = v vn2[j] = v } else { vn1[j] *= math.Sqrt(temp) // * } } } }