| // Code generated by "go generate gonum.org/v1/gonum/blas/gonum”; DO NOT EDIT. |
| |
| // 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. |
| |
| package gonum |
| |
| import ( |
| cmplx "gonum.org/v1/gonum/internal/cmplx64" |
| |
| "gonum.org/v1/gonum/blas" |
| "gonum.org/v1/gonum/internal/asm/c64" |
| ) |
| |
| var _ blas.Complex64Level2 = Implementation{} |
| |
| // Cgbmv performs one of the matrix-vector operations |
| // y = alpha * A * x + beta * y if trans = blas.NoTrans |
| // y = alpha * Aᵀ * x + beta * y if trans = blas.Trans |
| // y = alpha * Aᴴ * x + beta * y if trans = blas.ConjTrans |
| // where alpha and beta are scalars, x and y are vectors, and A is an m×n band matrix |
| // with kL sub-diagonals and kU super-diagonals. |
| // |
| // Complex64 implementations are autogenerated and not directly tested. |
| func (Implementation) Cgbmv(trans blas.Transpose, m, n, kL, kU int, alpha complex64, a []complex64, lda int, x []complex64, incX int, beta complex64, y []complex64, incY int) { |
| switch trans { |
| default: |
| panic(badTranspose) |
| case blas.NoTrans, blas.Trans, blas.ConjTrans: |
| } |
| if m < 0 { |
| panic(mLT0) |
| } |
| if n < 0 { |
| panic(nLT0) |
| } |
| if kL < 0 { |
| panic(kLLT0) |
| } |
| if kU < 0 { |
| panic(kULT0) |
| } |
| if lda < kL+kU+1 { |
| panic(badLdA) |
| } |
| if incX == 0 { |
| panic(zeroIncX) |
| } |
| if incY == 0 { |
| panic(zeroIncY) |
| } |
| |
| // Quick return if possible. |
| if m == 0 || n == 0 { |
| return |
| } |
| |
| // For zero matrix size the following slice length checks are trivially satisfied. |
| if len(a) < lda*(min(m, n+kL)-1)+kL+kU+1 { |
| panic(shortA) |
| } |
| var lenX, lenY int |
| if trans == blas.NoTrans { |
| lenX, lenY = n, m |
| } else { |
| lenX, lenY = m, n |
| } |
| if (incX > 0 && len(x) <= (lenX-1)*incX) || (incX < 0 && len(x) <= (1-lenX)*incX) { |
| panic(shortX) |
| } |
| if (incY > 0 && len(y) <= (lenY-1)*incY) || (incY < 0 && len(y) <= (1-lenY)*incY) { |
| panic(shortY) |
| } |
| |
| // Quick return if possible. |
| if alpha == 0 && beta == 1 { |
| return |
| } |
| |
| var kx int |
| if incX < 0 { |
| kx = (1 - lenX) * incX |
| } |
| var ky int |
| if incY < 0 { |
| ky = (1 - lenY) * incY |
| } |
| |
| // Form y = beta*y. |
| if beta != 1 { |
| if incY == 1 { |
| if beta == 0 { |
| for i := range y[:lenY] { |
| y[i] = 0 |
| } |
| } else { |
| c64.ScalUnitary(beta, y[:lenY]) |
| } |
| } else { |
| iy := ky |
| if beta == 0 { |
| for i := 0; i < lenY; i++ { |
| y[iy] = 0 |
| iy += incY |
| } |
| } else { |
| if incY > 0 { |
| c64.ScalInc(beta, y, uintptr(lenY), uintptr(incY)) |
| } else { |
| c64.ScalInc(beta, y, uintptr(lenY), uintptr(-incY)) |
| } |
| } |
| } |
| } |
| |
| nRow := min(m, n+kL) |
| nCol := kL + 1 + kU |
| switch trans { |
| case blas.NoTrans: |
| iy := ky |
| if incX == 1 { |
| for i := 0; i < nRow; i++ { |
| l := max(0, kL-i) |
| u := min(nCol, n+kL-i) |
| aRow := a[i*lda+l : i*lda+u] |
| off := max(0, i-kL) |
| xtmp := x[off : off+u-l] |
| var sum complex64 |
| for j, v := range aRow { |
| sum += xtmp[j] * v |
| } |
| y[iy] += alpha * sum |
| iy += incY |
| } |
| } else { |
| for i := 0; i < nRow; i++ { |
| l := max(0, kL-i) |
| u := min(nCol, n+kL-i) |
| aRow := a[i*lda+l : i*lda+u] |
| off := max(0, i-kL) * incX |
| jx := kx |
| var sum complex64 |
| for _, v := range aRow { |
| sum += x[off+jx] * v |
| jx += incX |
| } |
| y[iy] += alpha * sum |
| iy += incY |
| } |
| } |
| case blas.Trans: |
| if incX == 1 { |
| for i := 0; i < nRow; i++ { |
| l := max(0, kL-i) |
| u := min(nCol, n+kL-i) |
| aRow := a[i*lda+l : i*lda+u] |
| off := max(0, i-kL) * incY |
| alphaxi := alpha * x[i] |
| jy := ky |
| for _, v := range aRow { |
| y[off+jy] += alphaxi * v |
| jy += incY |
| } |
| } |
| } else { |
| ix := kx |
| for i := 0; i < nRow; i++ { |
| l := max(0, kL-i) |
| u := min(nCol, n+kL-i) |
| aRow := a[i*lda+l : i*lda+u] |
| off := max(0, i-kL) * incY |
| alphaxi := alpha * x[ix] |
| jy := ky |
| for _, v := range aRow { |
| y[off+jy] += alphaxi * v |
| jy += incY |
| } |
| ix += incX |
| } |
| } |
| case blas.ConjTrans: |
| if incX == 1 { |
| for i := 0; i < nRow; i++ { |
| l := max(0, kL-i) |
| u := min(nCol, n+kL-i) |
| aRow := a[i*lda+l : i*lda+u] |
| off := max(0, i-kL) * incY |
| alphaxi := alpha * x[i] |
| jy := ky |
| for _, v := range aRow { |
| y[off+jy] += alphaxi * cmplx.Conj(v) |
| jy += incY |
| } |
| } |
| } else { |
| ix := kx |
| for i := 0; i < nRow; i++ { |
| l := max(0, kL-i) |
| u := min(nCol, n+kL-i) |
| aRow := a[i*lda+l : i*lda+u] |
| off := max(0, i-kL) * incY |
| alphaxi := alpha * x[ix] |
| jy := ky |
| for _, v := range aRow { |
| y[off+jy] += alphaxi * cmplx.Conj(v) |
| jy += incY |
| } |
| ix += incX |
| } |
| } |
| } |
| } |
| |
| // Cgemv performs one of the matrix-vector operations |
| // y = alpha * A * x + beta * y if trans = blas.NoTrans |
| // y = alpha * Aᵀ * x + beta * y if trans = blas.Trans |
| // y = alpha * Aᴴ * x + beta * y if trans = blas.ConjTrans |
| // where alpha and beta are scalars, x and y are vectors, and A is an m×n dense matrix. |
| // |
| // Complex64 implementations are autogenerated and not directly tested. |
| func (Implementation) Cgemv(trans blas.Transpose, m, n int, alpha complex64, a []complex64, lda int, x []complex64, incX int, beta complex64, y []complex64, incY int) { |
| switch trans { |
| default: |
| panic(badTranspose) |
| case blas.NoTrans, blas.Trans, blas.ConjTrans: |
| } |
| if m < 0 { |
| panic(mLT0) |
| } |
| if n < 0 { |
| panic(nLT0) |
| } |
| if lda < max(1, n) { |
| panic(badLdA) |
| } |
| if incX == 0 { |
| panic(zeroIncX) |
| } |
| if incY == 0 { |
| panic(zeroIncY) |
| } |
| |
| // Quick return if possible. |
| if m == 0 || n == 0 { |
| return |
| } |
| |
| // For zero matrix size the following slice length checks are trivially satisfied. |
| var lenX, lenY int |
| if trans == blas.NoTrans { |
| lenX = n |
| lenY = m |
| } else { |
| lenX = m |
| lenY = n |
| } |
| if len(a) < lda*(m-1)+n { |
| panic(shortA) |
| } |
| if (incX > 0 && len(x) <= (lenX-1)*incX) || (incX < 0 && len(x) <= (1-lenX)*incX) { |
| panic(shortX) |
| } |
| if (incY > 0 && len(y) <= (lenY-1)*incY) || (incY < 0 && len(y) <= (1-lenY)*incY) { |
| panic(shortY) |
| } |
| |
| // Quick return if possible. |
| if alpha == 0 && beta == 1 { |
| return |
| } |
| |
| var kx int |
| if incX < 0 { |
| kx = (1 - lenX) * incX |
| } |
| var ky int |
| if incY < 0 { |
| ky = (1 - lenY) * incY |
| } |
| |
| // Form y = beta*y. |
| if beta != 1 { |
| if incY == 1 { |
| if beta == 0 { |
| for i := range y[:lenY] { |
| y[i] = 0 |
| } |
| } else { |
| c64.ScalUnitary(beta, y[:lenY]) |
| } |
| } else { |
| iy := ky |
| if beta == 0 { |
| for i := 0; i < lenY; i++ { |
| y[iy] = 0 |
| iy += incY |
| } |
| } else { |
| if incY > 0 { |
| c64.ScalInc(beta, y, uintptr(lenY), uintptr(incY)) |
| } else { |
| c64.ScalInc(beta, y, uintptr(lenY), uintptr(-incY)) |
| } |
| } |
| } |
| } |
| |
| if alpha == 0 { |
| return |
| } |
| |
| switch trans { |
| default: |
| // Form y = alpha*A*x + y. |
| iy := ky |
| if incX == 1 { |
| for i := 0; i < m; i++ { |
| y[iy] += alpha * c64.DotuUnitary(a[i*lda:i*lda+n], x[:n]) |
| iy += incY |
| } |
| return |
| } |
| for i := 0; i < m; i++ { |
| y[iy] += alpha * c64.DotuInc(a[i*lda:i*lda+n], x, uintptr(n), 1, uintptr(incX), 0, uintptr(kx)) |
| iy += incY |
| } |
| return |
| |
| case blas.Trans: |
| // Form y = alpha*Aᵀ*x + y. |
| ix := kx |
| if incY == 1 { |
| for i := 0; i < m; i++ { |
| c64.AxpyUnitary(alpha*x[ix], a[i*lda:i*lda+n], y[:n]) |
| ix += incX |
| } |
| return |
| } |
| for i := 0; i < m; i++ { |
| c64.AxpyInc(alpha*x[ix], a[i*lda:i*lda+n], y, uintptr(n), 1, uintptr(incY), 0, uintptr(ky)) |
| ix += incX |
| } |
| return |
| |
| case blas.ConjTrans: |
| // Form y = alpha*Aᴴ*x + y. |
| ix := kx |
| if incY == 1 { |
| for i := 0; i < m; i++ { |
| tmp := alpha * x[ix] |
| for j := 0; j < n; j++ { |
| y[j] += tmp * cmplx.Conj(a[i*lda+j]) |
| } |
| ix += incX |
| } |
| return |
| } |
| for i := 0; i < m; i++ { |
| tmp := alpha * x[ix] |
| jy := ky |
| for j := 0; j < n; j++ { |
| y[jy] += tmp * cmplx.Conj(a[i*lda+j]) |
| jy += incY |
| } |
| ix += incX |
| } |
| return |
| } |
| } |
| |
| // Cgerc performs the rank-one operation |
| // A += alpha * x * yᴴ |
| // where A is an m×n dense matrix, alpha is a scalar, x is an m element vector, |
| // and y is an n element vector. |
| // |
| // Complex64 implementations are autogenerated and not directly tested. |
| func (Implementation) Cgerc(m, n int, alpha complex64, x []complex64, incX int, y []complex64, incY int, a []complex64, lda int) { |
| if m < 0 { |
| panic(mLT0) |
| } |
| if n < 0 { |
| panic(nLT0) |
| } |
| if lda < max(1, n) { |
| panic(badLdA) |
| } |
| if incX == 0 { |
| panic(zeroIncX) |
| } |
| if incY == 0 { |
| panic(zeroIncY) |
| } |
| |
| // Quick return if possible. |
| if m == 0 || n == 0 { |
| return |
| } |
| |
| // For zero matrix size the following slice length checks are trivially satisfied. |
| if (incX > 0 && len(x) <= (m-1)*incX) || (incX < 0 && len(x) <= (1-m)*incX) { |
| panic(shortX) |
| } |
| if (incY > 0 && len(y) <= (n-1)*incY) || (incY < 0 && len(y) <= (1-n)*incY) { |
| panic(shortY) |
| } |
| if len(a) < lda*(m-1)+n { |
| panic(shortA) |
| } |
| |
| // Quick return if possible. |
| if alpha == 0 { |
| return |
| } |
| |
| var kx, jy int |
| if incX < 0 { |
| kx = (1 - m) * incX |
| } |
| if incY < 0 { |
| jy = (1 - n) * incY |
| } |
| for j := 0; j < n; j++ { |
| if y[jy] != 0 { |
| tmp := alpha * cmplx.Conj(y[jy]) |
| c64.AxpyInc(tmp, x, a[j:], uintptr(m), uintptr(incX), uintptr(lda), uintptr(kx), 0) |
| } |
| jy += incY |
| } |
| } |
| |
| // Cgeru performs the rank-one operation |
| // A += alpha * x * yᵀ |
| // where A is an m×n dense matrix, alpha is a scalar, x is an m element vector, |
| // and y is an n element vector. |
| // |
| // Complex64 implementations are autogenerated and not directly tested. |
| func (Implementation) Cgeru(m, n int, alpha complex64, x []complex64, incX int, y []complex64, incY int, a []complex64, lda int) { |
| if m < 0 { |
| panic(mLT0) |
| } |
| if n < 0 { |
| panic(nLT0) |
| } |
| if lda < max(1, n) { |
| panic(badLdA) |
| } |
| if incX == 0 { |
| panic(zeroIncX) |
| } |
| if incY == 0 { |
| panic(zeroIncY) |
| } |
| |
| // Quick return if possible. |
| if m == 0 || n == 0 { |
| return |
| } |
| |
| // For zero matrix size the following slice length checks are trivially satisfied. |
| if (incX > 0 && len(x) <= (m-1)*incX) || (incX < 0 && len(x) <= (1-m)*incX) { |
| panic(shortX) |
| } |
| if (incY > 0 && len(y) <= (n-1)*incY) || (incY < 0 && len(y) <= (1-n)*incY) { |
| panic(shortY) |
| } |
| if len(a) < lda*(m-1)+n { |
| panic(shortA) |
| } |
| |
| // Quick return if possible. |
| if alpha == 0 { |
| return |
| } |
| |
| var kx int |
| if incX < 0 { |
| kx = (1 - m) * incX |
| } |
| if incY == 1 { |
| for i := 0; i < m; i++ { |
| if x[kx] != 0 { |
| tmp := alpha * x[kx] |
| c64.AxpyUnitary(tmp, y[:n], a[i*lda:i*lda+n]) |
| } |
| kx += incX |
| } |
| return |
| } |
| var jy int |
| if incY < 0 { |
| jy = (1 - n) * incY |
| } |
| for i := 0; i < m; i++ { |
| if x[kx] != 0 { |
| tmp := alpha * x[kx] |
| c64.AxpyInc(tmp, y, a[i*lda:i*lda+n], uintptr(n), uintptr(incY), 1, uintptr(jy), 0) |
| } |
| kx += incX |
| } |
| } |
| |
| // Chbmv performs the matrix-vector operation |
| // y = alpha * A * x + beta * y |
| // where alpha and beta are scalars, x and y are vectors, and A is an n×n |
| // Hermitian band matrix with k super-diagonals. The imaginary parts of |
| // the diagonal elements of A are ignored and assumed to be zero. |
| // |
| // Complex64 implementations are autogenerated and not directly tested. |
| func (Implementation) Chbmv(uplo blas.Uplo, n, k int, alpha complex64, a []complex64, lda int, x []complex64, incX int, beta complex64, y []complex64, incY int) { |
| switch uplo { |
| default: |
| panic(badUplo) |
| case blas.Upper, blas.Lower: |
| } |
| if n < 0 { |
| panic(nLT0) |
| } |
| if k < 0 { |
| panic(kLT0) |
| } |
| if lda < k+1 { |
| panic(badLdA) |
| } |
| if incX == 0 { |
| panic(zeroIncX) |
| } |
| if incY == 0 { |
| panic(zeroIncY) |
| } |
| |
| // Quick return if possible. |
| if n == 0 { |
| return |
| } |
| |
| // For zero matrix size the following slice length checks are trivially satisfied. |
| if len(a) < lda*(n-1)+k+1 { |
| panic(shortA) |
| } |
| if (incX > 0 && len(x) <= (n-1)*incX) || (incX < 0 && len(x) <= (1-n)*incX) { |
| panic(shortX) |
| } |
| if (incY > 0 && len(y) <= (n-1)*incY) || (incY < 0 && len(y) <= (1-n)*incY) { |
| panic(shortY) |
| } |
| |
| // Quick return if possible. |
| if alpha == 0 && beta == 1 { |
| return |
| } |
| |
| // Set up the start indices in X and Y. |
| var kx int |
| if incX < 0 { |
| kx = (1 - n) * incX |
| } |
| var ky int |
| if incY < 0 { |
| ky = (1 - n) * incY |
| } |
| |
| // Form y = beta*y. |
| if beta != 1 { |
| if incY == 1 { |
| if beta == 0 { |
| for i := range y[:n] { |
| y[i] = 0 |
| } |
| } else { |
| for i, v := range y[:n] { |
| y[i] = beta * v |
| } |
| } |
| } else { |
| iy := ky |
| if beta == 0 { |
| for i := 0; i < n; i++ { |
| y[iy] = 0 |
| iy += incY |
| } |
| } else { |
| for i := 0; i < n; i++ { |
| y[iy] = beta * y[iy] |
| iy += incY |
| } |
| } |
| } |
| } |
| |
| if alpha == 0 { |
| return |
| } |
| |
| // The elements of A are accessed sequentially with one pass through a. |
| switch uplo { |
| case blas.Upper: |
| iy := ky |
| if incX == 1 { |
| for i := 0; i < n; i++ { |
| aRow := a[i*lda:] |
| alphaxi := alpha * x[i] |
| sum := alphaxi * complex(real(aRow[0]), 0) |
| u := min(k+1, n-i) |
| jy := incY |
| for j := 1; j < u; j++ { |
| v := aRow[j] |
| sum += alpha * x[i+j] * v |
| y[iy+jy] += alphaxi * cmplx.Conj(v) |
| jy += incY |
| } |
| y[iy] += sum |
| iy += incY |
| } |
| } else { |
| ix := kx |
| for i := 0; i < n; i++ { |
| aRow := a[i*lda:] |
| alphaxi := alpha * x[ix] |
| sum := alphaxi * complex(real(aRow[0]), 0) |
| u := min(k+1, n-i) |
| jx := incX |
| jy := incY |
| for j := 1; j < u; j++ { |
| v := aRow[j] |
| sum += alpha * x[ix+jx] * v |
| y[iy+jy] += alphaxi * cmplx.Conj(v) |
| jx += incX |
| jy += incY |
| } |
| y[iy] += sum |
| ix += incX |
| iy += incY |
| } |
| } |
| case blas.Lower: |
| iy := ky |
| if incX == 1 { |
| for i := 0; i < n; i++ { |
| l := max(0, k-i) |
| alphaxi := alpha * x[i] |
| jy := l * incY |
| aRow := a[i*lda:] |
| for j := l; j < k; j++ { |
| v := aRow[j] |
| y[iy] += alpha * v * x[i-k+j] |
| y[iy-k*incY+jy] += alphaxi * cmplx.Conj(v) |
| jy += incY |
| } |
| y[iy] += alphaxi * complex(real(aRow[k]), 0) |
| iy += incY |
| } |
| } else { |
| ix := kx |
| for i := 0; i < n; i++ { |
| l := max(0, k-i) |
| alphaxi := alpha * x[ix] |
| jx := l * incX |
| jy := l * incY |
| aRow := a[i*lda:] |
| for j := l; j < k; j++ { |
| v := aRow[j] |
| y[iy] += alpha * v * x[ix-k*incX+jx] |
| y[iy-k*incY+jy] += alphaxi * cmplx.Conj(v) |
| jx += incX |
| jy += incY |
| } |
| y[iy] += alphaxi * complex(real(aRow[k]), 0) |
| ix += incX |
| iy += incY |
| } |
| } |
| } |
| } |
| |
| // Chemv performs the matrix-vector operation |
| // y = alpha * A * x + beta * y |
| // where alpha and beta are scalars, x and y are vectors, and A is an n×n |
| // Hermitian matrix. The imaginary parts of the diagonal elements of A are |
| // ignored and assumed to be zero. |
| // |
| // Complex64 implementations are autogenerated and not directly tested. |
| func (Implementation) Chemv(uplo blas.Uplo, n int, alpha complex64, a []complex64, lda int, x []complex64, incX int, beta complex64, y []complex64, incY int) { |
| switch uplo { |
| default: |
| panic(badUplo) |
| case blas.Upper, blas.Lower: |
| } |
| if n < 0 { |
| panic(nLT0) |
| } |
| if lda < max(1, n) { |
| panic(badLdA) |
| } |
| if incX == 0 { |
| panic(zeroIncX) |
| } |
| if incY == 0 { |
| panic(zeroIncY) |
| } |
| |
| // Quick return if possible. |
| if n == 0 { |
| return |
| } |
| |
| // For zero matrix size the following slice length checks are trivially satisfied. |
| if len(a) < lda*(n-1)+n { |
| panic(shortA) |
| } |
| if (incX > 0 && len(x) <= (n-1)*incX) || (incX < 0 && len(x) <= (1-n)*incX) { |
| panic(shortX) |
| } |
| if (incY > 0 && len(y) <= (n-1)*incY) || (incY < 0 && len(y) <= (1-n)*incY) { |
| panic(shortY) |
| } |
| |
| // Quick return if possible. |
| if alpha == 0 && beta == 1 { |
| return |
| } |
| |
| // Set up the start indices in X and Y. |
| var kx int |
| if incX < 0 { |
| kx = (1 - n) * incX |
| } |
| var ky int |
| if incY < 0 { |
| ky = (1 - n) * incY |
| } |
| |
| // Form y = beta*y. |
| if beta != 1 { |
| if incY == 1 { |
| if beta == 0 { |
| for i := range y[:n] { |
| y[i] = 0 |
| } |
| } else { |
| for i, v := range y[:n] { |
| y[i] = beta * v |
| } |
| } |
| } else { |
| iy := ky |
| if beta == 0 { |
| for i := 0; i < n; i++ { |
| y[iy] = 0 |
| iy += incY |
| } |
| } else { |
| for i := 0; i < n; i++ { |
| y[iy] = beta * y[iy] |
| iy += incY |
| } |
| } |
| } |
| } |
| |
| if alpha == 0 { |
| return |
| } |
| |
| // The elements of A are accessed sequentially with one pass through |
| // the triangular part of A. |
| |
| if uplo == blas.Upper { |
| // Form y when A is stored in upper triangle. |
| if incX == 1 && incY == 1 { |
| for i := 0; i < n; i++ { |
| tmp1 := alpha * x[i] |
| var tmp2 complex64 |
| for j := i + 1; j < n; j++ { |
| y[j] += tmp1 * cmplx.Conj(a[i*lda+j]) |
| tmp2 += a[i*lda+j] * x[j] |
| } |
| aii := complex(real(a[i*lda+i]), 0) |
| y[i] += tmp1*aii + alpha*tmp2 |
| } |
| } else { |
| ix := kx |
| iy := ky |
| for i := 0; i < n; i++ { |
| tmp1 := alpha * x[ix] |
| var tmp2 complex64 |
| jx := ix |
| jy := iy |
| for j := i + 1; j < n; j++ { |
| jx += incX |
| jy += incY |
| y[jy] += tmp1 * cmplx.Conj(a[i*lda+j]) |
| tmp2 += a[i*lda+j] * x[jx] |
| } |
| aii := complex(real(a[i*lda+i]), 0) |
| y[iy] += tmp1*aii + alpha*tmp2 |
| ix += incX |
| iy += incY |
| } |
| } |
| return |
| } |
| |
| // Form y when A is stored in lower triangle. |
| if incX == 1 && incY == 1 { |
| for i := 0; i < n; i++ { |
| tmp1 := alpha * x[i] |
| var tmp2 complex64 |
| for j := 0; j < i; j++ { |
| y[j] += tmp1 * cmplx.Conj(a[i*lda+j]) |
| tmp2 += a[i*lda+j] * x[j] |
| } |
| aii := complex(real(a[i*lda+i]), 0) |
| y[i] += tmp1*aii + alpha*tmp2 |
| } |
| } else { |
| ix := kx |
| iy := ky |
| for i := 0; i < n; i++ { |
| tmp1 := alpha * x[ix] |
| var tmp2 complex64 |
| jx := kx |
| jy := ky |
| for j := 0; j < i; j++ { |
| y[jy] += tmp1 * cmplx.Conj(a[i*lda+j]) |
| tmp2 += a[i*lda+j] * x[jx] |
| jx += incX |
| jy += incY |
| } |
| aii := complex(real(a[i*lda+i]), 0) |
| y[iy] += tmp1*aii + alpha*tmp2 |
| ix += incX |
| iy += incY |
| } |
| } |
| } |
| |
| // Cher performs the Hermitian rank-one operation |
| // A += alpha * x * xᴴ |
| // where A is an n×n Hermitian matrix, alpha is a real scalar, and x is an n |
| // element vector. On entry, the imaginary parts of the diagonal elements of A |
| // are ignored and assumed to be zero, on return they will be set to zero. |
| // |
| // Complex64 implementations are autogenerated and not directly tested. |
| func (Implementation) Cher(uplo blas.Uplo, n int, alpha float32, x []complex64, incX int, a []complex64, lda int) { |
| switch uplo { |
| default: |
| panic(badUplo) |
| case blas.Upper, blas.Lower: |
| } |
| if n < 0 { |
| panic(nLT0) |
| } |
| if lda < max(1, n) { |
| panic(badLdA) |
| } |
| if incX == 0 { |
| panic(zeroIncX) |
| } |
| |
| // Quick return if possible. |
| if n == 0 { |
| return |
| } |
| |
| // For zero matrix size the following slice length checks are trivially satisfied. |
| if (incX > 0 && len(x) <= (n-1)*incX) || (incX < 0 && len(x) <= (1-n)*incX) { |
| panic(shortX) |
| } |
| if len(a) < lda*(n-1)+n { |
| panic(shortA) |
| } |
| |
| // Quick return if possible. |
| if alpha == 0 { |
| return |
| } |
| |
| var kx int |
| if incX < 0 { |
| kx = (1 - n) * incX |
| } |
| if uplo == blas.Upper { |
| if incX == 1 { |
| for i := 0; i < n; i++ { |
| if x[i] != 0 { |
| tmp := complex(alpha*real(x[i]), alpha*imag(x[i])) |
| aii := real(a[i*lda+i]) |
| xtmp := real(tmp * cmplx.Conj(x[i])) |
| a[i*lda+i] = complex(aii+xtmp, 0) |
| for j := i + 1; j < n; j++ { |
| a[i*lda+j] += tmp * cmplx.Conj(x[j]) |
| } |
| } else { |
| aii := real(a[i*lda+i]) |
| a[i*lda+i] = complex(aii, 0) |
| } |
| } |
| return |
| } |
| |
| ix := kx |
| for i := 0; i < n; i++ { |
| if x[ix] != 0 { |
| tmp := complex(alpha*real(x[ix]), alpha*imag(x[ix])) |
| aii := real(a[i*lda+i]) |
| xtmp := real(tmp * cmplx.Conj(x[ix])) |
| a[i*lda+i] = complex(aii+xtmp, 0) |
| jx := ix + incX |
| for j := i + 1; j < n; j++ { |
| a[i*lda+j] += tmp * cmplx.Conj(x[jx]) |
| jx += incX |
| } |
| } else { |
| aii := real(a[i*lda+i]) |
| a[i*lda+i] = complex(aii, 0) |
| } |
| ix += incX |
| } |
| return |
| } |
| |
| if incX == 1 { |
| for i := 0; i < n; i++ { |
| if x[i] != 0 { |
| tmp := complex(alpha*real(x[i]), alpha*imag(x[i])) |
| for j := 0; j < i; j++ { |
| a[i*lda+j] += tmp * cmplx.Conj(x[j]) |
| } |
| aii := real(a[i*lda+i]) |
| xtmp := real(tmp * cmplx.Conj(x[i])) |
| a[i*lda+i] = complex(aii+xtmp, 0) |
| } else { |
| aii := real(a[i*lda+i]) |
| a[i*lda+i] = complex(aii, 0) |
| } |
| } |
| return |
| } |
| |
| ix := kx |
| for i := 0; i < n; i++ { |
| if x[ix] != 0 { |
| tmp := complex(alpha*real(x[ix]), alpha*imag(x[ix])) |
| jx := kx |
| for j := 0; j < i; j++ { |
| a[i*lda+j] += tmp * cmplx.Conj(x[jx]) |
| jx += incX |
| } |
| aii := real(a[i*lda+i]) |
| xtmp := real(tmp * cmplx.Conj(x[ix])) |
| a[i*lda+i] = complex(aii+xtmp, 0) |
| |
| } else { |
| aii := real(a[i*lda+i]) |
| a[i*lda+i] = complex(aii, 0) |
| } |
| ix += incX |
| } |
| } |
| |
| // Cher2 performs the Hermitian rank-two operation |
| // A += alpha * x * yᴴ + conj(alpha) * y * xᴴ |
| // where alpha is a scalar, x and y are n element vectors and A is an n×n |
| // Hermitian matrix. On entry, the imaginary parts of the diagonal elements are |
| // ignored and assumed to be zero. On return they will be set to zero. |
| // |
| // Complex64 implementations are autogenerated and not directly tested. |
| func (Implementation) Cher2(uplo blas.Uplo, n int, alpha complex64, x []complex64, incX int, y []complex64, incY int, a []complex64, lda int) { |
| switch uplo { |
| default: |
| panic(badUplo) |
| case blas.Upper, blas.Lower: |
| } |
| if n < 0 { |
| panic(nLT0) |
| } |
| if lda < max(1, n) { |
| panic(badLdA) |
| } |
| if incX == 0 { |
| panic(zeroIncX) |
| } |
| if incY == 0 { |
| panic(zeroIncY) |
| } |
| |
| // Quick return if possible. |
| if n == 0 { |
| return |
| } |
| |
| // For zero matrix size the following slice length checks are trivially satisfied. |
| if (incX > 0 && len(x) <= (n-1)*incX) || (incX < 0 && len(x) <= (1-n)*incX) { |
| panic(shortX) |
| } |
| if (incY > 0 && len(y) <= (n-1)*incY) || (incY < 0 && len(y) <= (1-n)*incY) { |
| panic(shortY) |
| } |
| if len(a) < lda*(n-1)+n { |
| panic(shortA) |
| } |
| |
| // Quick return if possible. |
| if alpha == 0 { |
| return |
| } |
| |
| var kx, ky int |
| var ix, iy int |
| if incX != 1 || incY != 1 { |
| if incX < 0 { |
| kx = (1 - n) * incX |
| } |
| if incY < 0 { |
| ky = (1 - n) * incY |
| } |
| ix = kx |
| iy = ky |
| } |
| if uplo == blas.Upper { |
| if incX == 1 && incY == 1 { |
| for i := 0; i < n; i++ { |
| if x[i] != 0 || y[i] != 0 { |
| tmp1 := alpha * x[i] |
| tmp2 := cmplx.Conj(alpha) * y[i] |
| aii := real(a[i*lda+i]) + real(tmp1*cmplx.Conj(y[i])) + real(tmp2*cmplx.Conj(x[i])) |
| a[i*lda+i] = complex(aii, 0) |
| for j := i + 1; j < n; j++ { |
| a[i*lda+j] += tmp1*cmplx.Conj(y[j]) + tmp2*cmplx.Conj(x[j]) |
| } |
| } else { |
| aii := real(a[i*lda+i]) |
| a[i*lda+i] = complex(aii, 0) |
| } |
| } |
| return |
| } |
| for i := 0; i < n; i++ { |
| if x[ix] != 0 || y[iy] != 0 { |
| tmp1 := alpha * x[ix] |
| tmp2 := cmplx.Conj(alpha) * y[iy] |
| aii := real(a[i*lda+i]) + real(tmp1*cmplx.Conj(y[iy])) + real(tmp2*cmplx.Conj(x[ix])) |
| a[i*lda+i] = complex(aii, 0) |
| jx := ix + incX |
| jy := iy + incY |
| for j := i + 1; j < n; j++ { |
| a[i*lda+j] += tmp1*cmplx.Conj(y[jy]) + tmp2*cmplx.Conj(x[jx]) |
| jx += incX |
| jy += incY |
| } |
| } else { |
| aii := real(a[i*lda+i]) |
| a[i*lda+i] = complex(aii, 0) |
| } |
| ix += incX |
| iy += incY |
| } |
| return |
| } |
| |
| if incX == 1 && incY == 1 { |
| for i := 0; i < n; i++ { |
| if x[i] != 0 || y[i] != 0 { |
| tmp1 := alpha * x[i] |
| tmp2 := cmplx.Conj(alpha) * y[i] |
| for j := 0; j < i; j++ { |
| a[i*lda+j] += tmp1*cmplx.Conj(y[j]) + tmp2*cmplx.Conj(x[j]) |
| } |
| aii := real(a[i*lda+i]) + real(tmp1*cmplx.Conj(y[i])) + real(tmp2*cmplx.Conj(x[i])) |
| a[i*lda+i] = complex(aii, 0) |
| } else { |
| aii := real(a[i*lda+i]) |
| a[i*lda+i] = complex(aii, 0) |
| } |
| } |
| return |
| } |
| for i := 0; i < n; i++ { |
| if x[ix] != 0 || y[iy] != 0 { |
| tmp1 := alpha * x[ix] |
| tmp2 := cmplx.Conj(alpha) * y[iy] |
| jx := kx |
| jy := ky |
| for j := 0; j < i; j++ { |
| a[i*lda+j] += tmp1*cmplx.Conj(y[jy]) + tmp2*cmplx.Conj(x[jx]) |
| jx += incX |
| jy += incY |
| } |
| aii := real(a[i*lda+i]) + real(tmp1*cmplx.Conj(y[iy])) + real(tmp2*cmplx.Conj(x[ix])) |
| a[i*lda+i] = complex(aii, 0) |
| } else { |
| aii := real(a[i*lda+i]) |
| a[i*lda+i] = complex(aii, 0) |
| } |
| ix += incX |
| iy += incY |
| } |
| } |
| |
| // Chpmv performs the matrix-vector operation |
| // y = alpha * A * x + beta * y |
| // where alpha and beta are scalars, x and y are vectors, and A is an n×n |
| // Hermitian matrix in packed form. The imaginary parts of the diagonal |
| // elements of A are ignored and assumed to be zero. |
| // |
| // Complex64 implementations are autogenerated and not directly tested. |
| func (Implementation) Chpmv(uplo blas.Uplo, n int, alpha complex64, ap []complex64, x []complex64, incX int, beta complex64, y []complex64, incY int) { |
| switch uplo { |
| default: |
| panic(badUplo) |
| case blas.Upper, blas.Lower: |
| } |
| if n < 0 { |
| panic(nLT0) |
| } |
| if incX == 0 { |
| panic(zeroIncX) |
| } |
| if incY == 0 { |
| panic(zeroIncY) |
| } |
| |
| // Quick return if possible. |
| if n == 0 { |
| return |
| } |
| |
| // For zero matrix size the following slice length checks are trivially satisfied. |
| if len(ap) < n*(n+1)/2 { |
| panic(shortAP) |
| } |
| if (incX > 0 && len(x) <= (n-1)*incX) || (incX < 0 && len(x) <= (1-n)*incX) { |
| panic(shortX) |
| } |
| if (incY > 0 && len(y) <= (n-1)*incY) || (incY < 0 && len(y) <= (1-n)*incY) { |
| panic(shortY) |
| } |
| |
| // Quick return if possible. |
| if alpha == 0 && beta == 1 { |
| return |
| } |
| |
| // Set up the start indices in X and Y. |
| var kx int |
| if incX < 0 { |
| kx = (1 - n) * incX |
| } |
| var ky int |
| if incY < 0 { |
| ky = (1 - n) * incY |
| } |
| |
| // Form y = beta*y. |
| if beta != 1 { |
| if incY == 1 { |
| if beta == 0 { |
| for i := range y[:n] { |
| y[i] = 0 |
| } |
| } else { |
| for i, v := range y[:n] { |
| y[i] = beta * v |
| } |
| } |
| } else { |
| iy := ky |
| if beta == 0 { |
| for i := 0; i < n; i++ { |
| y[iy] = 0 |
| iy += incY |
| } |
| } else { |
| for i := 0; i < n; i++ { |
| y[iy] *= beta |
| iy += incY |
| } |
| } |
| } |
| } |
| |
| if alpha == 0 { |
| return |
| } |
| |
| // The elements of A are accessed sequentially with one pass through ap. |
| |
| var kk int |
| if uplo == blas.Upper { |
| // Form y when ap contains the upper triangle. |
| // Here, kk points to the current diagonal element in ap. |
| if incX == 1 && incY == 1 { |
| for i := 0; i < n; i++ { |
| tmp1 := alpha * x[i] |
| y[i] += tmp1 * complex(real(ap[kk]), 0) |
| var tmp2 complex64 |
| k := kk + 1 |
| for j := i + 1; j < n; j++ { |
| y[j] += tmp1 * cmplx.Conj(ap[k]) |
| tmp2 += ap[k] * x[j] |
| k++ |
| } |
| y[i] += alpha * tmp2 |
| kk += n - i |
| } |
| } else { |
| ix := kx |
| iy := ky |
| for i := 0; i < n; i++ { |
| tmp1 := alpha * x[ix] |
| y[iy] += tmp1 * complex(real(ap[kk]), 0) |
| var tmp2 complex64 |
| jx := ix |
| jy := iy |
| for k := kk + 1; k < kk+n-i; k++ { |
| jx += incX |
| jy += incY |
| y[jy] += tmp1 * cmplx.Conj(ap[k]) |
| tmp2 += ap[k] * x[jx] |
| } |
| y[iy] += alpha * tmp2 |
| ix += incX |
| iy += incY |
| kk += n - i |
| } |
| } |
| return |
| } |
| |
| // Form y when ap contains the lower triangle. |
| // Here, kk points to the beginning of current row in ap. |
| if incX == 1 && incY == 1 { |
| for i := 0; i < n; i++ { |
| tmp1 := alpha * x[i] |
| var tmp2 complex64 |
| k := kk |
| for j := 0; j < i; j++ { |
| y[j] += tmp1 * cmplx.Conj(ap[k]) |
| tmp2 += ap[k] * x[j] |
| k++ |
| } |
| aii := complex(real(ap[kk+i]), 0) |
| y[i] += tmp1*aii + alpha*tmp2 |
| kk += i + 1 |
| } |
| } else { |
| ix := kx |
| iy := ky |
| for i := 0; i < n; i++ { |
| tmp1 := alpha * x[ix] |
| var tmp2 complex64 |
| jx := kx |
| jy := ky |
| for k := kk; k < kk+i; k++ { |
| y[jy] += tmp1 * cmplx.Conj(ap[k]) |
| tmp2 += ap[k] * x[jx] |
| jx += incX |
| jy += incY |
| } |
| aii := complex(real(ap[kk+i]), 0) |
| y[iy] += tmp1*aii + alpha*tmp2 |
| ix += incX |
| iy += incY |
| kk += i + 1 |
| } |
| } |
| } |
| |
| // Chpr performs the Hermitian rank-1 operation |
| // A += alpha * x * xᴴ |
| // where alpha is a real scalar, x is a vector, and A is an n×n hermitian matrix |
| // in packed form. On entry, the imaginary parts of the diagonal elements are |
| // assumed to be zero, and on return they are set to zero. |
| // |
| // Complex64 implementations are autogenerated and not directly tested. |
| func (Implementation) Chpr(uplo blas.Uplo, n int, alpha float32, x []complex64, incX int, ap []complex64) { |
| switch uplo { |
| default: |
| panic(badUplo) |
| case blas.Upper, blas.Lower: |
| } |
| if n < 0 { |
| panic(nLT0) |
| } |
| if incX == 0 { |
| panic(zeroIncX) |
| } |
| |
| // Quick return if possible. |
| if n == 0 { |
| return |
| } |
| |
| // For zero matrix size the following slice length checks are trivially satisfied. |
| if (incX > 0 && len(x) <= (n-1)*incX) || (incX < 0 && len(x) <= (1-n)*incX) { |
| panic(shortX) |
| } |
| if len(ap) < n*(n+1)/2 { |
| panic(shortAP) |
| } |
| |
| // Quick return if possible. |
| if alpha == 0 { |
| return |
| } |
| |
| // Set up start index in X. |
| var kx int |
| if incX < 0 { |
| kx = (1 - n) * incX |
| } |
| |
| // The elements of A are accessed sequentially with one pass through ap. |
| |
| var kk int |
| if uplo == blas.Upper { |
| // Form A when upper triangle is stored in AP. |
| // Here, kk points to the current diagonal element in ap. |
| if incX == 1 { |
| for i := 0; i < n; i++ { |
| xi := x[i] |
| if xi != 0 { |
| aii := real(ap[kk]) + alpha*real(cmplx.Conj(xi)*xi) |
| ap[kk] = complex(aii, 0) |
| |
| tmp := complex(alpha, 0) * xi |
| a := ap[kk+1 : kk+n-i] |
| x := x[i+1 : n] |
| for j, v := range x { |
| a[j] += tmp * cmplx.Conj(v) |
| } |
| } else { |
| ap[kk] = complex(real(ap[kk]), 0) |
| } |
| kk += n - i |
| } |
| } else { |
| ix := kx |
| for i := 0; i < n; i++ { |
| xi := x[ix] |
| if xi != 0 { |
| aii := real(ap[kk]) + alpha*real(cmplx.Conj(xi)*xi) |
| ap[kk] = complex(aii, 0) |
| |
| tmp := complex(alpha, 0) * xi |
| jx := ix + incX |
| a := ap[kk+1 : kk+n-i] |
| for k := range a { |
| a[k] += tmp * cmplx.Conj(x[jx]) |
| jx += incX |
| } |
| } else { |
| ap[kk] = complex(real(ap[kk]), 0) |
| } |
| ix += incX |
| kk += n - i |
| } |
| } |
| return |
| } |
| |
| // Form A when lower triangle is stored in AP. |
| // Here, kk points to the beginning of current row in ap. |
| if incX == 1 { |
| for i := 0; i < n; i++ { |
| xi := x[i] |
| if xi != 0 { |
| tmp := complex(alpha, 0) * xi |
| a := ap[kk : kk+i] |
| for j, v := range x[:i] { |
| a[j] += tmp * cmplx.Conj(v) |
| } |
| |
| aii := real(ap[kk+i]) + alpha*real(cmplx.Conj(xi)*xi) |
| ap[kk+i] = complex(aii, 0) |
| } else { |
| ap[kk+i] = complex(real(ap[kk+i]), 0) |
| } |
| kk += i + 1 |
| } |
| } else { |
| ix := kx |
| for i := 0; i < n; i++ { |
| xi := x[ix] |
| if xi != 0 { |
| tmp := complex(alpha, 0) * xi |
| a := ap[kk : kk+i] |
| jx := kx |
| for k := range a { |
| a[k] += tmp * cmplx.Conj(x[jx]) |
| jx += incX |
| } |
| |
| aii := real(ap[kk+i]) + alpha*real(cmplx.Conj(xi)*xi) |
| ap[kk+i] = complex(aii, 0) |
| } else { |
| ap[kk+i] = complex(real(ap[kk+i]), 0) |
| } |
| ix += incX |
| kk += i + 1 |
| } |
| } |
| } |
| |
| // Chpr2 performs the Hermitian rank-2 operation |
| // A += alpha * x * yᴴ + conj(alpha) * y * xᴴ |
| // where alpha is a complex scalar, x and y are n element vectors, and A is an |
| // n×n Hermitian matrix, supplied in packed form. On entry, the imaginary parts |
| // of the diagonal elements are assumed to be zero, and on return they are set to zero. |
| // |
| // Complex64 implementations are autogenerated and not directly tested. |
| func (Implementation) Chpr2(uplo blas.Uplo, n int, alpha complex64, x []complex64, incX int, y []complex64, incY int, ap []complex64) { |
| switch uplo { |
| default: |
| panic(badUplo) |
| case blas.Upper, blas.Lower: |
| } |
| if n < 0 { |
| panic(nLT0) |
| } |
| if incX == 0 { |
| panic(zeroIncX) |
| } |
| if incY == 0 { |
| panic(zeroIncY) |
| } |
| |
| // Quick return if possible. |
| if n == 0 { |
| return |
| } |
| |
| // For zero matrix size the following slice length checks are trivially satisfied. |
| if (incX > 0 && len(x) <= (n-1)*incX) || (incX < 0 && len(x) <= (1-n)*incX) { |
| panic(shortX) |
| } |
| if (incY > 0 && len(y) <= (n-1)*incY) || (incY < 0 && len(y) <= (1-n)*incY) { |
| panic(shortY) |
| } |
| if len(ap) < n*(n+1)/2 { |
| panic(shortAP) |
| } |
| |
| // Quick return if possible. |
| if alpha == 0 { |
| return |
| } |
| |
| // Set up start indices in X and Y. |
| var kx int |
| if incX < 0 { |
| kx = (1 - n) * incX |
| } |
| var ky int |
| if incY < 0 { |
| ky = (1 - n) * incY |
| } |
| |
| // The elements of A are accessed sequentially with one pass through ap. |
| |
| var kk int |
| if uplo == blas.Upper { |
| // Form A when upper triangle is stored in AP. |
| // Here, kk points to the current diagonal element in ap. |
| if incX == 1 && incY == 1 { |
| for i := 0; i < n; i++ { |
| if x[i] != 0 || y[i] != 0 { |
| tmp1 := alpha * x[i] |
| tmp2 := cmplx.Conj(alpha) * y[i] |
| aii := real(ap[kk]) + real(tmp1*cmplx.Conj(y[i])) + real(tmp2*cmplx.Conj(x[i])) |
| ap[kk] = complex(aii, 0) |
| k := kk + 1 |
| for j := i + 1; j < n; j++ { |
| ap[k] += tmp1*cmplx.Conj(y[j]) + tmp2*cmplx.Conj(x[j]) |
| k++ |
| } |
| } else { |
| ap[kk] = complex(real(ap[kk]), 0) |
| } |
| kk += n - i |
| } |
| } else { |
| ix := kx |
| iy := ky |
| for i := 0; i < n; i++ { |
| if x[ix] != 0 || y[iy] != 0 { |
| tmp1 := alpha * x[ix] |
| tmp2 := cmplx.Conj(alpha) * y[iy] |
| aii := real(ap[kk]) + real(tmp1*cmplx.Conj(y[iy])) + real(tmp2*cmplx.Conj(x[ix])) |
| ap[kk] = complex(aii, 0) |
| jx := ix + incX |
| jy := iy + incY |
| for k := kk + 1; k < kk+n-i; k++ { |
| ap[k] += tmp1*cmplx.Conj(y[jy]) + tmp2*cmplx.Conj(x[jx]) |
| jx += incX |
| jy += incY |
| } |
| } else { |
| ap[kk] = complex(real(ap[kk]), 0) |
| } |
| ix += incX |
| iy += incY |
| kk += n - i |
| } |
| } |
| return |
| } |
| |
| // Form A when lower triangle is stored in AP. |
| // Here, kk points to the beginning of current row in ap. |
| if incX == 1 && incY == 1 { |
| for i := 0; i < n; i++ { |
| if x[i] != 0 || y[i] != 0 { |
| tmp1 := alpha * x[i] |
| tmp2 := cmplx.Conj(alpha) * y[i] |
| k := kk |
| for j := 0; j < i; j++ { |
| ap[k] += tmp1*cmplx.Conj(y[j]) + tmp2*cmplx.Conj(x[j]) |
| k++ |
| } |
| aii := real(ap[kk+i]) + real(tmp1*cmplx.Conj(y[i])) + real(tmp2*cmplx.Conj(x[i])) |
| ap[kk+i] = complex(aii, 0) |
| } else { |
| ap[kk+i] = complex(real(ap[kk+i]), 0) |
| } |
| kk += i + 1 |
| } |
| } else { |
| ix := kx |
| iy := ky |
| for i := 0; i < n; i++ { |
| if x[ix] != 0 || y[iy] != 0 { |
| tmp1 := alpha * x[ix] |
| tmp2 := cmplx.Conj(alpha) * y[iy] |
| jx := kx |
| jy := ky |
| for k := kk; k < kk+i; k++ { |
| ap[k] += tmp1*cmplx.Conj(y[jy]) + tmp2*cmplx.Conj(x[jx]) |
| jx += incX |
| jy += incY |
| } |
| aii := real(ap[kk+i]) + real(tmp1*cmplx.Conj(y[iy])) + real(tmp2*cmplx.Conj(x[ix])) |
| ap[kk+i] = complex(aii, 0) |
| } else { |
| ap[kk+i] = complex(real(ap[kk+i]), 0) |
| } |
| ix += incX |
| iy += incY |
| kk += i + 1 |
| } |
| } |
| } |
| |
| // Ctbmv performs one of the matrix-vector operations |
| // x = A * x if trans = blas.NoTrans |
| // x = Aᵀ * x if trans = blas.Trans |
| // x = Aᴴ * x if trans = blas.ConjTrans |
| // where x is an n element vector and A is an n×n triangular band matrix, with |
| // (k+1) diagonals. |
| // |
| // Complex64 implementations are autogenerated and not directly tested. |
| func (Implementation) Ctbmv(uplo blas.Uplo, trans blas.Transpose, diag blas.Diag, n, k int, a []complex64, lda int, x []complex64, incX int) { |
| switch trans { |
| default: |
| panic(badTranspose) |
| case blas.NoTrans, blas.Trans, blas.ConjTrans: |
| } |
| switch uplo { |
| default: |
| panic(badUplo) |
| case blas.Upper, blas.Lower: |
| } |
| switch diag { |
| default: |
| panic(badDiag) |
| case blas.NonUnit, blas.Unit: |
| } |
| if n < 0 { |
| panic(nLT0) |
| } |
| if k < 0 { |
| panic(kLT0) |
| } |
| if lda < k+1 { |
| panic(badLdA) |
| } |
| if incX == 0 { |
| panic(zeroIncX) |
| } |
| |
| // Quick return if possible. |
| if n == 0 { |
| return |
| } |
| |
| // For zero matrix size the following slice length checks are trivially satisfied. |
| if len(a) < lda*(n-1)+k+1 { |
| panic(shortA) |
| } |
| if (incX > 0 && len(x) <= (n-1)*incX) || (incX < 0 && len(x) <= (1-n)*incX) { |
| panic(shortX) |
| } |
| |
| // Set up start index in X. |
| var kx int |
| if incX < 0 { |
| kx = (1 - n) * incX |
| } |
| |
| switch trans { |
| case blas.NoTrans: |
| if uplo == blas.Upper { |
| if incX == 1 { |
| for i := 0; i < n; i++ { |
| xi := x[i] |
| if diag == blas.NonUnit { |
| xi *= a[i*lda] |
| } |
| kk := min(k, n-i-1) |
| for j, aij := range a[i*lda+1 : i*lda+kk+1] { |
| xi += x[i+j+1] * aij |
| } |
| x[i] = xi |
| } |
| } else { |
| ix := kx |
| for i := 0; i < n; i++ { |
| xi := x[ix] |
| if diag == blas.NonUnit { |
| xi *= a[i*lda] |
| } |
| kk := min(k, n-i-1) |
| jx := ix + incX |
| for _, aij := range a[i*lda+1 : i*lda+kk+1] { |
| xi += x[jx] * aij |
| jx += incX |
| } |
| x[ix] = xi |
| ix += incX |
| } |
| } |
| } else { |
| if incX == 1 { |
| for i := n - 1; i >= 0; i-- { |
| xi := x[i] |
| if diag == blas.NonUnit { |
| xi *= a[i*lda+k] |
| } |
| kk := min(k, i) |
| for j, aij := range a[i*lda+k-kk : i*lda+k] { |
| xi += x[i-kk+j] * aij |
| } |
| x[i] = xi |
| } |
| } else { |
| ix := kx + (n-1)*incX |
| for i := n - 1; i >= 0; i-- { |
| xi := x[ix] |
| if diag == blas.NonUnit { |
| xi *= a[i*lda+k] |
| } |
| kk := min(k, i) |
| jx := ix - kk*incX |
| for _, aij := range a[i*lda+k-kk : i*lda+k] { |
| xi += x[jx] * aij |
| jx += incX |
| } |
| x[ix] = xi |
| ix -= incX |
| } |
| } |
| } |
| case blas.Trans: |
| if uplo == blas.Upper { |
| if incX == 1 { |
| for i := n - 1; i >= 0; i-- { |
| kk := min(k, n-i-1) |
| xi := x[i] |
| for j, aij := range a[i*lda+1 : i*lda+kk+1] { |
| x[i+j+1] += xi * aij |
| } |
| if diag == blas.NonUnit { |
| x[i] *= a[i*lda] |
| } |
| } |
| } else { |
| ix := kx + (n-1)*incX |
| for i := n - 1; i >= 0; i-- { |
| kk := min(k, n-i-1) |
| jx := ix + incX |
| xi := x[ix] |
| for _, aij := range a[i*lda+1 : i*lda+kk+1] { |
| x[jx] += xi * aij |
| jx += incX |
| } |
| if diag == blas.NonUnit { |
| x[ix] *= a[i*lda] |
| } |
| ix -= incX |
| } |
| } |
| } else { |
| if incX == 1 { |
| for i := 0; i < n; i++ { |
| kk := min(k, i) |
| xi := x[i] |
| for j, aij := range a[i*lda+k-kk : i*lda+k] { |
| x[i-kk+j] += xi * aij |
| } |
| if diag == blas.NonUnit { |
| x[i] *= a[i*lda+k] |
| } |
| } |
| } else { |
| ix := kx |
| for i := 0; i < n; i++ { |
| kk := min(k, i) |
| jx := ix - kk*incX |
| xi := x[ix] |
| for _, aij := range a[i*lda+k-kk : i*lda+k] { |
| x[jx] += xi * aij |
| jx += incX |
| } |
| if diag == blas.NonUnit { |
| x[ix] *= a[i*lda+k] |
| } |
| ix += incX |
| } |
| } |
| } |
| case blas.ConjTrans: |
| if uplo == blas.Upper { |
| if incX == 1 { |
| for i := n - 1; i >= 0; i-- { |
| kk := min(k, n-i-1) |
| xi := x[i] |
| for j, aij := range a[i*lda+1 : i*lda+kk+1] { |
| x[i+j+1] += xi * cmplx.Conj(aij) |
| } |
| if diag == blas.NonUnit { |
| x[i] *= cmplx.Conj(a[i*lda]) |
| } |
| } |
| } else { |
| ix := kx + (n-1)*incX |
| for i := n - 1; i >= 0; i-- { |
| kk := min(k, n-i-1) |
| jx := ix + incX |
| xi := x[ix] |
| for _, aij := range a[i*lda+1 : i*lda+kk+1] { |
| x[jx] += xi * cmplx.Conj(aij) |
| jx += incX |
| } |
| if diag == blas.NonUnit { |
| x[ix] *= cmplx.Conj(a[i*lda]) |
| } |
| ix -= incX |
| } |
| } |
| } else { |
| if incX == 1 { |
| for i := 0; i < n; i++ { |
| kk := min(k, i) |
| xi := x[i] |
| for j, aij := range a[i*lda+k-kk : i*lda+k] { |
| x[i-kk+j] += xi * cmplx.Conj(aij) |
| } |
| if diag == blas.NonUnit { |
| x[i] *= cmplx.Conj(a[i*lda+k]) |
| } |
| } |
| } else { |
| ix := kx |
| for i := 0; i < n; i++ { |
| kk := min(k, i) |
| jx := ix - kk*incX |
| xi := x[ix] |
| for _, aij := range a[i*lda+k-kk : i*lda+k] { |
| x[jx] += xi * cmplx.Conj(aij) |
| jx += incX |
| } |
| if diag == blas.NonUnit { |
| x[ix] *= cmplx.Conj(a[i*lda+k]) |
| } |
| ix += incX |
| } |
| } |
| } |
| } |
| } |
| |
| // Ctbsv solves one of the systems of equations |
| // A * x = b if trans == blas.NoTrans |
| // Aᵀ * x = b if trans == blas.Trans |
| // Aᴴ * x = b if trans == blas.ConjTrans |
| // where b and x are n element vectors and A is an n×n triangular band matrix |
| // with (k+1) diagonals. |
| // |
| // On entry, x contains the values of b, and the solution is |
| // stored in-place into x. |
| // |
| // No test for singularity or near-singularity is included in this |
| // routine. Such tests must be performed before calling this routine. |
| // |
| // Complex64 implementations are autogenerated and not directly tested. |
| func (Implementation) Ctbsv(uplo blas.Uplo, trans blas.Transpose, diag blas.Diag, n, k int, a []complex64, lda int, x []complex64, incX int) { |
| switch trans { |
| default: |
| panic(badTranspose) |
| case blas.NoTrans, blas.Trans, blas.ConjTrans: |
| } |
| switch uplo { |
| default: |
| panic(badUplo) |
| case blas.Upper, blas.Lower: |
| } |
| switch diag { |
| default: |
| panic(badDiag) |
| case blas.NonUnit, blas.Unit: |
| } |
| if n < 0 { |
| panic(nLT0) |
| } |
| if k < 0 { |
| panic(kLT0) |
| } |
| if lda < k+1 { |
| panic(badLdA) |
| } |
| if incX == 0 { |
| panic(zeroIncX) |
| } |
| |
| // Quick return if possible. |
| if n == 0 { |
| return |
| } |
| |
| // For zero matrix size the following slice length checks are trivially satisfied. |
| if len(a) < lda*(n-1)+k+1 { |
| panic(shortA) |
| } |
| if (incX > 0 && len(x) <= (n-1)*incX) || (incX < 0 && len(x) <= (1-n)*incX) { |
| panic(shortX) |
| } |
| |
| // Set up start index in X. |
| var kx int |
| if incX < 0 { |
| kx = (1 - n) * incX |
| } |
| |
| switch trans { |
| case blas.NoTrans: |
| if uplo == blas.Upper { |
| if incX == 1 { |
| for i := n - 1; i >= 0; i-- { |
| kk := min(k, n-i-1) |
| var sum complex64 |
| for j, aij := range a[i*lda+1 : i*lda+kk+1] { |
| sum += x[i+1+j] * aij |
| } |
| x[i] -= sum |
| if diag == blas.NonUnit { |
| x[i] /= a[i*lda] |
| } |
| } |
| } else { |
| ix := kx + (n-1)*incX |
| for i := n - 1; i >= 0; i-- { |
| kk := min(k, n-i-1) |
| var sum complex64 |
| jx := ix + incX |
| for _, aij := range a[i*lda+1 : i*lda+kk+1] { |
| sum += x[jx] * aij |
| jx += incX |
| } |
| x[ix] -= sum |
| if diag == blas.NonUnit { |
| x[ix] /= a[i*lda] |
| } |
| ix -= incX |
| } |
| } |
| } else { |
| if incX == 1 { |
| for i := 0; i < n; i++ { |
| kk := min(k, i) |
| var sum complex64 |
| for j, aij := range a[i*lda+k-kk : i*lda+k] { |
| sum += x[i-kk+j] * aij |
| } |
| x[i] -= sum |
| if diag == blas.NonUnit { |
| x[i] /= a[i*lda+k] |
| } |
| } |
| } else { |
| ix := kx |
| for i := 0; i < n; i++ { |
| kk := min(k, i) |
| var sum complex64 |
| jx := ix - kk*incX |
| for _, aij := range a[i*lda+k-kk : i*lda+k] { |
| sum += x[jx] * aij |
| jx += incX |
| } |
| x[ix] -= sum |
| if diag == blas.NonUnit { |
| x[ix] /= a[i*lda+k] |
| } |
| ix += incX |
| } |
| } |
| } |
| case blas.Trans: |
| if uplo == blas.Upper { |
| if incX == 1 { |
| for i := 0; i < n; i++ { |
| if diag == blas.NonUnit { |
| x[i] /= a[i*lda] |
| } |
| kk := min(k, n-i-1) |
| xi := x[i] |
| for j, aij := range a[i*lda+1 : i*lda+kk+1] { |
| x[i+1+j] -= xi * aij |
| } |
| } |
| } else { |
| ix := kx |
| for i := 0; i < n; i++ { |
| if diag == blas.NonUnit { |
| x[ix] /= a[i*lda] |
| } |
| kk := min(k, n-i-1) |
| xi := x[ix] |
| jx := ix + incX |
| for _, aij := range a[i*lda+1 : i*lda+kk+1] { |
| x[jx] -= xi * aij |
| jx += incX |
| } |
| ix += incX |
| } |
| } |
| } else { |
| if incX == 1 { |
| for i := n - 1; i >= 0; i-- { |
| if diag == blas.NonUnit { |
| x[i] /= a[i*lda+k] |
| } |
| kk := min(k, i) |
| xi := x[i] |
| for j, aij := range a[i*lda+k-kk : i*lda+k] { |
| x[i-kk+j] -= xi * aij |
| } |
| } |
| } else { |
| ix := kx + (n-1)*incX |
| for i := n - 1; i >= 0; i-- { |
| if diag == blas.NonUnit { |
| x[ix] /= a[i*lda+k] |
| } |
| kk := min(k, i) |
| xi := x[ix] |
| jx := ix - kk*incX |
| for _, aij := range a[i*lda+k-kk : i*lda+k] { |
| x[jx] -= xi * aij |
| jx += incX |
| } |
| ix -= incX |
| } |
| } |
| } |
| case blas.ConjTrans: |
| if uplo == blas.Upper { |
| if incX == 1 { |
| for i := 0; i < n; i++ { |
| if diag == blas.NonUnit { |
| x[i] /= cmplx.Conj(a[i*lda]) |
| } |
| kk := min(k, n-i-1) |
| xi := x[i] |
| for j, aij := range a[i*lda+1 : i*lda+kk+1] { |
| x[i+1+j] -= xi * cmplx.Conj(aij) |
| } |
| } |
| } else { |
| ix := kx |
| for i := 0; i < n; i++ { |
| if diag == blas.NonUnit { |
| x[ix] /= cmplx.Conj(a[i*lda]) |
| } |
| kk := min(k, n-i-1) |
| xi := x[ix] |
| jx := ix + incX |
| for _, aij := range a[i*lda+1 : i*lda+kk+1] { |
| x[jx] -= xi * cmplx.Conj(aij) |
| jx += incX |
| } |
| ix += incX |
| } |
| } |
| } else { |
| if incX == 1 { |
| for i := n - 1; i >= 0; i-- { |
| if diag == blas.NonUnit { |
| x[i] /= cmplx.Conj(a[i*lda+k]) |
| } |
| kk := min(k, i) |
| xi := x[i] |
| for j, aij := range a[i*lda+k-kk : i*lda+k] { |
| x[i-kk+j] -= xi * cmplx.Conj(aij) |
| } |
| } |
| } else { |
| ix := kx + (n-1)*incX |
| for i := n - 1; i >= 0; i-- { |
| if diag == blas.NonUnit { |
| x[ix] /= cmplx.Conj(a[i*lda+k]) |
| } |
| kk := min(k, i) |
| xi := x[ix] |
| jx := ix - kk*incX |
| for _, aij := range a[i*lda+k-kk : i*lda+k] { |
| x[jx] -= xi * cmplx.Conj(aij) |
| jx += incX |
| } |
| ix -= incX |
| } |
| } |
| } |
| } |
| } |
| |
| // Ctpmv performs one of the matrix-vector operations |
| // x = A * x if trans = blas.NoTrans |
| // x = Aᵀ * x if trans = blas.Trans |
| // x = Aᴴ * x if trans = blas.ConjTrans |
| // where x is an n element vector and A is an n×n triangular matrix, supplied in |
| // packed form. |
| // |
| // Complex64 implementations are autogenerated and not directly tested. |
| func (Implementation) Ctpmv(uplo blas.Uplo, trans blas.Transpose, diag blas.Diag, n int, ap []complex64, x []complex64, incX int) { |
| switch uplo { |
| default: |
| panic(badUplo) |
| case blas.Upper, blas.Lower: |
| } |
| switch trans { |
| default: |
| panic(badTranspose) |
| case blas.NoTrans, blas.Trans, blas.ConjTrans: |
| } |
| switch diag { |
| default: |
| panic(badDiag) |
| case blas.NonUnit, blas.Unit: |
| } |
| if n < 0 { |
| panic(nLT0) |
| } |
| if incX == 0 { |
| panic(zeroIncX) |
| } |
| |
| // Quick return if possible. |
| if n == 0 { |
| return |
| } |
| |
| // For zero matrix size the following slice length checks are trivially satisfied. |
| if len(ap) < n*(n+1)/2 { |
| panic(shortAP) |
| } |
| if (incX > 0 && len(x) <= (n-1)*incX) || (incX < 0 && len(x) <= (1-n)*incX) { |
| panic(shortX) |
| } |
| |
| // Set up start index in X. |
| var kx int |
| if incX < 0 { |
| kx = (1 - n) * incX |
| } |
| |
| // The elements of A are accessed sequentially with one pass through A. |
| |
| if trans == blas.NoTrans { |
| // Form x = A*x. |
| if uplo == blas.Upper { |
| // kk points to the current diagonal element in ap. |
| kk := 0 |
| if incX == 1 { |
| x = x[:n] |
| for i := range x { |
| if diag == blas.NonUnit { |
| x[i] *= ap[kk] |
| } |
| if n-i-1 > 0 { |
| x[i] += c64.DotuUnitary(ap[kk+1:kk+n-i], x[i+1:]) |
| } |
| kk += n - i |
| } |
| } else { |
| ix := kx |
| for i := 0; i < n; i++ { |
| if diag == blas.NonUnit { |
| x[ix] *= ap[kk] |
| } |
| if n-i-1 > 0 { |
| x[ix] += c64.DotuInc(ap[kk+1:kk+n-i], x, uintptr(n-i-1), 1, uintptr(incX), 0, uintptr(ix+incX)) |
| } |
| ix += incX |
| kk += n - i |
| } |
| } |
| } else { |
| // kk points to the beginning of current row in ap. |
| kk := n*(n+1)/2 - n |
| if incX == 1 { |
| for i := n - 1; i >= 0; i-- { |
| if diag == blas.NonUnit { |
| x[i] *= ap[kk+i] |
| } |
| if i > 0 { |
| x[i] += c64.DotuUnitary(ap[kk:kk+i], x[:i]) |
| } |
| kk -= i |
| } |
| } else { |
| ix := kx + (n-1)*incX |
| for i := n - 1; i >= 0; i-- { |
| if diag == blas.NonUnit { |
| x[ix] *= ap[kk+i] |
| } |
| if i > 0 { |
| x[ix] += c64.DotuInc(ap[kk:kk+i], x, uintptr(i), 1, uintptr(incX), 0, uintptr(kx)) |
| } |
| ix -= incX |
| kk -= i |
| } |
| } |
| } |
| return |
| } |
| |
| if trans == blas.Trans { |
| // Form x = Aᵀ*x. |
| if uplo == blas.Upper { |
| // kk points to the current diagonal element in ap. |
| kk := n*(n+1)/2 - 1 |
| if incX == 1 { |
| for i := n - 1; i >= 0; i-- { |
| xi := x[i] |
| if diag == blas.NonUnit { |
| x[i] *= ap[kk] |
| } |
| if n-i-1 > 0 { |
| c64.AxpyUnitary(xi, ap[kk+1:kk+n-i], x[i+1:n]) |
| } |
| kk -= n - i + 1 |
| } |
| } else { |
| ix := kx + (n-1)*incX |
| for i := n - 1; i >= 0; i-- { |
| xi := x[ix] |
| if diag == blas.NonUnit { |
| x[ix] *= ap[kk] |
| } |
| if n-i-1 > 0 { |
| c64.AxpyInc(xi, ap[kk+1:kk+n-i], x, uintptr(n-i-1), 1, uintptr(incX), 0, uintptr(ix+incX)) |
| } |
| ix -= incX |
| kk -= n - i + 1 |
| } |
| } |
| } else { |
| // kk points to the beginning of current row in ap. |
| kk := 0 |
| if incX == 1 { |
| x = x[:n] |
| for i := range x { |
| if i > 0 { |
| c64.AxpyUnitary(x[i], ap[kk:kk+i], x[:i]) |
| } |
| if diag == blas.NonUnit { |
| x[i] *= ap[kk+i] |
| } |
| kk += i + 1 |
| } |
| } else { |
| ix := kx |
| for i := 0; i < n; i++ { |
| if i > 0 { |
| c64.AxpyInc(x[ix], ap[kk:kk+i], x, uintptr(i), 1, uintptr(incX), 0, uintptr(kx)) |
| } |
| if diag == blas.NonUnit { |
| x[ix] *= ap[kk+i] |
| } |
| ix += incX |
| kk += i + 1 |
| } |
| } |
| } |
| return |
| } |
| |
| // Form x = Aᴴ*x. |
| if uplo == blas.Upper { |
| // kk points to the current diagonal element in ap. |
| kk := n*(n+1)/2 - 1 |
| if incX == 1 { |
| for i := n - 1; i >= 0; i-- { |
| xi := x[i] |
| if diag == blas.NonUnit { |
| x[i] *= cmplx.Conj(ap[kk]) |
| } |
| k := kk + 1 |
| for j := i + 1; j < n; j++ { |
| x[j] += xi * cmplx.Conj(ap[k]) |
| k++ |
| } |
| kk -= n - i + 1 |
| } |
| } else { |
| ix := kx + (n-1)*incX |
| for i := n - 1; i >= 0; i-- { |
| xi := x[ix] |
| if diag == blas.NonUnit { |
| x[ix] *= cmplx.Conj(ap[kk]) |
| } |
| jx := ix + incX |
| k := kk + 1 |
| for j := i + 1; j < n; j++ { |
| x[jx] += xi * cmplx.Conj(ap[k]) |
| jx += incX |
| k++ |
| } |
| ix -= incX |
| kk -= n - i + 1 |
| } |
| } |
| } else { |
| // kk points to the beginning of current row in ap. |
| kk := 0 |
| if incX == 1 { |
| x = x[:n] |
| for i, xi := range x { |
| for j := 0; j < i; j++ { |
| x[j] += xi * cmplx.Conj(ap[kk+j]) |
| } |
| if diag == blas.NonUnit { |
| x[i] *= cmplx.Conj(ap[kk+i]) |
| } |
| kk += i + 1 |
| } |
| } else { |
| ix := kx |
| for i := 0; i < n; i++ { |
| xi := x[ix] |
| jx := kx |
| for j := 0; j < i; j++ { |
| x[jx] += xi * cmplx.Conj(ap[kk+j]) |
| jx += incX |
| } |
| if diag == blas.NonUnit { |
| x[ix] *= cmplx.Conj(ap[kk+i]) |
| } |
| ix += incX |
| kk += i + 1 |
| } |
| } |
| } |
| } |
| |
| // Ctpsv solves one of the systems of equations |
| // A * x = b if trans == blas.NoTrans |
| // Aᵀ * x = b if trans == blas.Trans |
| // Aᴴ * x = b if trans == blas.ConjTrans |
| // where b and x are n element vectors and A is an n×n triangular matrix in |
| // packed form. |
| // |
| // On entry, x contains the values of b, and the solution is |
| // stored in-place into x. |
| // |
| // No test for singularity or near-singularity is included in this |
| // routine. Such tests must be performed before calling this routine. |
| // |
| // Complex64 implementations are autogenerated and not directly tested. |
| func (Implementation) Ctpsv(uplo blas.Uplo, trans blas.Transpose, diag blas.Diag, n int, ap []complex64, x []complex64, incX int) { |
| switch uplo { |
| default: |
| panic(badUplo) |
| case blas.Upper, blas.Lower: |
| } |
| switch trans { |
| default: |
| panic(badTranspose) |
| case blas.NoTrans, blas.Trans, blas.ConjTrans: |
| } |
| switch diag { |
| default: |
| panic(badDiag) |
| case blas.NonUnit, blas.Unit: |
| } |
| if n < 0 { |
| panic(nLT0) |
| } |
| if incX == 0 { |
| panic(zeroIncX) |
| } |
| |
| // Quick return if possible. |
| if n == 0 { |
| return |
| } |
| |
| // For zero matrix size the following slice length checks are trivially satisfied. |
| if len(ap) < n*(n+1)/2 { |
| panic(shortAP) |
| } |
| if (incX > 0 && len(x) <= (n-1)*incX) || (incX < 0 && len(x) <= (1-n)*incX) { |
| panic(shortX) |
| } |
| |
| // Set up start index in X. |
| var kx int |
| if incX < 0 { |
| kx = (1 - n) * incX |
| } |
| |
| // The elements of A are accessed sequentially with one pass through ap. |
| |
| if trans == blas.NoTrans { |
| // Form x = inv(A)*x. |
| if uplo == blas.Upper { |
| kk := n*(n+1)/2 - 1 |
| if incX == 1 { |
| for i := n - 1; i >= 0; i-- { |
| aii := ap[kk] |
| if n-i-1 > 0 { |
| x[i] -= c64.DotuUnitary(x[i+1:n], ap[kk+1:kk+n-i]) |
| } |
| if diag == blas.NonUnit { |
| x[i] /= aii |
| } |
| kk -= n - i + 1 |
| } |
| } else { |
| ix := kx + (n-1)*incX |
| for i := n - 1; i >= 0; i-- { |
| aii := ap[kk] |
| if n-i-1 > 0 { |
| x[ix] -= c64.DotuInc(x, ap[kk+1:kk+n-i], uintptr(n-i-1), uintptr(incX), 1, uintptr(ix+incX), 0) |
| } |
| if diag == blas.NonUnit { |
| x[ix] /= aii |
| } |
| ix -= incX |
| kk -= n - i + 1 |
| } |
| } |
| } else { |
| kk := 0 |
| if incX == 1 { |
| for i := 0; i < n; i++ { |
| if i > 0 { |
| x[i] -= c64.DotuUnitary(x[:i], ap[kk:kk+i]) |
| } |
| if diag == blas.NonUnit { |
| x[i] /= ap[kk+i] |
| } |
| kk += i + 1 |
| } |
| } else { |
| ix := kx |
| for i := 0; i < n; i++ { |
| if i > 0 { |
| x[ix] -= c64.DotuInc(x, ap[kk:kk+i], uintptr(i), uintptr(incX), 1, uintptr(kx), 0) |
| } |
| if diag == blas.NonUnit { |
| x[ix] /= ap[kk+i] |
| } |
| ix += incX |
| kk += i + 1 |
| } |
| } |
| } |
| return |
| } |
| |
| if trans == blas.Trans { |
| // Form x = inv(Aᵀ)*x. |
| if uplo == blas.Upper { |
| kk := 0 |
| if incX == 1 { |
| for j := 0; j < n; j++ { |
| if diag == blas.NonUnit { |
| x[j] /= ap[kk] |
| } |
| if n-j-1 > 0 { |
| c64.AxpyUnitary(-x[j], ap[kk+1:kk+n-j], x[j+1:n]) |
| } |
| kk += n - j |
| } |
| } else { |
| jx := kx |
| for j := 0; j < n; j++ { |
| if diag == blas.NonUnit { |
| x[jx] /= ap[kk] |
| } |
| if n-j-1 > 0 { |
| c64.AxpyInc(-x[jx], ap[kk+1:kk+n-j], x, uintptr(n-j-1), 1, uintptr(incX), 0, uintptr(jx+incX)) |
| } |
| jx += incX |
| kk += n - j |
| } |
| } |
| } else { |
| kk := n*(n+1)/2 - n |
| if incX == 1 { |
| for j := n - 1; j >= 0; j-- { |
| if diag == blas.NonUnit { |
| x[j] /= ap[kk+j] |
| } |
| if j > 0 { |
| c64.AxpyUnitary(-x[j], ap[kk:kk+j], x[:j]) |
| } |
| kk -= j |
| } |
| } else { |
| jx := kx + (n-1)*incX |
| for j := n - 1; j >= 0; j-- { |
| if diag == blas.NonUnit { |
| x[jx] /= ap[kk+j] |
| } |
| if j > 0 { |
| c64.AxpyInc(-x[jx], ap[kk:kk+j], x, uintptr(j), 1, uintptr(incX), 0, uintptr(kx)) |
| } |
| jx -= incX |
| kk -= j |
| } |
| } |
| } |
| return |
| } |
| |
| // Form x = inv(Aᴴ)*x. |
| if uplo == blas.Upper { |
| kk := 0 |
| if incX == 1 { |
| for j := 0; j < n; j++ { |
| if diag == blas.NonUnit { |
| x[j] /= cmplx.Conj(ap[kk]) |
| } |
| xj := x[j] |
| k := kk + 1 |
| for i := j + 1; i < n; i++ { |
| x[i] -= xj * cmplx.Conj(ap[k]) |
| k++ |
| } |
| kk += n - j |
| } |
| } else { |
| jx := kx |
| for j := 0; j < n; j++ { |
| if diag == blas.NonUnit { |
| x[jx] /= cmplx.Conj(ap[kk]) |
| } |
| xj := x[jx] |
| ix := jx + incX |
| k := kk + 1 |
| for i := j + 1; i < n; i++ { |
| x[ix] -= xj * cmplx.Conj(ap[k]) |
| ix += incX |
| k++ |
| } |
| jx += incX |
| kk += n - j |
| } |
| } |
| } else { |
| kk := n*(n+1)/2 - n |
| if incX == 1 { |
| for j := n - 1; j >= 0; j-- { |
| if diag == blas.NonUnit { |
| x[j] /= cmplx.Conj(ap[kk+j]) |
| } |
| xj := x[j] |
| for i := 0; i < j; i++ { |
| x[i] -= xj * cmplx.Conj(ap[kk+i]) |
| } |
| kk -= j |
| } |
| } else { |
| jx := kx + (n-1)*incX |
| for j := n - 1; j >= 0; j-- { |
| if diag == blas.NonUnit { |
| x[jx] /= cmplx.Conj(ap[kk+j]) |
| } |
| xj := x[jx] |
| ix := kx |
| for i := 0; i < j; i++ { |
| x[ix] -= xj * cmplx.Conj(ap[kk+i]) |
| ix += incX |
| } |
| jx -= incX |
| kk -= j |
| } |
| } |
| } |
| } |
| |
| // Ctrmv performs one of the matrix-vector operations |
| // x = A * x if trans = blas.NoTrans |
| // x = Aᵀ * x if trans = blas.Trans |
| // x = Aᴴ * x if trans = blas.ConjTrans |
| // where x is a vector, and A is an n×n triangular matrix. |
| // |
| // Complex64 implementations are autogenerated and not directly tested. |
| func (Implementation) Ctrmv(uplo blas.Uplo, trans blas.Transpose, diag blas.Diag, n int, a []complex64, lda int, x []complex64, incX int) { |
| switch trans { |
| default: |
| panic(badTranspose) |
| case blas.NoTrans, blas.Trans, blas.ConjTrans: |
| } |
| switch uplo { |
| default: |
| panic(badUplo) |
| case blas.Upper, blas.Lower: |
| } |
| switch diag { |
| default: |
| panic(badDiag) |
| case blas.NonUnit, blas.Unit: |
| } |
| if n < 0 { |
| panic(nLT0) |
| } |
| if lda < max(1, n) { |
| panic(badLdA) |
| } |
| if incX == 0 { |
| panic(zeroIncX) |
| } |
| |
| // Quick return if possible. |
| if n == 0 { |
| return |
| } |
| |
| // For zero matrix size the following slice length checks are trivially satisfied. |
| if len(a) < lda*(n-1)+n { |
| panic(shortA) |
| } |
| if (incX > 0 && len(x) <= (n-1)*incX) || (incX < 0 && len(x) <= (1-n)*incX) { |
| panic(shortX) |
| } |
| |
| // Set up start index in X. |
| var kx int |
| if incX < 0 { |
| kx = (1 - n) * incX |
| } |
| |
| // The elements of A are accessed sequentially with one pass through A. |
| |
| if trans == blas.NoTrans { |
| // Form x = A*x. |
| if uplo == blas.Upper { |
| if incX == 1 { |
| for i := 0; i < n; i++ { |
| if diag == blas.NonUnit { |
| x[i] *= a[i*lda+i] |
| } |
| if n-i-1 > 0 { |
| x[i] += c64.DotuUnitary(a[i*lda+i+1:i*lda+n], x[i+1:n]) |
| } |
| } |
| } else { |
| ix := kx |
| for i := 0; i < n; i++ { |
| if diag == blas.NonUnit { |
| x[ix] *= a[i*lda+i] |
| } |
| if n-i-1 > 0 { |
| x[ix] += c64.DotuInc(a[i*lda+i+1:i*lda+n], x, uintptr(n-i-1), 1, uintptr(incX), 0, uintptr(ix+incX)) |
| } |
| ix += incX |
| } |
| } |
| } else { |
| if incX == 1 { |
| for i := n - 1; i >= 0; i-- { |
| if diag == blas.NonUnit { |
| x[i] *= a[i*lda+i] |
| } |
| if i > 0 { |
| x[i] += c64.DotuUnitary(a[i*lda:i*lda+i], x[:i]) |
| } |
| } |
| } else { |
| ix := kx + (n-1)*incX |
| for i := n - 1; i >= 0; i-- { |
| if diag == blas.NonUnit { |
| x[ix] *= a[i*lda+i] |
| } |
| if i > 0 { |
| x[ix] += c64.DotuInc(a[i*lda:i*lda+i], x, uintptr(i), 1, uintptr(incX), 0, uintptr(kx)) |
| } |
| ix -= incX |
| } |
| } |
| } |
| return |
| } |
| |
| if trans == blas.Trans { |
| // Form x = Aᵀ*x. |
| if uplo == blas.Upper { |
| if incX == 1 { |
| for i := n - 1; i >= 0; i-- { |
| xi := x[i] |
| if diag == blas.NonUnit { |
| x[i] *= a[i*lda+i] |
| } |
| if n-i-1 > 0 { |
| c64.AxpyUnitary(xi, a[i*lda+i+1:i*lda+n], x[i+1:n]) |
| } |
| } |
| } else { |
| ix := kx + (n-1)*incX |
| for i := n - 1; i >= 0; i-- { |
| xi := x[ix] |
| if diag == blas.NonUnit { |
| x[ix] *= a[i*lda+i] |
| } |
| if n-i-1 > 0 { |
| c64.AxpyInc(xi, a[i*lda+i+1:i*lda+n], x, uintptr(n-i-1), 1, uintptr(incX), 0, uintptr(ix+incX)) |
| } |
| ix -= incX |
| } |
| } |
| } else { |
| if incX == 1 { |
| for i := 0; i < n; i++ { |
| if i > 0 { |
| c64.AxpyUnitary(x[i], a[i*lda:i*lda+i], x[:i]) |
| } |
| if diag == blas.NonUnit { |
| x[i] *= a[i*lda+i] |
| } |
| } |
| } else { |
| ix := kx |
| for i := 0; i < n; i++ { |
| if i > 0 { |
| c64.AxpyInc(x[ix], a[i*lda:i*lda+i], x, uintptr(i), 1, uintptr(incX), 0, uintptr(kx)) |
| } |
| if diag == blas.NonUnit { |
| x[ix] *= a[i*lda+i] |
| } |
| ix += incX |
| } |
| } |
| } |
| return |
| } |
| |
| // Form x = Aᴴ*x. |
| if uplo == blas.Upper { |
| if incX == 1 { |
| for i := n - 1; i >= 0; i-- { |
| xi := x[i] |
| if diag == blas.NonUnit { |
| x[i] *= cmplx.Conj(a[i*lda+i]) |
| } |
| for j := i + 1; j < n; j++ { |
| x[j] += xi * cmplx.Conj(a[i*lda+j]) |
| } |
| } |
| } else { |
| ix := kx + (n-1)*incX |
| for i := n - 1; i >= 0; i-- { |
| xi := x[ix] |
| if diag == blas.NonUnit { |
| x[ix] *= cmplx.Conj(a[i*lda+i]) |
| } |
| jx := ix + incX |
| for j := i + 1; j < n; j++ { |
| x[jx] += xi * cmplx.Conj(a[i*lda+j]) |
| jx += incX |
| } |
| ix -= incX |
| } |
| } |
| } else { |
| if incX == 1 { |
| for i := 0; i < n; i++ { |
| for j := 0; j < i; j++ { |
| x[j] += x[i] * cmplx.Conj(a[i*lda+j]) |
| } |
| if diag == blas.NonUnit { |
| x[i] *= cmplx.Conj(a[i*lda+i]) |
| } |
| } |
| } else { |
| ix := kx |
| for i := 0; i < n; i++ { |
| jx := kx |
| for j := 0; j < i; j++ { |
| x[jx] += x[ix] * cmplx.Conj(a[i*lda+j]) |
| jx += incX |
| } |
| if diag == blas.NonUnit { |
| x[ix] *= cmplx.Conj(a[i*lda+i]) |
| } |
| ix += incX |
| } |
| } |
| } |
| } |
| |
| // Ctrsv solves one of the systems of equations |
| // A * x = b if trans == blas.NoTrans |
| // Aᵀ * x = b if trans == blas.Trans |
| // Aᴴ * x = b if trans == blas.ConjTrans |
| // where b and x are n element vectors and A is an n×n triangular matrix. |
| // |
| // On entry, x contains the values of b, and the solution is |
| // stored in-place into x. |
| // |
| // No test for singularity or near-singularity is included in this |
| // routine. Such tests must be performed before calling this routine. |
| // |
| // Complex64 implementations are autogenerated and not directly tested. |
| func (Implementation) Ctrsv(uplo blas.Uplo, trans blas.Transpose, diag blas.Diag, n int, a []complex64, lda int, x []complex64, incX int) { |
| switch trans { |
| default: |
| panic(badTranspose) |
| case blas.NoTrans, blas.Trans, blas.ConjTrans: |
| } |
| switch uplo { |
| default: |
| panic(badUplo) |
| case blas.Upper, blas.Lower: |
| } |
| switch diag { |
| default: |
| panic(badDiag) |
| case blas.NonUnit, blas.Unit: |
| } |
| if n < 0 { |
| panic(nLT0) |
| } |
| if lda < max(1, n) { |
| panic(badLdA) |
| } |
| if incX == 0 { |
| panic(zeroIncX) |
| } |
| |
| // Quick return if possible. |
| if n == 0 { |
| return |
| } |
| |
| // For zero matrix size the following slice length checks are trivially satisfied. |
| if len(a) < lda*(n-1)+n { |
| panic(shortA) |
| } |
| if (incX > 0 && len(x) <= (n-1)*incX) || (incX < 0 && len(x) <= (1-n)*incX) { |
| panic(shortX) |
| } |
| |
| // Set up start index in X. |
| var kx int |
| if incX < 0 { |
| kx = (1 - n) * incX |
| } |
| |
| // The elements of A are accessed sequentially with one pass through A. |
| |
| if trans == blas.NoTrans { |
| // Form x = inv(A)*x. |
| if uplo == blas.Upper { |
| if incX == 1 { |
| for i := n - 1; i >= 0; i-- { |
| aii := a[i*lda+i] |
| if n-i-1 > 0 { |
| x[i] -= c64.DotuUnitary(x[i+1:n], a[i*lda+i+1:i*lda+n]) |
| } |
| if diag == blas.NonUnit { |
| x[i] /= aii |
| } |
| } |
| } else { |
| ix := kx + (n-1)*incX |
| for i := n - 1; i >= 0; i-- { |
| aii := a[i*lda+i] |
| if n-i-1 > 0 { |
| x[ix] -= c64.DotuInc(x, a[i*lda+i+1:i*lda+n], uintptr(n-i-1), uintptr(incX), 1, uintptr(ix+incX), 0) |
| } |
| if diag == blas.NonUnit { |
| x[ix] /= aii |
| } |
| ix -= incX |
| } |
| } |
| } else { |
| if incX == 1 { |
| for i := 0; i < n; i++ { |
| if i > 0 { |
| x[i] -= c64.DotuUnitary(x[:i], a[i*lda:i*lda+i]) |
| } |
| if diag == blas.NonUnit { |
| x[i] /= a[i*lda+i] |
| } |
| } |
| } else { |
| ix := kx |
| for i := 0; i < n; i++ { |
| if i > 0 { |
| x[ix] -= c64.DotuInc(x, a[i*lda:i*lda+i], uintptr(i), uintptr(incX), 1, uintptr(kx), 0) |
| } |
| if diag == blas.NonUnit { |
| x[ix] /= a[i*lda+i] |
| } |
| ix += incX |
| } |
| } |
| } |
| return |
| } |
| |
| if trans == blas.Trans { |
| // Form x = inv(Aᵀ)*x. |
| if uplo == blas.Upper { |
| if incX == 1 { |
| for j := 0; j < n; j++ { |
| if diag == blas.NonUnit { |
| x[j] /= a[j*lda+j] |
| } |
| if n-j-1 > 0 { |
| c64.AxpyUnitary(-x[j], a[j*lda+j+1:j*lda+n], x[j+1:n]) |
| } |
| } |
| } else { |
| jx := kx |
| for j := 0; j < n; j++ { |
| if diag == blas.NonUnit { |
| x[jx] /= a[j*lda+j] |
| } |
| if n-j-1 > 0 { |
| c64.AxpyInc(-x[jx], a[j*lda+j+1:j*lda+n], x, uintptr(n-j-1), 1, uintptr(incX), 0, uintptr(jx+incX)) |
| } |
| jx += incX |
| } |
| } |
| } else { |
| if incX == 1 { |
| for j := n - 1; j >= 0; j-- { |
| if diag == blas.NonUnit { |
| x[j] /= a[j*lda+j] |
| } |
| xj := x[j] |
| if j > 0 { |
| c64.AxpyUnitary(-xj, a[j*lda:j*lda+j], x[:j]) |
| } |
| } |
| } else { |
| jx := kx + (n-1)*incX |
| for j := n - 1; j >= 0; j-- { |
| if diag == blas.NonUnit { |
| x[jx] /= a[j*lda+j] |
| } |
| if j > 0 { |
| c64.AxpyInc(-x[jx], a[j*lda:j*lda+j], x, uintptr(j), 1, uintptr(incX), 0, uintptr(kx)) |
| } |
| jx -= incX |
| } |
| } |
| } |
| return |
| } |
| |
| // Form x = inv(Aᴴ)*x. |
| if uplo == blas.Upper { |
| if incX == 1 { |
| for j := 0; j < n; j++ { |
| if diag == blas.NonUnit { |
| x[j] /= cmplx.Conj(a[j*lda+j]) |
| } |
| xj := x[j] |
| for i := j + 1; i < n; i++ { |
| x[i] -= xj * cmplx.Conj(a[j*lda+i]) |
| } |
| } |
| } else { |
| jx := kx |
| for j := 0; j < n; j++ { |
| if diag == blas.NonUnit { |
| x[jx] /= cmplx.Conj(a[j*lda+j]) |
| } |
| xj := x[jx] |
| ix := jx + incX |
| for i := j + 1; i < n; i++ { |
| x[ix] -= xj * cmplx.Conj(a[j*lda+i]) |
| ix += incX |
| } |
| jx += incX |
| } |
| } |
| } else { |
| if incX == 1 { |
| for j := n - 1; j >= 0; j-- { |
| if diag == blas.NonUnit { |
| x[j] /= cmplx.Conj(a[j*lda+j]) |
| } |
| xj := x[j] |
| for i := 0; i < j; i++ { |
| x[i] -= xj * cmplx.Conj(a[j*lda+i]) |
| } |
| } |
| } else { |
| jx := kx + (n-1)*incX |
| for j := n - 1; j >= 0; j-- { |
| if diag == blas.NonUnit { |
| x[jx] /= cmplx.Conj(a[j*lda+j]) |
| } |
| xj := x[jx] |
| ix := kx |
| for i := 0; i < j; i++ { |
| x[ix] -= xj * cmplx.Conj(a[j*lda+i]) |
| ix += incX |
| } |
| jx -= incX |
| } |
| } |
| } |
| } |