1// Copyright ©2017 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. 4 5// +build !amd64 noasm gccgo safe 6 7package f64 8 9// Ger performs the rank-one operation 10// A += alpha * x * yᵀ 11// where A is an m×n dense matrix, x and y are vectors, and alpha is a scalar. 12func Ger(m, n uintptr, alpha float64, x []float64, incX uintptr, y []float64, incY uintptr, a []float64, lda uintptr) { 13 if incX == 1 && incY == 1 { 14 x = x[:m] 15 y = y[:n] 16 for i, xv := range x { 17 AxpyUnitary(alpha*xv, y, a[uintptr(i)*lda:uintptr(i)*lda+n]) 18 } 19 return 20 } 21 22 var ky, kx uintptr 23 if int(incY) < 0 { 24 ky = uintptr(-int(n-1) * int(incY)) 25 } 26 if int(incX) < 0 { 27 kx = uintptr(-int(m-1) * int(incX)) 28 } 29 30 ix := kx 31 for i := 0; i < int(m); i++ { 32 AxpyInc(alpha*x[ix], y, a[uintptr(i)*lda:uintptr(i)*lda+n], n, incY, 1, ky, 0) 33 ix += incX 34 } 35} 36 37// GemvN computes 38// y = alpha * A * x + beta * y 39// where A is an m×n dense matrix, x and y are vectors, and alpha and beta are scalars. 40func GemvN(m, n uintptr, alpha float64, a []float64, lda uintptr, x []float64, incX uintptr, beta float64, y []float64, incY uintptr) { 41 var kx, ky, i uintptr 42 if int(incX) < 0 { 43 kx = uintptr(-int(n-1) * int(incX)) 44 } 45 if int(incY) < 0 { 46 ky = uintptr(-int(m-1) * int(incY)) 47 } 48 49 if incX == 1 && incY == 1 { 50 if beta == 0 { 51 for i = 0; i < m; i++ { 52 y[i] = alpha * DotUnitary(a[lda*i:lda*i+n], x) 53 } 54 return 55 } 56 for i = 0; i < m; i++ { 57 y[i] = y[i]*beta + alpha*DotUnitary(a[lda*i:lda*i+n], x) 58 } 59 return 60 } 61 iy := ky 62 if beta == 0 { 63 for i = 0; i < m; i++ { 64 y[iy] = alpha * DotInc(x, a[lda*i:lda*i+n], n, incX, 1, kx, 0) 65 iy += incY 66 } 67 return 68 } 69 for i = 0; i < m; i++ { 70 y[iy] = y[iy]*beta + alpha*DotInc(x, a[lda*i:lda*i+n], n, incX, 1, kx, 0) 71 iy += incY 72 } 73} 74 75// GemvT computes 76// y = alpha * Aᵀ * x + beta * y 77// where A is an m×n dense matrix, x and y are vectors, and alpha and beta are scalars. 78func GemvT(m, n uintptr, alpha float64, a []float64, lda uintptr, x []float64, incX uintptr, beta float64, y []float64, incY uintptr) { 79 var kx, ky, i uintptr 80 if int(incX) < 0 { 81 kx = uintptr(-int(m-1) * int(incX)) 82 } 83 if int(incY) < 0 { 84 ky = uintptr(-int(n-1) * int(incY)) 85 } 86 switch { 87 case beta == 0: // beta == 0 is special-cased to memclear 88 if incY == 1 { 89 for i := range y { 90 y[i] = 0 91 } 92 } else { 93 iy := ky 94 for i := 0; i < int(n); i++ { 95 y[iy] = 0 96 iy += incY 97 } 98 } 99 case int(incY) < 0: 100 ScalInc(beta, y, n, uintptr(int(-incY))) 101 case incY == 1: 102 ScalUnitary(beta, y[:n]) 103 default: 104 ScalInc(beta, y, n, incY) 105 } 106 107 if incX == 1 && incY == 1 { 108 for i = 0; i < m; i++ { 109 AxpyUnitaryTo(y, alpha*x[i], a[lda*i:lda*i+n], y) 110 } 111 return 112 } 113 ix := kx 114 for i = 0; i < m; i++ { 115 AxpyInc(alpha*x[ix], a[lda*i:lda*i+n], y, n, 1, incY, 0, ky) 116 ix += incX 117 } 118} 119