2942 lines
62 KiB
Go
2942 lines
62 KiB
Go
// 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^T * x + beta * y if trans = blas.Trans
|
||
// y = alpha * A^H * 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^T * x + beta * y if trans = blas.Trans
|
||
// y = alpha * A^H * 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^T*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^H*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^H
|
||
// 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^T
|
||
// 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^H
|
||
// 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^H + conj(alpha) * y * x^H
|
||
// 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^H
|
||
// 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^H + conj(alpha) * y * x^H
|
||
// 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^T * x if trans = blas.Trans
|
||
// x = A^H * 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^T * x = b if trans == blas.Trans
|
||
// A^H * 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^T * x if trans = blas.Trans
|
||
// x = A^H * 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^T*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^H*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^T * x = b if trans == blas.Trans
|
||
// A^H * 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^T)*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^H)*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^T * x if trans = blas.Trans
|
||
// x = A^H * 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^T*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^H*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^T * x = b if trans == blas.Trans
|
||
// A^H * 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^T)*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^H)*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
|
||
}
|
||
}
|
||
}
|
||
}
|