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
|
// 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.
//go:build !amd64 || noasm || gccgo || safe
// +build !amd64 noasm gccgo safe
package f64
// Ger performs the rank-one operation
//
// A += alpha * x * yᵀ
//
// where A is an m×n dense matrix, x and y are vectors, and alpha is a scalar.
func Ger(m, n uintptr, alpha float64, x []float64, incX uintptr, y []float64, incY uintptr, a []float64, lda uintptr) {
if incX == 1 && incY == 1 {
x = x[:m]
y = y[:n]
for i, xv := range x {
AxpyUnitary(alpha*xv, y, a[uintptr(i)*lda:uintptr(i)*lda+n])
}
return
}
var ky, kx uintptr
if int(incY) < 0 {
ky = uintptr(-int(n-1) * int(incY))
}
if int(incX) < 0 {
kx = uintptr(-int(m-1) * int(incX))
}
ix := kx
for i := 0; i < int(m); i++ {
AxpyInc(alpha*x[ix], y, a[uintptr(i)*lda:uintptr(i)*lda+n], n, incY, 1, ky, 0)
ix += incX
}
}
// GemvN computes
//
// y = alpha * A * x + beta * y
//
// where A is an m×n dense matrix, x and y are vectors, and alpha and beta are scalars.
func GemvN(m, n uintptr, alpha float64, a []float64, lda uintptr, x []float64, incX uintptr, beta float64, y []float64, incY uintptr) {
var kx, ky, i uintptr
if int(incX) < 0 {
kx = uintptr(-int(n-1) * int(incX))
}
if int(incY) < 0 {
ky = uintptr(-int(m-1) * int(incY))
}
if incX == 1 && incY == 1 {
if beta == 0 {
for i = 0; i < m; i++ {
y[i] = alpha * DotUnitary(a[lda*i:lda*i+n], x)
}
return
}
for i = 0; i < m; i++ {
y[i] = y[i]*beta + alpha*DotUnitary(a[lda*i:lda*i+n], x)
}
return
}
iy := ky
if beta == 0 {
for i = 0; i < m; i++ {
y[iy] = alpha * DotInc(x, a[lda*i:lda*i+n], n, incX, 1, kx, 0)
iy += incY
}
return
}
for i = 0; i < m; i++ {
y[iy] = y[iy]*beta + alpha*DotInc(x, a[lda*i:lda*i+n], n, incX, 1, kx, 0)
iy += incY
}
}
// GemvT computes
//
// y = alpha * Aᵀ * x + beta * y
//
// where A is an m×n dense matrix, x and y are vectors, and alpha and beta are scalars.
func GemvT(m, n uintptr, alpha float64, a []float64, lda uintptr, x []float64, incX uintptr, beta float64, y []float64, incY uintptr) {
var kx, ky, i uintptr
if int(incX) < 0 {
kx = uintptr(-int(m-1) * int(incX))
}
if int(incY) < 0 {
ky = uintptr(-int(n-1) * int(incY))
}
switch {
case beta == 0: // beta == 0 is special-cased to memclear
if incY == 1 {
for i := range y {
y[i] = 0
}
} else {
iy := ky
for i := 0; i < int(n); i++ {
y[iy] = 0
iy += incY
}
}
case int(incY) < 0:
ScalInc(beta, y, n, uintptr(int(-incY)))
case incY == 1:
ScalUnitary(beta, y[:n])
default:
ScalInc(beta, y, n, incY)
}
if incX == 1 && incY == 1 {
for i = 0; i < m; i++ {
AxpyUnitaryTo(y, alpha*x[i], a[lda*i:lda*i+n], y)
}
return
}
ix := kx
for i = 0; i < m; i++ {
AxpyInc(alpha*x[ix], a[lda*i:lda*i+n], y, n, 1, incY, 0, ky)
ix += incX
}
}
|
