1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
|
// 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) // *
}
}
}
}
|
