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