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
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
|
// Copyright ©2015 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 (
"gonum.org/v1/gonum/blas"
"gonum.org/v1/gonum/blas/blas64"
)
// Dlabrd reduces the first NB rows and columns of a real general m×n matrix
// A to upper or lower bidiagonal form by an orthogonal transformation
//
// Q**T * A * P
//
// If m >= n, A is reduced to upper bidiagonal form and upon exit the elements
// on and below the diagonal in the first nb columns represent the elementary
// reflectors, and the elements above the diagonal in the first nb rows represent
// the matrix P. If m < n, A is reduced to lower bidiagonal form and the elements
// P is instead stored above the diagonal.
//
// The reduction to bidiagonal form is stored in d and e, where d are the diagonal
// elements, and e are the off-diagonal elements.
//
// The matrices Q and P are products of elementary reflectors
//
// Q = H_0 * H_1 * ... * H_{nb-1}
// P = G_0 * G_1 * ... * G_{nb-1}
//
// where
//
// H_i = I - tauQ[i] * v_i * v_iᵀ
// G_i = I - tauP[i] * u_i * u_iᵀ
//
// As an example, on exit the entries of A when m = 6, n = 5, and nb = 2
//
// [ 1 1 u1 u1 u1]
// [v1 1 1 u2 u2]
// [v1 v2 a a a]
// [v1 v2 a a a]
// [v1 v2 a a a]
// [v1 v2 a a a]
//
// and when m = 5, n = 6, and nb = 2
//
// [ 1 u1 u1 u1 u1 u1]
// [ 1 1 u2 u2 u2 u2]
// [v1 1 a a a a]
// [v1 v2 a a a a]
// [v1 v2 a a a a]
//
// Dlabrd also returns the matrices X and Y which are used with U and V to
// apply the transformation to the unreduced part of the matrix
//
// A := A - V*Yᵀ - X*Uᵀ
//
// and returns the matrices X and Y which are needed to apply the
// transformation to the unreduced part of A.
//
// X is an m×nb matrix, Y is an n×nb matrix. d, e, taup, and tauq must all have
// length at least nb. Dlabrd will panic if these size constraints are violated.
//
// Dlabrd is an internal routine. It is exported for testing purposes.
func (impl Implementation) Dlabrd(m, n, nb int, a []float64, lda int, d, e, tauQ, tauP, x []float64, ldx int, y []float64, ldy int) {
switch {
case m < 0:
panic(mLT0)
case n < 0:
panic(nLT0)
case nb < 0:
panic(nbLT0)
case nb > n:
panic(nbGTN)
case nb > m:
panic(nbGTM)
case lda < max(1, n):
panic(badLdA)
case ldx < max(1, nb):
panic(badLdX)
case ldy < max(1, nb):
panic(badLdY)
}
if m == 0 || n == 0 || nb == 0 {
return
}
switch {
case len(a) < (m-1)*lda+n:
panic(shortA)
case len(d) < nb:
panic(shortD)
case len(e) < nb:
panic(shortE)
case len(tauQ) < nb:
panic(shortTauQ)
case len(tauP) < nb:
panic(shortTauP)
case len(x) < (m-1)*ldx+nb:
panic(shortX)
case len(y) < (n-1)*ldy+nb:
panic(shortY)
}
bi := blas64.Implementation()
if m >= n {
// Reduce to upper bidiagonal form.
for i := 0; i < nb; i++ {
bi.Dgemv(blas.NoTrans, m-i, i, -1, a[i*lda:], lda, y[i*ldy:], 1, 1, a[i*lda+i:], lda)
bi.Dgemv(blas.NoTrans, m-i, i, -1, x[i*ldx:], ldx, a[i:], lda, 1, a[i*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]
if i < n-1 {
// Compute Y[i+1:n, i].
a[i*lda+i] = 1
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)
bi.Dgemv(blas.Trans, m-i, i, 1, a[i*lda:], lda, a[i*lda+i:], lda, 0, y[i:], ldy)
bi.Dgemv(blas.NoTrans, n-i-1, i, -1, y[(i+1)*ldy:], ldy, y[i:], ldy, 1, y[(i+1)*ldy+i:], ldy)
bi.Dgemv(blas.Trans, m-i, i, 1, x[i*ldx:], ldx, a[i*lda+i:], lda, 0, y[i:], ldy)
bi.Dgemv(blas.Trans, i, n-i-1, -1, a[i+1:], lda, y[i:], ldy, 1, y[(i+1)*ldy+i:], ldy)
bi.Dscal(n-i-1, tauQ[i], y[(i+1)*ldy+i:], ldy)
// Update A[i, i+1:n].
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)
bi.Dgemv(blas.Trans, i, n-i-1, -1, a[i+1:], lda, x[i*ldx:], 1, 1, a[i*lda+i+1:], 1)
// Generate reflection P[i] to annihilate A[i, i+2:n].
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)
e[i] = a[i*lda+i+1]
a[i*lda+i+1] = 1
// Compute X[i+1:m, i].
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)
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)
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)
bi.Dgemv(blas.NoTrans, i, n-i-1, 1, a[i+1:], lda, a[i*lda+i+1:], 1, 0, x[i:], ldx)
bi.Dgemv(blas.NoTrans, m-i-1, i, -1, x[(i+1)*ldx:], ldx, x[i:], ldx, 1, x[(i+1)*ldx+i:], ldx)
bi.Dscal(m-i-1, tauP[i], x[(i+1)*ldx+i:], ldx)
}
}
return
}
// Reduce to lower bidiagonal form.
for i := 0; i < nb; i++ {
// Update A[i,i:n]
bi.Dgemv(blas.NoTrans, n-i, i, -1, y[i*ldy:], ldy, a[i*lda:], 1, 1, a[i*lda+i:], 1)
bi.Dgemv(blas.Trans, i, n-i, -1, a[i:], lda, x[i*ldx:], 1, 1, a[i*lda+i:], 1)
// Generate reflection P[i] to annihilate A[i, i+1:n]
a[i*lda+i], tauP[i] = impl.Dlarfg(n-i, a[i*lda+i], a[i*lda+min(i+1, n-1):], 1)
d[i] = a[i*lda+i]
if i < m-1 {
a[i*lda+i] = 1
// Compute X[i+1:m, i].
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)
bi.Dgemv(blas.Trans, n-i, i, 1, y[i*ldy:], ldy, a[i*lda+i:], 1, 0, x[i:], ldx)
bi.Dgemv(blas.NoTrans, m-i-1, i, -1, a[(i+1)*lda:], lda, x[i:], ldx, 1, x[(i+1)*ldx+i:], ldx)
bi.Dgemv(blas.NoTrans, i, n-i, 1, a[i:], lda, a[i*lda+i:], 1, 0, x[i:], ldx)
bi.Dgemv(blas.NoTrans, m-i-1, i, -1, x[(i+1)*ldx:], ldx, x[i:], ldx, 1, x[(i+1)*ldx+i:], ldx)
bi.Dscal(m-i-1, tauP[i], x[(i+1)*ldx+i:], ldx)
// Update A[i+1:m, i].
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)
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)
// Generate reflection Q[i] to annihilate A[i+2:m, i].
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)
e[i] = a[(i+1)*lda+i]
a[(i+1)*lda+i] = 1
// Compute Y[i+1:n, i].
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)
bi.Dgemv(blas.Trans, m-i-1, i, 1, a[(i+1)*lda:], lda, a[(i+1)*lda+i:], lda, 0, y[i:], ldy)
bi.Dgemv(blas.NoTrans, n-i-1, i, -1, y[(i+1)*ldy:], ldy, y[i:], ldy, 1, y[(i+1)*ldy+i:], ldy)
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)
bi.Dgemv(blas.Trans, i+1, n-i-1, -1, a[i+1:], lda, y[i:], ldy, 1, y[(i+1)*ldy+i:], ldy)
bi.Dscal(n-i-1, tauQ[i], y[(i+1)*ldy+i:], ldy)
}
}
}
|
