| // Copyright ©2016 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" |
| |
| // Dlarfx applies an elementary reflector H to a real m×n matrix C, from either |
| // the left or the right, with loop unrolling when the reflector has order less |
| // than 11. |
| // |
| // H is represented in the form |
| // H = I - tau * v * v^T, |
| // where tau is a real scalar and v is a real vector. If tau = 0, then H is |
| // taken to be the identity matrix. |
| // |
| // v must have length equal to m if side == blas.Left, and equal to n if side == |
| // blas.Right, otherwise Dlarfx will panic. |
| // |
| // c and ldc represent the m×n matrix C. On return, C is overwritten by the |
| // matrix H * C if side == blas.Left, or C * H if side == blas.Right. |
| // |
| // work must have length at least n if side == blas.Left, and at least m if side |
| // == blas.Right, otherwise Dlarfx will panic. work is not referenced if H has |
| // order < 11. |
| // |
| // Dlarfx is an internal routine. It is exported for testing purposes. |
| func (impl Implementation) Dlarfx(side blas.Side, m, n int, v []float64, tau float64, c []float64, ldc int, work []float64) { |
| checkMatrix(m, n, c, ldc) |
| switch side { |
| case blas.Left: |
| checkVector(m, v, 1) |
| if m > 10 && len(work) < n { |
| panic(badWork) |
| } |
| case blas.Right: |
| checkVector(n, v, 1) |
| if n > 10 && len(work) < m { |
| panic(badWork) |
| } |
| default: |
| panic(badSide) |
| } |
| |
| if tau == 0 { |
| return |
| } |
| |
| if side == blas.Left { |
| // Form H * C, where H has order m. |
| switch m { |
| default: // Code for general m. |
| impl.Dlarf(side, m, n, v, 1, tau, c, ldc, work) |
| return |
| |
| case 0: // No-op for zero size matrix. |
| return |
| |
| case 1: // Special code for 1×1 Householder matrix. |
| t0 := 1 - tau*v[0]*v[0] |
| for j := 0; j < n; j++ { |
| c[j] *= t0 |
| } |
| return |
| |
| case 2: // Special code for 2×2 Householder matrix. |
| v0 := v[0] |
| t0 := tau * v0 |
| v1 := v[1] |
| t1 := tau * v1 |
| for j := 0; j < n; j++ { |
| sum := v0*c[j] + v1*c[ldc+j] |
| c[j] -= sum * t0 |
| c[ldc+j] -= sum * t1 |
| } |
| return |
| |
| case 3: // Special code for 3×3 Householder matrix. |
| v0 := v[0] |
| t0 := tau * v0 |
| v1 := v[1] |
| t1 := tau * v1 |
| v2 := v[2] |
| t2 := tau * v2 |
| for j := 0; j < n; j++ { |
| sum := v0*c[j] + v1*c[ldc+j] + v2*c[2*ldc+j] |
| c[j] -= sum * t0 |
| c[ldc+j] -= sum * t1 |
| c[2*ldc+j] -= sum * t2 |
| } |
| return |
| |
| case 4: // Special code for 4×4 Householder matrix. |
| v0 := v[0] |
| t0 := tau * v0 |
| v1 := v[1] |
| t1 := tau * v1 |
| v2 := v[2] |
| t2 := tau * v2 |
| v3 := v[3] |
| t3 := tau * v3 |
| for j := 0; j < n; j++ { |
| sum := v0*c[j] + v1*c[ldc+j] + v2*c[2*ldc+j] + v3*c[3*ldc+j] |
| c[j] -= sum * t0 |
| c[ldc+j] -= sum * t1 |
| c[2*ldc+j] -= sum * t2 |
| c[3*ldc+j] -= sum * t3 |
| } |
| return |
| |
| case 5: // Special code for 5×5 Householder matrix. |
| v0 := v[0] |
| t0 := tau * v0 |
| v1 := v[1] |
| t1 := tau * v1 |
| v2 := v[2] |
| t2 := tau * v2 |
| v3 := v[3] |
| t3 := tau * v3 |
| v4 := v[4] |
| t4 := tau * v4 |
| for j := 0; j < n; j++ { |
| sum := v0*c[j] + v1*c[ldc+j] + v2*c[2*ldc+j] + v3*c[3*ldc+j] + v4*c[4*ldc+j] |
| c[j] -= sum * t0 |
| c[ldc+j] -= sum * t1 |
| c[2*ldc+j] -= sum * t2 |
| c[3*ldc+j] -= sum * t3 |
| c[4*ldc+j] -= sum * t4 |
| } |
| return |
| |
| case 6: // Special code for 6×6 Householder matrix. |
| v0 := v[0] |
| t0 := tau * v0 |
| v1 := v[1] |
| t1 := tau * v1 |
| v2 := v[2] |
| t2 := tau * v2 |
| v3 := v[3] |
| t3 := tau * v3 |
| v4 := v[4] |
| t4 := tau * v4 |
| v5 := v[5] |
| t5 := tau * v5 |
| for j := 0; j < n; j++ { |
| sum := v0*c[j] + v1*c[ldc+j] + v2*c[2*ldc+j] + v3*c[3*ldc+j] + v4*c[4*ldc+j] + |
| v5*c[5*ldc+j] |
| c[j] -= sum * t0 |
| c[ldc+j] -= sum * t1 |
| c[2*ldc+j] -= sum * t2 |
| c[3*ldc+j] -= sum * t3 |
| c[4*ldc+j] -= sum * t4 |
| c[5*ldc+j] -= sum * t5 |
| } |
| return |
| |
| case 7: // Special code for 7×7 Householder matrix. |
| v0 := v[0] |
| t0 := tau * v0 |
| v1 := v[1] |
| t1 := tau * v1 |
| v2 := v[2] |
| t2 := tau * v2 |
| v3 := v[3] |
| t3 := tau * v3 |
| v4 := v[4] |
| t4 := tau * v4 |
| v5 := v[5] |
| t5 := tau * v5 |
| v6 := v[6] |
| t6 := tau * v6 |
| for j := 0; j < n; j++ { |
| sum := v0*c[j] + v1*c[ldc+j] + v2*c[2*ldc+j] + v3*c[3*ldc+j] + v4*c[4*ldc+j] + |
| v5*c[5*ldc+j] + v6*c[6*ldc+j] |
| c[j] -= sum * t0 |
| c[ldc+j] -= sum * t1 |
| c[2*ldc+j] -= sum * t2 |
| c[3*ldc+j] -= sum * t3 |
| c[4*ldc+j] -= sum * t4 |
| c[5*ldc+j] -= sum * t5 |
| c[6*ldc+j] -= sum * t6 |
| } |
| return |
| |
| case 8: // Special code for 8×8 Householder matrix. |
| v0 := v[0] |
| t0 := tau * v0 |
| v1 := v[1] |
| t1 := tau * v1 |
| v2 := v[2] |
| t2 := tau * v2 |
| v3 := v[3] |
| t3 := tau * v3 |
| v4 := v[4] |
| t4 := tau * v4 |
| v5 := v[5] |
| t5 := tau * v5 |
| v6 := v[6] |
| t6 := tau * v6 |
| v7 := v[7] |
| t7 := tau * v7 |
| for j := 0; j < n; j++ { |
| sum := v0*c[j] + v1*c[ldc+j] + v2*c[2*ldc+j] + v3*c[3*ldc+j] + v4*c[4*ldc+j] + |
| v5*c[5*ldc+j] + v6*c[6*ldc+j] + v7*c[7*ldc+j] |
| c[j] -= sum * t0 |
| c[ldc+j] -= sum * t1 |
| c[2*ldc+j] -= sum * t2 |
| c[3*ldc+j] -= sum * t3 |
| c[4*ldc+j] -= sum * t4 |
| c[5*ldc+j] -= sum * t5 |
| c[6*ldc+j] -= sum * t6 |
| c[7*ldc+j] -= sum * t7 |
| } |
| return |
| |
| case 9: // Special code for 9×9 Householder matrix. |
| v0 := v[0] |
| t0 := tau * v0 |
| v1 := v[1] |
| t1 := tau * v1 |
| v2 := v[2] |
| t2 := tau * v2 |
| v3 := v[3] |
| t3 := tau * v3 |
| v4 := v[4] |
| t4 := tau * v4 |
| v5 := v[5] |
| t5 := tau * v5 |
| v6 := v[6] |
| t6 := tau * v6 |
| v7 := v[7] |
| t7 := tau * v7 |
| v8 := v[8] |
| t8 := tau * v8 |
| for j := 0; j < n; j++ { |
| sum := v0*c[j] + v1*c[ldc+j] + v2*c[2*ldc+j] + v3*c[3*ldc+j] + v4*c[4*ldc+j] + |
| v5*c[5*ldc+j] + v6*c[6*ldc+j] + v7*c[7*ldc+j] + v8*c[8*ldc+j] |
| c[j] -= sum * t0 |
| c[ldc+j] -= sum * t1 |
| c[2*ldc+j] -= sum * t2 |
| c[3*ldc+j] -= sum * t3 |
| c[4*ldc+j] -= sum * t4 |
| c[5*ldc+j] -= sum * t5 |
| c[6*ldc+j] -= sum * t6 |
| c[7*ldc+j] -= sum * t7 |
| c[8*ldc+j] -= sum * t8 |
| } |
| return |
| |
| case 10: // Special code for 10×10 Householder matrix. |
| v0 := v[0] |
| t0 := tau * v0 |
| v1 := v[1] |
| t1 := tau * v1 |
| v2 := v[2] |
| t2 := tau * v2 |
| v3 := v[3] |
| t3 := tau * v3 |
| v4 := v[4] |
| t4 := tau * v4 |
| v5 := v[5] |
| t5 := tau * v5 |
| v6 := v[6] |
| t6 := tau * v6 |
| v7 := v[7] |
| t7 := tau * v7 |
| v8 := v[8] |
| t8 := tau * v8 |
| v9 := v[9] |
| t9 := tau * v9 |
| for j := 0; j < n; j++ { |
| sum := v0*c[j] + v1*c[ldc+j] + v2*c[2*ldc+j] + v3*c[3*ldc+j] + v4*c[4*ldc+j] + |
| 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] |
| c[j] -= sum * t0 |
| c[ldc+j] -= sum * t1 |
| c[2*ldc+j] -= sum * t2 |
| c[3*ldc+j] -= sum * t3 |
| c[4*ldc+j] -= sum * t4 |
| c[5*ldc+j] -= sum * t5 |
| c[6*ldc+j] -= sum * t6 |
| c[7*ldc+j] -= sum * t7 |
| c[8*ldc+j] -= sum * t8 |
| c[9*ldc+j] -= sum * t9 |
| } |
| return |
| } |
| } |
| |
| // Form C * H, where H has order n. |
| switch n { |
| default: // Code for general n. |
| impl.Dlarf(side, m, n, v, 1, tau, c, ldc, work) |
| return |
| |
| case 0: // No-op for zero size matrix. |
| return |
| |
| case 1: // Special code for 1×1 Householder matrix. |
| t0 := 1 - tau*v[0]*v[0] |
| for j := 0; j < m; j++ { |
| c[j*ldc] *= t0 |
| } |
| return |
| |
| case 2: // Special code for 2×2 Householder matrix. |
| v0 := v[0] |
| t0 := tau * v0 |
| v1 := v[1] |
| t1 := tau * v1 |
| for j := 0; j < m; j++ { |
| cs := c[j*ldc:] |
| sum := v0*cs[0] + v1*cs[1] |
| cs[0] -= sum * t0 |
| cs[1] -= sum * t1 |
| } |
| return |
| |
| case 3: // Special code for 3×3 Householder matrix. |
| v0 := v[0] |
| t0 := tau * v0 |
| v1 := v[1] |
| t1 := tau * v1 |
| v2 := v[2] |
| t2 := tau * v2 |
| for j := 0; j < m; j++ { |
| cs := c[j*ldc:] |
| sum := v0*cs[0] + v1*cs[1] + v2*cs[2] |
| cs[0] -= sum * t0 |
| cs[1] -= sum * t1 |
| cs[2] -= sum * t2 |
| } |
| return |
| |
| case 4: // Special code for 4×4 Householder matrix. |
| v0 := v[0] |
| t0 := tau * v0 |
| v1 := v[1] |
| t1 := tau * v1 |
| v2 := v[2] |
| t2 := tau * v2 |
| v3 := v[3] |
| t3 := tau * v3 |
| for j := 0; j < m; j++ { |
| cs := c[j*ldc:] |
| sum := v0*cs[0] + v1*cs[1] + v2*cs[2] + v3*cs[3] |
| cs[0] -= sum * t0 |
| cs[1] -= sum * t1 |
| cs[2] -= sum * t2 |
| cs[3] -= sum * t3 |
| } |
| return |
| |
| case 5: // Special code for 5×5 Householder matrix. |
| v0 := v[0] |
| t0 := tau * v0 |
| v1 := v[1] |
| t1 := tau * v1 |
| v2 := v[2] |
| t2 := tau * v2 |
| v3 := v[3] |
| t3 := tau * v3 |
| v4 := v[4] |
| t4 := tau * v4 |
| for j := 0; j < m; j++ { |
| cs := c[j*ldc:] |
| sum := v0*cs[0] + v1*cs[1] + v2*cs[2] + v3*cs[3] + v4*cs[4] |
| cs[0] -= sum * t0 |
| cs[1] -= sum * t1 |
| cs[2] -= sum * t2 |
| cs[3] -= sum * t3 |
| cs[4] -= sum * t4 |
| } |
| return |
| |
| case 6: // Special code for 6×6 Householder matrix. |
| v0 := v[0] |
| t0 := tau * v0 |
| v1 := v[1] |
| t1 := tau * v1 |
| v2 := v[2] |
| t2 := tau * v2 |
| v3 := v[3] |
| t3 := tau * v3 |
| v4 := v[4] |
| t4 := tau * v4 |
| v5 := v[5] |
| t5 := tau * v5 |
| for j := 0; j < m; j++ { |
| cs := c[j*ldc:] |
| sum := v0*cs[0] + v1*cs[1] + v2*cs[2] + v3*cs[3] + v4*cs[4] + v5*cs[5] |
| cs[0] -= sum * t0 |
| cs[1] -= sum * t1 |
| cs[2] -= sum * t2 |
| cs[3] -= sum * t3 |
| cs[4] -= sum * t4 |
| cs[5] -= sum * t5 |
| } |
| return |
| |
| case 7: // Special code for 7×7 Householder matrix. |
| v0 := v[0] |
| t0 := tau * v0 |
| v1 := v[1] |
| t1 := tau * v1 |
| v2 := v[2] |
| t2 := tau * v2 |
| v3 := v[3] |
| t3 := tau * v3 |
| v4 := v[4] |
| t4 := tau * v4 |
| v5 := v[5] |
| t5 := tau * v5 |
| v6 := v[6] |
| t6 := tau * v6 |
| for j := 0; j < m; j++ { |
| cs := c[j*ldc:] |
| sum := v0*cs[0] + v1*cs[1] + v2*cs[2] + v3*cs[3] + v4*cs[4] + |
| v5*cs[5] + v6*cs[6] |
| cs[0] -= sum * t0 |
| cs[1] -= sum * t1 |
| cs[2] -= sum * t2 |
| cs[3] -= sum * t3 |
| cs[4] -= sum * t4 |
| cs[5] -= sum * t5 |
| cs[6] -= sum * t6 |
| } |
| return |
| |
| case 8: // Special code for 8×8 Householder matrix. |
| v0 := v[0] |
| t0 := tau * v0 |
| v1 := v[1] |
| t1 := tau * v1 |
| v2 := v[2] |
| t2 := tau * v2 |
| v3 := v[3] |
| t3 := tau * v3 |
| v4 := v[4] |
| t4 := tau * v4 |
| v5 := v[5] |
| t5 := tau * v5 |
| v6 := v[6] |
| t6 := tau * v6 |
| v7 := v[7] |
| t7 := tau * v7 |
| for j := 0; j < m; j++ { |
| cs := c[j*ldc:] |
| sum := v0*cs[0] + v1*cs[1] + v2*cs[2] + v3*cs[3] + v4*cs[4] + |
| v5*cs[5] + v6*cs[6] + v7*cs[7] |
| cs[0] -= sum * t0 |
| cs[1] -= sum * t1 |
| cs[2] -= sum * t2 |
| cs[3] -= sum * t3 |
| cs[4] -= sum * t4 |
| cs[5] -= sum * t5 |
| cs[6] -= sum * t6 |
| cs[7] -= sum * t7 |
| } |
| return |
| |
| case 9: // Special code for 9×9 Householder matrix. |
| v0 := v[0] |
| t0 := tau * v0 |
| v1 := v[1] |
| t1 := tau * v1 |
| v2 := v[2] |
| t2 := tau * v2 |
| v3 := v[3] |
| t3 := tau * v3 |
| v4 := v[4] |
| t4 := tau * v4 |
| v5 := v[5] |
| t5 := tau * v5 |
| v6 := v[6] |
| t6 := tau * v6 |
| v7 := v[7] |
| t7 := tau * v7 |
| v8 := v[8] |
| t8 := tau * v8 |
| for j := 0; j < m; j++ { |
| cs := c[j*ldc:] |
| sum := v0*cs[0] + v1*cs[1] + v2*cs[2] + v3*cs[3] + v4*cs[4] + |
| v5*cs[5] + v6*cs[6] + v7*cs[7] + v8*cs[8] |
| cs[0] -= sum * t0 |
| cs[1] -= sum * t1 |
| cs[2] -= sum * t2 |
| cs[3] -= sum * t3 |
| cs[4] -= sum * t4 |
| cs[5] -= sum * t5 |
| cs[6] -= sum * t6 |
| cs[7] -= sum * t7 |
| cs[8] -= sum * t8 |
| } |
| return |
| |
| case 10: // Special code for 10×10 Householder matrix. |
| v0 := v[0] |
| t0 := tau * v0 |
| v1 := v[1] |
| t1 := tau * v1 |
| v2 := v[2] |
| t2 := tau * v2 |
| v3 := v[3] |
| t3 := tau * v3 |
| v4 := v[4] |
| t4 := tau * v4 |
| v5 := v[5] |
| t5 := tau * v5 |
| v6 := v[6] |
| t6 := tau * v6 |
| v7 := v[7] |
| t7 := tau * v7 |
| v8 := v[8] |
| t8 := tau * v8 |
| v9 := v[9] |
| t9 := tau * v9 |
| for j := 0; j < m; j++ { |
| cs := c[j*ldc:] |
| sum := v0*cs[0] + v1*cs[1] + v2*cs[2] + v3*cs[3] + v4*cs[4] + |
| v5*cs[5] + v6*cs[6] + v7*cs[7] + v8*cs[8] + v9*cs[9] |
| cs[0] -= sum * t0 |
| cs[1] -= sum * t1 |
| cs[2] -= sum * t2 |
| cs[3] -= sum * t3 |
| cs[4] -= sum * t4 |
| cs[5] -= sum * t5 |
| cs[6] -= sum * t6 |
| cs[7] -= sum * t7 |
| cs[8] -= sum * t8 |
| cs[9] -= sum * t9 |
| } |
| return |
| } |
| } |