1 // Copyright ©2016 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.
7 import "gonum.org/v1/gonum/blas"
9 // Dlarfx applies an elementary reflector H to a real m×n matrix C, from either
10 // the left or the right, with loop unrolling when the reflector has order less
13 // H is represented in the form
14 // H = I - tau * v * v^T,
15 // where tau is a real scalar and v is a real vector. If tau = 0, then H is
16 // taken to be the identity matrix.
18 // v must have length equal to m if side == blas.Left, and equal to n if side ==
19 // blas.Right, otherwise Dlarfx will panic.
21 // c and ldc represent the m×n matrix C. On return, C is overwritten by the
22 // matrix H * C if side == blas.Left, or C * H if side == blas.Right.
24 // work must have length at least n if side == blas.Left, and at least m if side
25 // == blas.Right, otherwise Dlarfx will panic. work is not referenced if H has
28 // Dlarfx is an internal routine. It is exported for testing purposes.
29 func (impl Implementation) Dlarfx(side blas.Side, m, n int, v []float64, tau float64, c []float64, ldc int, work []float64) {
30 checkMatrix(m, n, c, ldc)
34 if m > 10 && len(work) < n {
39 if n > 10 && len(work) < m {
50 if side == blas.Left {
51 // Form H * C, where H has order m.
53 default: // Code for general m.
54 impl.Dlarf(side, m, n, v, 1, tau, c, ldc, work)
57 case 0: // No-op for zero size matrix.
60 case 1: // Special code for 1×1 Householder matrix.
61 t0 := 1 - tau*v[0]*v[0]
62 for j := 0; j < n; j++ {
67 case 2: // Special code for 2×2 Householder matrix.
72 for j := 0; j < n; j++ {
73 sum := v0*c[j] + v1*c[ldc+j]
79 case 3: // Special code for 3×3 Householder matrix.
86 for j := 0; j < n; j++ {
87 sum := v0*c[j] + v1*c[ldc+j] + v2*c[2*ldc+j]
90 c[2*ldc+j] -= sum * t2
94 case 4: // Special code for 4×4 Householder matrix.
103 for j := 0; j < n; j++ {
104 sum := v0*c[j] + v1*c[ldc+j] + v2*c[2*ldc+j] + v3*c[3*ldc+j]
107 c[2*ldc+j] -= sum * t2
108 c[3*ldc+j] -= sum * t3
112 case 5: // Special code for 5×5 Householder matrix.
123 for j := 0; j < n; j++ {
124 sum := v0*c[j] + v1*c[ldc+j] + v2*c[2*ldc+j] + v3*c[3*ldc+j] + v4*c[4*ldc+j]
127 c[2*ldc+j] -= sum * t2
128 c[3*ldc+j] -= sum * t3
129 c[4*ldc+j] -= sum * t4
133 case 6: // Special code for 6×6 Householder matrix.
146 for j := 0; j < n; j++ {
147 sum := v0*c[j] + v1*c[ldc+j] + v2*c[2*ldc+j] + v3*c[3*ldc+j] + v4*c[4*ldc+j] +
151 c[2*ldc+j] -= sum * t2
152 c[3*ldc+j] -= sum * t3
153 c[4*ldc+j] -= sum * t4
154 c[5*ldc+j] -= sum * t5
158 case 7: // Special code for 7×7 Householder matrix.
173 for j := 0; j < n; j++ {
174 sum := v0*c[j] + v1*c[ldc+j] + v2*c[2*ldc+j] + v3*c[3*ldc+j] + v4*c[4*ldc+j] +
175 v5*c[5*ldc+j] + v6*c[6*ldc+j]
178 c[2*ldc+j] -= sum * t2
179 c[3*ldc+j] -= sum * t3
180 c[4*ldc+j] -= sum * t4
181 c[5*ldc+j] -= sum * t5
182 c[6*ldc+j] -= sum * t6
186 case 8: // Special code for 8×8 Householder matrix.
203 for j := 0; j < n; j++ {
204 sum := v0*c[j] + v1*c[ldc+j] + v2*c[2*ldc+j] + v3*c[3*ldc+j] + v4*c[4*ldc+j] +
205 v5*c[5*ldc+j] + v6*c[6*ldc+j] + v7*c[7*ldc+j]
208 c[2*ldc+j] -= sum * t2
209 c[3*ldc+j] -= sum * t3
210 c[4*ldc+j] -= sum * t4
211 c[5*ldc+j] -= sum * t5
212 c[6*ldc+j] -= sum * t6
213 c[7*ldc+j] -= sum * t7
217 case 9: // Special code for 9×9 Householder matrix.
236 for j := 0; j < n; j++ {
237 sum := v0*c[j] + v1*c[ldc+j] + v2*c[2*ldc+j] + v3*c[3*ldc+j] + v4*c[4*ldc+j] +
238 v5*c[5*ldc+j] + v6*c[6*ldc+j] + v7*c[7*ldc+j] + v8*c[8*ldc+j]
241 c[2*ldc+j] -= sum * t2
242 c[3*ldc+j] -= sum * t3
243 c[4*ldc+j] -= sum * t4
244 c[5*ldc+j] -= sum * t5
245 c[6*ldc+j] -= sum * t6
246 c[7*ldc+j] -= sum * t7
247 c[8*ldc+j] -= sum * t8
251 case 10: // Special code for 10×10 Householder matrix.
272 for j := 0; j < n; j++ {
273 sum := v0*c[j] + v1*c[ldc+j] + v2*c[2*ldc+j] + v3*c[3*ldc+j] + v4*c[4*ldc+j] +
274 v5*c[5*ldc+j] + v6*c[6*ldc+j] + v7*c[7*ldc+j] + v8*c[8*ldc+j] + v9*c[9*ldc+j]
277 c[2*ldc+j] -= sum * t2
278 c[3*ldc+j] -= sum * t3
279 c[4*ldc+j] -= sum * t4
280 c[5*ldc+j] -= sum * t5
281 c[6*ldc+j] -= sum * t6
282 c[7*ldc+j] -= sum * t7
283 c[8*ldc+j] -= sum * t8
284 c[9*ldc+j] -= sum * t9
290 // Form C * H, where H has order n.
292 default: // Code for general n.
293 impl.Dlarf(side, m, n, v, 1, tau, c, ldc, work)
296 case 0: // No-op for zero size matrix.
299 case 1: // Special code for 1×1 Householder matrix.
300 t0 := 1 - tau*v[0]*v[0]
301 for j := 0; j < m; j++ {
306 case 2: // Special code for 2×2 Householder matrix.
311 for j := 0; j < m; j++ {
313 sum := v0*cs[0] + v1*cs[1]
319 case 3: // Special code for 3×3 Householder matrix.
326 for j := 0; j < m; j++ {
328 sum := v0*cs[0] + v1*cs[1] + v2*cs[2]
335 case 4: // Special code for 4×4 Householder matrix.
344 for j := 0; j < m; j++ {
346 sum := v0*cs[0] + v1*cs[1] + v2*cs[2] + v3*cs[3]
354 case 5: // Special code for 5×5 Householder matrix.
365 for j := 0; j < m; j++ {
367 sum := v0*cs[0] + v1*cs[1] + v2*cs[2] + v3*cs[3] + v4*cs[4]
376 case 6: // Special code for 6×6 Householder matrix.
389 for j := 0; j < m; j++ {
391 sum := v0*cs[0] + v1*cs[1] + v2*cs[2] + v3*cs[3] + v4*cs[4] + v5*cs[5]
401 case 7: // Special code for 7×7 Householder matrix.
416 for j := 0; j < m; j++ {
418 sum := v0*cs[0] + v1*cs[1] + v2*cs[2] + v3*cs[3] + v4*cs[4] +
430 case 8: // Special code for 8×8 Householder matrix.
447 for j := 0; j < m; j++ {
449 sum := v0*cs[0] + v1*cs[1] + v2*cs[2] + v3*cs[3] + v4*cs[4] +
450 v5*cs[5] + v6*cs[6] + v7*cs[7]
462 case 9: // Special code for 9×9 Householder matrix.
481 for j := 0; j < m; j++ {
483 sum := v0*cs[0] + v1*cs[1] + v2*cs[2] + v3*cs[3] + v4*cs[4] +
484 v5*cs[5] + v6*cs[6] + v7*cs[7] + v8*cs[8]
497 case 10: // Special code for 10×10 Householder matrix.
518 for j := 0; j < m; j++ {
520 sum := v0*cs[0] + v1*cs[1] + v2*cs[2] + v3*cs[3] + v4*cs[4] +
521 v5*cs[5] + v6*cs[6] + v7*cs[7] + v8*cs[8] + v9*cs[9]