github.com/gopherd/gonum@v0.0.4/internal/asm/f32/gemv.go (about)

     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  package f32
     6  
     7  // GemvN computes
     8  //  y = alpha * A * x + beta * y
     9  // where A is an m×n dense matrix, x and y are vectors, and alpha and beta are scalars.
    10  func GemvN(m, n uintptr, alpha float32, a []float32, lda uintptr, x []float32, incX uintptr, beta float32, y []float32, incY uintptr) {
    11  	var kx, ky, i uintptr
    12  	if int(incX) < 0 {
    13  		kx = uintptr(-int(n-1) * int(incX))
    14  	}
    15  	if int(incY) < 0 {
    16  		ky = uintptr(-int(m-1) * int(incY))
    17  	}
    18  
    19  	if incX == 1 && incY == 1 {
    20  		if beta == 0 {
    21  			for i = 0; i < m; i++ {
    22  				y[i] = alpha * DotUnitary(a[lda*i:lda*i+n], x)
    23  			}
    24  			return
    25  		}
    26  		for i = 0; i < m; i++ {
    27  			y[i] = y[i]*beta + alpha*DotUnitary(a[lda*i:lda*i+n], x)
    28  		}
    29  		return
    30  	}
    31  	iy := ky
    32  	if beta == 0 {
    33  		for i = 0; i < m; i++ {
    34  			y[iy] = alpha * DotInc(x, a[lda*i:lda*i+n], n, incX, 1, kx, 0)
    35  			iy += incY
    36  		}
    37  		return
    38  	}
    39  	for i = 0; i < m; i++ {
    40  		y[iy] = y[iy]*beta + alpha*DotInc(x, a[lda*i:lda*i+n], n, incX, 1, kx, 0)
    41  		iy += incY
    42  	}
    43  }
    44  
    45  // GemvT computes
    46  //  y = alpha * Aᵀ * x + beta * y
    47  // where A is an m×n dense matrix, x and y are vectors, and alpha and beta are scalars.
    48  func GemvT(m, n uintptr, alpha float32, a []float32, lda uintptr, x []float32, incX uintptr, beta float32, y []float32, incY uintptr) {
    49  	var kx, ky, i uintptr
    50  	if int(incX) < 0 {
    51  		kx = uintptr(-int(m-1) * int(incX))
    52  	}
    53  	if int(incY) < 0 {
    54  		ky = uintptr(-int(n-1) * int(incY))
    55  	}
    56  	switch {
    57  	case beta == 0: // beta == 0 is special-cased to memclear
    58  		if incY == 1 {
    59  			for i := range y {
    60  				y[i] = 0
    61  			}
    62  		} else {
    63  			iy := ky
    64  			for i := 0; i < int(n); i++ {
    65  				y[iy] = 0
    66  				iy += incY
    67  			}
    68  		}
    69  	case int(incY) < 0:
    70  		ScalInc(beta, y, n, uintptr(int(-incY)))
    71  	case incY == 1:
    72  		ScalUnitary(beta, y[:n])
    73  	default:
    74  		ScalInc(beta, y, n, incY)
    75  	}
    76  
    77  	if incX == 1 && incY == 1 {
    78  		for i = 0; i < m; i++ {
    79  			AxpyUnitaryTo(y, alpha*x[i], a[lda*i:lda*i+n], y)
    80  		}
    81  		return
    82  	}
    83  	ix := kx
    84  	for i = 0; i < m; i++ {
    85  		AxpyInc(alpha*x[ix], a[lda*i:lda*i+n], y, n, 1, incY, 0, ky)
    86  		ix += incX
    87  	}
    88  }