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