// Copyright ©2014 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 (
	"gonum.org/v1/gonum/blas"
	"gonum.org/v1/gonum/internal/asm/f64"
)

var _ blas.Float64Level2 = Implementation{}

// Dgemv computes
//  y = alpha * a * x + beta * y if tA = blas.NoTrans
//  y = alpha * A^T * x + beta * y if tA = blas.Trans or blas.ConjTrans
// where A is an m×n dense matrix, x and y are vectors, and alpha is a scalar.
func (Implementation) Dgemv(tA blas.Transpose, m, n int, alpha float64, a []float64, lda int, x []float64, incX int, beta float64, y []float64, incY int) {
	if tA != blas.NoTrans && tA != blas.Trans && tA != blas.ConjTrans {
		panic(badTranspose)
	}
	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)
	}
	// Set up indexes
	lenX := m
	lenY := n
	if tA == blas.NoTrans {
		lenX = n
		lenY = m
	}
	if (incX > 0 && (lenX-1)*incX >= len(x)) || (incX < 0 && (1-lenX)*incX >= len(x)) {
		panic(badX)
	}
	if (incY > 0 && (lenY-1)*incY >= len(y)) || (incY < 0 && (1-lenY)*incY >= len(y)) {
		panic(badY)
	}
	if lda*(m-1)+n > len(a) || lda < max(1, n) {
		panic(badLdA)
	}

	// Quick return if possible
	if m == 0 || n == 0 || (alpha == 0 && beta == 1) {
		return
	}

	var kx, ky int
	if incX > 0 {
		kx = 0
	} else {
		kx = -(lenX - 1) * incX
	}
	if incY > 0 {
		ky = 0
	} else {
		ky = -(lenY - 1) * incY
	}

	// First form y := beta * y
	if incY > 0 {
		Implementation{}.Dscal(lenY, beta, y, incY)
	} else {
		Implementation{}.Dscal(lenY, beta, y, -incY)
	}

	if alpha == 0 {
		return
	}

	// Form y := alpha * A * x + y
	if tA == blas.NoTrans {
		if incX == 1 && incY == 1 {
			for i := 0; i < m; i++ {
				y[i] += alpha * f64.DotUnitary(a[lda*i:lda*i+n], x)
			}
			return
		}
		iy := ky
		for i := 0; i < m; i++ {
			y[iy] += alpha * f64.DotInc(x, a[lda*i:lda*i+n], uintptr(n), uintptr(incX), 1, uintptr(kx), 0)
			iy += incY
		}
		return
	}
	// Cases where a is transposed.
	if incX == 1 && incY == 1 {
		for i := 0; i < m; i++ {
			tmp := alpha * x[i]
			if tmp != 0 {
				f64.AxpyUnitaryTo(y, tmp, a[lda*i:lda*i+n], y)
			}
		}
		return
	}
	ix := kx
	for i := 0; i < m; i++ {
		tmp := alpha * x[ix]
		if tmp != 0 {
			f64.AxpyInc(tmp, a[lda*i:lda*i+n], y, uintptr(n), 1, uintptr(incY), 0, uintptr(ky))
		}
		ix += incX
	}
}

// Dger performs the rank-one operation
//  A += alpha * x * y^T
// where A is an m×n dense matrix, x and y are vectors, and alpha is a scalar.
func (Implementation) Dger(m, n int, alpha float64, x []float64, incX int, y []float64, incY int, a []float64, lda int) {
	// Check inputs
	if m < 0 {
		panic("m < 0")
	}
	if n < 0 {
		panic(negativeN)
	}
	if incX == 0 {
		panic(zeroIncX)
	}
	if incY == 0 {
		panic(zeroIncY)
	}
	if (incX > 0 && (m-1)*incX >= len(x)) || (incX < 0 && (1-m)*incX >= len(x)) {
		panic(badX)
	}
	if (incY > 0 && (n-1)*incY >= len(y)) || (incY < 0 && (1-n)*incY >= len(y)) {
		panic(badY)
	}
	if lda*(m-1)+n > len(a) || lda < max(1, n) {
		panic(badLdA)
	}
	if lda < max(1, n) {
		panic(badLdA)
	}

	// Quick return if possible
	if m == 0 || n == 0 || alpha == 0 {
		return
	}

	var ky, kx int
	if incY > 0 {
		ky = 0
	} else {
		ky = -(n - 1) * incY
	}

	if incX > 0 {
		kx = 0
	} else {
		kx = -(m - 1) * incX
	}

	if incX == 1 && incY == 1 {
		x = x[:m]
		y = y[:n]
		for i, xv := range x {
			tmp := alpha * xv
			if tmp != 0 {
				atmp := a[i*lda : i*lda+n]
				f64.AxpyUnitaryTo(atmp, tmp, y, atmp)
			}
		}
		return
	}

	ix := kx
	for i := 0; i < m; i++ {
		tmp := alpha * x[ix]
		if tmp != 0 {
			f64.AxpyInc(tmp, y, a[i*lda:i*lda+n], uintptr(n), uintptr(incY), 1, uintptr(ky), 0)
		}
		ix += incX
	}
}

// Dgbmv computes
//  y = alpha * A * x + beta * y if tA == blas.NoTrans
//  y = alpha * A^T * x + beta * y if tA == blas.Trans or blas.ConjTrans
// where a is an m×n band matrix kL subdiagonals and kU super-diagonals, and
// m and n refer to the size of the full dense matrix it represents.
// x and y are vectors, and alpha and beta are scalars.
func (Implementation) Dgbmv(tA blas.Transpose, m, n, kL, kU int, alpha float64, a []float64, lda int, x []float64, incX int, beta float64, y []float64, incY int) {
	if tA != blas.NoTrans && tA != blas.Trans && tA != blas.ConjTrans {
		panic(badTranspose)
	}
	if m < 0 {
		panic(mLT0)
	}
	if n < 0 {
		panic(nLT0)
	}
	if kL < 0 {
		panic(kLLT0)
	}
	if kL < 0 {
		panic(kULT0)
	}
	if lda < kL+kU+1 {
		panic(badLdA)
	}
	if incX == 0 {
		panic(zeroIncX)
	}
	if incY == 0 {
		panic(zeroIncY)
	}
	// Set up indexes
	lenX := m
	lenY := n
	if tA == blas.NoTrans {
		lenX = n
		lenY = m
	}
	if (incX > 0 && (lenX-1)*incX >= len(x)) || (incX < 0 && (1-lenX)*incX >= len(x)) {
		panic(badX)
	}
	if (incY > 0 && (lenY-1)*incY >= len(y)) || (incY < 0 && (1-lenY)*incY >= len(y)) {
		panic(badY)
	}
	if lda*(min(m, n+kL)-1)+kL+kU+1 > len(a) || lda < kL+kU+1 {
		panic(badLdA)
	}

	// Quick return if possible
	if m == 0 || n == 0 || (alpha == 0 && beta == 1) {
		return
	}

	var kx, ky int
	if incX > 0 {
		kx = 0
	} else {
		kx = -(lenX - 1) * incX
	}
	if incY > 0 {
		ky = 0
	} else {
		ky = -(lenY - 1) * incY
	}

	// First form y := beta * y
	if incY > 0 {
		Implementation{}.Dscal(lenY, beta, y, incY)
	} else {
		Implementation{}.Dscal(lenY, beta, y, -incY)
	}

	if alpha == 0 {
		return
	}

	// i and j are indices of the compacted banded matrix.
	// off is the offset into the dense matrix (off + j = densej)
	ld := min(m, n)
	nCol := kU + 1 + kL
	if tA == blas.NoTrans {
		iy := ky
		if incX == 1 {
			for i := 0; i < min(m, n+kL); i++ {
				l := max(0, kL-i)
				u := min(nCol, ld+kL-i)
				off := max(0, i-kL)
				atmp := a[i*lda+l : i*lda+u]
				xtmp := x[off : off+u-l]
				var sum float64
				for j, v := range atmp {
					sum += xtmp[j] * v
				}
				y[iy] += sum * alpha
				iy += incY
			}
			return
		}
		for i := 0; i < min(m, n+kL); i++ {
			l := max(0, kL-i)
			u := min(nCol, ld+kL-i)
			off := max(0, i-kL)
			atmp := a[i*lda+l : i*lda+u]
			jx := kx
			var sum float64
			for _, v := range atmp {
				sum += x[off*incX+jx] * v
				jx += incX
			}
			y[iy] += sum * alpha
			iy += incY
		}
		return
	}
	if incX == 1 {
		for i := 0; i < min(m, n+kL); i++ {
			l := max(0, kL-i)
			u := min(nCol, ld+kL-i)
			off := max(0, i-kL)
			atmp := a[i*lda+l : i*lda+u]
			tmp := alpha * x[i]
			jy := ky
			for _, v := range atmp {
				y[jy+off*incY] += tmp * v
				jy += incY
			}
		}
		return
	}
	ix := kx
	for i := 0; i < min(m, n+kL); i++ {
		l := max(0, kL-i)
		u := min(nCol, ld+kL-i)
		off := max(0, i-kL)
		atmp := a[i*lda+l : i*lda+u]
		tmp := alpha * x[ix]
		jy := ky
		for _, v := range atmp {
			y[jy+off*incY] += tmp * v
			jy += incY
		}
		ix += incX
	}
}

// Dtrmv computes
//  x = A * x if tA == blas.NoTrans
//  x = A^T * x if tA == blas.Trans or blas.ConjTrans
// A is an n×n Triangular matrix and x is a vector.
func (Implementation) Dtrmv(ul blas.Uplo, tA blas.Transpose, d blas.Diag, n int, a []float64, lda int, x []float64, incX int) {
	if ul != blas.Lower && ul != blas.Upper {
		panic(badUplo)
	}
	if tA != blas.NoTrans && tA != blas.Trans && tA != blas.ConjTrans {
		panic(badTranspose)
	}
	if d != blas.NonUnit && d != blas.Unit {
		panic(badDiag)
	}
	if n < 0 {
		panic(nLT0)
	}
	if lda < n {
		panic(badLdA)
	}
	if incX == 0 {
		panic(zeroIncX)
	}
	if (incX > 0 && (n-1)*incX >= len(x)) || (incX < 0 && (1-n)*incX >= len(x)) {
		panic(badX)
	}
	if lda*(n-1)+n > len(a) || lda < max(1, n) {
		panic(badLdA)
	}
	if n == 0 {
		return
	}
	nonUnit := d != blas.Unit
	if n == 1 {
		if nonUnit {
			x[0] *= a[0]
		}
		return
	}
	var kx int
	if incX <= 0 {
		kx = -(n - 1) * incX
	}
	if tA == blas.NoTrans {
		if ul == blas.Upper {
			if incX == 1 {
				for i := 0; i < n; i++ {
					ilda := i * lda
					var tmp float64
					if nonUnit {
						tmp = a[ilda+i] * x[i]
					} else {
						tmp = x[i]
					}
					xtmp := x[i+1:]
					x[i] = tmp + f64.DotUnitary(a[ilda+i+1:ilda+n], xtmp)
				}
				return
			}
			ix := kx
			for i := 0; i < n; i++ {
				ilda := i * lda
				var tmp float64
				if nonUnit {
					tmp = a[ilda+i] * x[ix]
				} else {
					tmp = x[ix]
				}
				x[ix] = tmp + f64.DotInc(x, a[ilda+i+1:ilda+n], uintptr(n-i-1), uintptr(incX), 1, uintptr(ix+incX), 0)
				ix += incX
			}
			return
		}
		if incX == 1 {
			for i := n - 1; i >= 0; i-- {
				ilda := i * lda
				var tmp float64
				if nonUnit {
					tmp += a[ilda+i] * x[i]
				} else {
					tmp = x[i]
				}
				x[i] = tmp + f64.DotUnitary(a[ilda:ilda+i], x)
			}
			return
		}
		ix := kx + (n-1)*incX
		for i := n - 1; i >= 0; i-- {
			ilda := i * lda
			var tmp float64
			if nonUnit {
				tmp = a[ilda+i] * x[ix]
			} else {
				tmp = x[ix]
			}
			x[ix] = tmp + f64.DotInc(x, a[ilda:ilda+i], uintptr(i), uintptr(incX), 1, uintptr(kx), 0)
			ix -= incX
		}
		return
	}
	// Cases where a is transposed.
	if ul == blas.Upper {
		if incX == 1 {
			for i := n - 1; i >= 0; i-- {
				ilda := i * lda
				xi := x[i]
				f64.AxpyUnitary(xi, a[ilda+i+1:ilda+n], x[i+1:n])
				if nonUnit {
					x[i] *= a[ilda+i]
				}
			}
			return
		}
		ix := kx + (n-1)*incX
		for i := n - 1; i >= 0; i-- {
			ilda := i * lda
			xi := x[ix]
			f64.AxpyInc(xi, a[ilda+i+1:ilda+n], x, uintptr(n-i-1), 1, uintptr(incX), 0, uintptr(kx+(i+1)*incX))
			if nonUnit {
				x[ix] *= a[ilda+i]
			}
			ix -= incX
		}
		return
	}
	if incX == 1 {
		for i := 0; i < n; i++ {
			ilda := i * lda
			xi := x[i]
			f64.AxpyUnitary(xi, a[ilda:ilda+i], x)
			if nonUnit {
				x[i] *= a[i*lda+i]
			}
		}
		return
	}
	ix := kx
	for i := 0; i < n; i++ {
		ilda := i * lda
		xi := x[ix]
		f64.AxpyInc(xi, a[ilda:ilda+i], x, uintptr(i), 1, uintptr(incX), 0, uintptr(kx))
		if nonUnit {
			x[ix] *= a[ilda+i]
		}
		ix += incX
	}
}

// Dtrsv solves
//  A * x = b if tA == blas.NoTrans
//  A^T * x = b if tA == blas.Trans or blas.ConjTrans
// A is an n×n triangular matrix and x is a vector.
// At entry to the function, x contains the values of b, and the result 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.
func (Implementation) Dtrsv(ul blas.Uplo, tA blas.Transpose, d blas.Diag, n int, a []float64, lda int, x []float64, incX int) {
	// Test the input parameters
	// Verify inputs
	if ul != blas.Lower && ul != blas.Upper {
		panic(badUplo)
	}
	if tA != blas.NoTrans && tA != blas.Trans && tA != blas.ConjTrans {
		panic(badTranspose)
	}
	if d != blas.NonUnit && d != blas.Unit {
		panic(badDiag)
	}
	if n < 0 {
		panic(nLT0)
	}
	if lda*(n-1)+n > len(a) || lda < max(1, n) {
		panic(badLdA)
	}
	if incX == 0 {
		panic(zeroIncX)
	}
	if (incX > 0 && (n-1)*incX >= len(x)) || (incX < 0 && (1-n)*incX >= len(x)) {
		panic(badX)
	}
	// Quick return if possible
	if n == 0 {
		return
	}
	if n == 1 {
		if d == blas.NonUnit {
			x[0] /= a[0]
		}
		return
	}

	var kx int
	if incX < 0 {
		kx = -(n - 1) * incX
	}
	nonUnit := d == blas.NonUnit
	if tA == blas.NoTrans {
		if ul == blas.Upper {
			if incX == 1 {
				for i := n - 1; i >= 0; i-- {
					var sum float64
					atmp := a[i*lda+i+1 : i*lda+n]
					for j, v := range atmp {
						jv := i + j + 1
						sum += x[jv] * v
					}
					x[i] -= sum
					if nonUnit {
						x[i] /= a[i*lda+i]
					}
				}
				return
			}
			ix := kx + (n-1)*incX
			for i := n - 1; i >= 0; i-- {
				var sum float64
				jx := ix + incX
				atmp := a[i*lda+i+1 : i*lda+n]
				for _, v := range atmp {
					sum += x[jx] * v
					jx += incX
				}
				x[ix] -= sum
				if nonUnit {
					x[ix] /= a[i*lda+i]
				}
				ix -= incX
			}
			return
		}
		if incX == 1 {
			for i := 0; i < n; i++ {
				var sum float64
				atmp := a[i*lda : i*lda+i]
				for j, v := range atmp {
					sum += x[j] * v
				}
				x[i] -= sum
				if nonUnit {
					x[i] /= a[i*lda+i]
				}
			}
			return
		}
		ix := kx
		for i := 0; i < n; i++ {
			jx := kx
			var sum float64
			atmp := a[i*lda : i*lda+i]
			for _, v := range atmp {
				sum += x[jx] * v
				jx += incX
			}
			x[ix] -= sum
			if nonUnit {
				x[ix] /= a[i*lda+i]
			}
			ix += incX
		}
		return
	}
	// Cases where a is transposed.
	if ul == blas.Upper {
		if incX == 1 {
			for i := 0; i < n; i++ {
				if nonUnit {
					x[i] /= a[i*lda+i]
				}
				xi := x[i]
				atmp := a[i*lda+i+1 : i*lda+n]
				for j, v := range atmp {
					jv := j + i + 1
					x[jv] -= v * xi
				}
			}
			return
		}
		ix := kx
		for i := 0; i < n; i++ {
			if nonUnit {
				x[ix] /= a[i*lda+i]
			}
			xi := x[ix]
			jx := kx + (i+1)*incX
			atmp := a[i*lda+i+1 : i*lda+n]
			for _, v := range atmp {
				x[jx] -= v * xi
				jx += incX
			}
			ix += incX
		}
		return
	}
	if incX == 1 {
		for i := n - 1; i >= 0; i-- {
			if nonUnit {
				x[i] /= a[i*lda+i]
			}
			xi := x[i]
			atmp := a[i*lda : i*lda+i]
			for j, v := range atmp {
				x[j] -= v * xi
			}
		}
		return
	}
	ix := kx + (n-1)*incX
	for i := n - 1; i >= 0; i-- {
		if nonUnit {
			x[ix] /= a[i*lda+i]
		}
		xi := x[ix]
		jx := kx
		atmp := a[i*lda : i*lda+i]
		for _, v := range atmp {
			x[jx] -= v * xi
			jx += incX
		}
		ix -= incX
	}
}

// Dsymv computes
//    y = alpha * A * x + beta * y,
// where a is an n×n symmetric matrix, x and y are vectors, and alpha and
// beta are scalars.
func (Implementation) Dsymv(ul blas.Uplo, n int, alpha float64, a []float64, lda int, x []float64, incX int, beta float64, y []float64, incY int) {
	// Check inputs
	if ul != blas.Lower && ul != blas.Upper {
		panic(badUplo)
	}
	if n < 0 {
		panic(negativeN)
	}
	if lda > 1 && lda < n {
		panic(badLdA)
	}
	if incX == 0 {
		panic(zeroIncX)
	}
	if incY == 0 {
		panic(zeroIncY)
	}
	if (incX > 0 && (n-1)*incX >= len(x)) || (incX < 0 && (1-n)*incX >= len(x)) {
		panic(badX)
	}
	if (incY > 0 && (n-1)*incY >= len(y)) || (incY < 0 && (1-n)*incY >= len(y)) {
		panic(badY)
	}
	if lda*(n-1)+n > len(a) || lda < max(1, n) {
		panic(badLdA)
	}
	// Quick return if possible
	if n == 0 || (alpha == 0 && beta == 1) {
		return
	}

	// Set up start points
	var kx, ky int
	if incX > 0 {
		kx = 0
	} else {
		kx = -(n - 1) * incX
	}
	if incY > 0 {
		ky = 0
	} else {
		ky = -(n - 1) * incY
	}

	// Form y = beta * y
	if beta != 1 {
		if incY > 0 {
			Implementation{}.Dscal(n, beta, y, incY)
		} else {
			Implementation{}.Dscal(n, beta, y, -incY)
		}
	}

	if alpha == 0 {
		return
	}

	if n == 1 {
		y[0] += alpha * a[0] * x[0]
		return
	}

	if ul == blas.Upper {
		if incX == 1 {
			iy := ky
			for i := 0; i < n; i++ {
				xv := x[i] * alpha
				sum := x[i] * a[i*lda+i]
				jy := ky + (i+1)*incY
				atmp := a[i*lda+i+1 : i*lda+n]
				for j, v := range atmp {
					jp := j + i + 1
					sum += x[jp] * v
					y[jy] += xv * v
					jy += incY
				}
				y[iy] += alpha * sum
				iy += incY
			}
			return
		}
		ix := kx
		iy := ky
		for i := 0; i < n; i++ {
			xv := x[ix] * alpha
			sum := x[ix] * a[i*lda+i]
			jx := kx + (i+1)*incX
			jy := ky + (i+1)*incY
			atmp := a[i*lda+i+1 : i*lda+n]
			for _, v := range atmp {
				sum += x[jx] * v
				y[jy] += xv * v
				jx += incX
				jy += incY
			}
			y[iy] += alpha * sum
			ix += incX
			iy += incY
		}
		return
	}
	// Cases where a is lower triangular.
	if incX == 1 {
		iy := ky
		for i := 0; i < n; i++ {
			jy := ky
			xv := alpha * x[i]
			atmp := a[i*lda : i*lda+i]
			var sum float64
			for j, v := range atmp {
				sum += x[j] * v
				y[jy] += xv * v
				jy += incY
			}
			sum += x[i] * a[i*lda+i]
			sum *= alpha
			y[iy] += sum
			iy += incY
		}
		return
	}
	ix := kx
	iy := ky
	for i := 0; i < n; i++ {
		jx := kx
		jy := ky
		xv := alpha * x[ix]
		atmp := a[i*lda : i*lda+i]
		var sum float64
		for _, v := range atmp {
			sum += x[jx] * v
			y[jy] += xv * v
			jx += incX
			jy += incY
		}
		sum += x[ix] * a[i*lda+i]
		sum *= alpha
		y[iy] += sum
		ix += incX
		iy += incY
	}
}

// Dtbmv computes
//  x = A * x if tA == blas.NoTrans
//  x = A^T * x if tA == blas.Trans or blas.ConjTrans
// where A is an n×n triangular banded matrix with k diagonals, and x is a vector.
func (Implementation) Dtbmv(ul blas.Uplo, tA blas.Transpose, d blas.Diag, n, k int, a []float64, lda int, x []float64, incX int) {
	if ul != blas.Lower && ul != blas.Upper {
		panic(badUplo)
	}
	if tA != blas.NoTrans && tA != blas.Trans && tA != blas.ConjTrans {
		panic(badTranspose)
	}
	if d != blas.NonUnit && d != blas.Unit {
		panic(badDiag)
	}
	if n < 0 {
		panic(nLT0)
	}
	if k < 0 {
		panic(kLT0)
	}
	if lda*(n-1)+k+1 > len(a) || lda < k+1 {
		panic(badLdA)
	}
	if incX == 0 {
		panic(zeroIncX)
	}
	if (incX > 0 && (n-1)*incX >= len(x)) || (incX < 0 && (1-n)*incX >= len(x)) {
		panic(badX)
	}
	if n == 0 {
		return
	}
	var kx int
	if incX <= 0 {
		kx = -(n - 1) * incX
	} else if incX != 1 {
		kx = 0
	}

	nonunit := d != blas.Unit

	if tA == blas.NoTrans {
		if ul == blas.Upper {
			if incX == 1 {
				for i := 0; i < n; i++ {
					u := min(1+k, n-i)
					var sum float64
					atmp := a[i*lda:]
					xtmp := x[i:]
					for j := 1; j < u; j++ {
						sum += xtmp[j] * atmp[j]
					}
					if nonunit {
						sum += xtmp[0] * atmp[0]
					} else {
						sum += xtmp[0]
					}
					x[i] = sum
				}
				return
			}
			ix := kx
			for i := 0; i < n; i++ {
				u := min(1+k, n-i)
				var sum float64
				atmp := a[i*lda:]
				jx := incX
				for j := 1; j < u; j++ {
					sum += x[ix+jx] * atmp[j]
					jx += incX
				}
				if nonunit {
					sum += x[ix] * atmp[0]
				} else {
					sum += x[ix]
				}
				x[ix] = sum
				ix += incX
			}
			return
		}
		if incX == 1 {
			for i := n - 1; i >= 0; i-- {
				l := max(0, k-i)
				atmp := a[i*lda:]
				var sum float64
				for j := l; j < k; j++ {
					sum += x[i-k+j] * atmp[j]
				}
				if nonunit {
					sum += x[i] * atmp[k]
				} else {
					sum += x[i]
				}
				x[i] = sum
			}
			return
		}
		ix := kx + (n-1)*incX
		for i := n - 1; i >= 0; i-- {
			l := max(0, k-i)
			atmp := a[i*lda:]
			var sum float64
			jx := l * incX
			for j := l; j < k; j++ {
				sum += x[ix-k*incX+jx] * atmp[j]
				jx += incX
			}
			if nonunit {
				sum += x[ix] * atmp[k]
			} else {
				sum += x[ix]
			}
			x[ix] = sum
			ix -= incX
		}
		return
	}
	if ul == blas.Upper {
		if incX == 1 {
			for i := n - 1; i >= 0; i-- {
				u := k + 1
				if i < u {
					u = i + 1
				}
				var sum float64
				for j := 1; j < u; j++ {
					sum += x[i-j] * a[(i-j)*lda+j]
				}
				if nonunit {
					sum += x[i] * a[i*lda]
				} else {
					sum += x[i]
				}
				x[i] = sum
			}
			return
		}
		ix := kx + (n-1)*incX
		for i := n - 1; i >= 0; i-- {
			u := k + 1
			if i < u {
				u = i + 1
			}
			var sum float64
			jx := incX
			for j := 1; j < u; j++ {
				sum += x[ix-jx] * a[(i-j)*lda+j]
				jx += incX
			}
			if nonunit {
				sum += x[ix] * a[i*lda]
			} else {
				sum += x[ix]
			}
			x[ix] = sum
			ix -= incX
		}
		return
	}
	if incX == 1 {
		for i := 0; i < n; i++ {
			u := k
			if i+k >= n {
				u = n - i - 1
			}
			var sum float64
			for j := 0; j < u; j++ {
				sum += x[i+j+1] * a[(i+j+1)*lda+k-j-1]
			}
			if nonunit {
				sum += x[i] * a[i*lda+k]
			} else {
				sum += x[i]
			}
			x[i] = sum
		}
		return
	}
	ix := kx
	for i := 0; i < n; i++ {
		u := k
		if i+k >= n {
			u = n - i - 1
		}
		var (
			sum float64
			jx  int
		)
		for j := 0; j < u; j++ {
			sum += x[ix+jx+incX] * a[(i+j+1)*lda+k-j-1]
			jx += incX
		}
		if nonunit {
			sum += x[ix] * a[i*lda+k]
		} else {
			sum += x[ix]
		}
		x[ix] = sum
		ix += incX
	}
}

// Dtpmv computes
//  x = A * x if tA == blas.NoTrans
//  x = A^T * x if tA == blas.Trans or blas.ConjTrans
// where A is an n×n unit triangular matrix in packed format, and x is a vector.
func (Implementation) Dtpmv(ul blas.Uplo, tA blas.Transpose, d blas.Diag, n int, ap []float64, x []float64, incX int) {
	// Verify inputs
	if ul != blas.Lower && ul != blas.Upper {
		panic(badUplo)
	}
	if tA != blas.NoTrans && tA != blas.Trans && tA != blas.ConjTrans {
		panic(badTranspose)
	}
	if d != blas.NonUnit && d != blas.Unit {
		panic(badDiag)
	}
	if n < 0 {
		panic(nLT0)
	}
	if len(ap) < (n*(n+1))/2 {
		panic(badLdA)
	}
	if incX == 0 {
		panic(zeroIncX)
	}
	if (incX > 0 && (n-1)*incX >= len(x)) || (incX < 0 && (1-n)*incX >= len(x)) {
		panic(badX)
	}
	if n == 0 {
		return
	}
	var kx int
	if incX <= 0 {
		kx = -(n - 1) * incX
	}

	nonUnit := d == blas.NonUnit
	var offset int // Offset is the index of (i,i)
	if tA == blas.NoTrans {
		if ul == blas.Upper {
			if incX == 1 {
				for i := 0; i < n; i++ {
					xi := x[i]
					if nonUnit {
						xi *= ap[offset]
					}
					atmp := ap[offset+1 : offset+n-i]
					xtmp := x[i+1:]
					for j, v := range atmp {
						xi += v * xtmp[j]
					}
					x[i] = xi
					offset += n - i
				}
				return
			}
			ix := kx
			for i := 0; i < n; i++ {
				xix := x[ix]
				if nonUnit {
					xix *= ap[offset]
				}
				atmp := ap[offset+1 : offset+n-i]
				jx := kx + (i+1)*incX
				for _, v := range atmp {
					xix += v * x[jx]
					jx += incX
				}
				x[ix] = xix
				offset += n - i
				ix += incX
			}
			return
		}
		if incX == 1 {
			offset = n*(n+1)/2 - 1
			for i := n - 1; i >= 0; i-- {
				xi := x[i]
				if nonUnit {
					xi *= ap[offset]
				}
				atmp := ap[offset-i : offset]
				for j, v := range atmp {
					xi += v * x[j]
				}
				x[i] = xi
				offset -= i + 1
			}
			return
		}
		ix := kx + (n-1)*incX
		offset = n*(n+1)/2 - 1
		for i := n - 1; i >= 0; i-- {
			xix := x[ix]
			if nonUnit {
				xix *= ap[offset]
			}
			atmp := ap[offset-i : offset]
			jx := kx
			for _, v := range atmp {
				xix += v * x[jx]
				jx += incX
			}
			x[ix] = xix
			offset -= i + 1
			ix -= incX
		}
		return
	}
	// Cases where ap is transposed.
	if ul == blas.Upper {
		if incX == 1 {
			offset = n*(n+1)/2 - 1
			for i := n - 1; i >= 0; i-- {
				xi := x[i]
				atmp := ap[offset+1 : offset+n-i]
				xtmp := x[i+1:]
				for j, v := range atmp {
					xtmp[j] += v * xi
				}
				if nonUnit {
					x[i] *= ap[offset]
				}
				offset -= n - i + 1
			}
			return
		}
		ix := kx + (n-1)*incX
		offset = n*(n+1)/2 - 1
		for i := n - 1; i >= 0; i-- {
			xix := x[ix]
			jx := kx + (i+1)*incX
			atmp := ap[offset+1 : offset+n-i]
			for _, v := range atmp {
				x[jx] += v * xix
				jx += incX
			}
			if nonUnit {
				x[ix] *= ap[offset]
			}
			offset -= n - i + 1
			ix -= incX
		}
		return
	}
	if incX == 1 {
		for i := 0; i < n; i++ {
			xi := x[i]
			atmp := ap[offset-i : offset]
			for j, v := range atmp {
				x[j] += v * xi
			}
			if nonUnit {
				x[i] *= ap[offset]
			}
			offset += i + 2
		}
		return
	}
	ix := kx
	for i := 0; i < n; i++ {
		xix := x[ix]
		jx := kx
		atmp := ap[offset-i : offset]
		for _, v := range atmp {
			x[jx] += v * xix
			jx += incX
		}
		if nonUnit {
			x[ix] *= ap[offset]
		}
		ix += incX
		offset += i + 2
	}
}

// Dtbsv solves
//  A * x = b
// where A is an n×n triangular banded matrix with k diagonals in packed format,
// and x is a vector.
// At entry to the function, x contains the values of b, and the result 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.
func (Implementation) Dtbsv(ul blas.Uplo, tA blas.Transpose, d blas.Diag, n, k int, a []float64, lda int, x []float64, incX int) {
	if ul != blas.Lower && ul != blas.Upper {
		panic(badUplo)
	}
	if tA != blas.NoTrans && tA != blas.Trans && tA != blas.ConjTrans {
		panic(badTranspose)
	}
	if d != blas.NonUnit && d != blas.Unit {
		panic(badDiag)
	}
	if n < 0 {
		panic(nLT0)
	}
	if lda*(n-1)+k+1 > len(a) || lda < k+1 {
		panic(badLdA)
	}
	if incX == 0 {
		panic(zeroIncX)
	}
	if (incX > 0 && (n-1)*incX >= len(x)) || (incX < 0 && (1-n)*incX >= len(x)) {
		panic(badX)
	}
	if n == 0 {
		return
	}
	var kx int
	if incX < 0 {
		kx = -(n - 1) * incX
	} else {
		kx = 0
	}
	nonUnit := d == blas.NonUnit
	// Form x = A^-1 x.
	// Several cases below use subslices for speed improvement.
	// The incX != 1 cases usually do not because incX may be negative.
	if tA == blas.NoTrans {
		if ul == blas.Upper {
			if incX == 1 {
				for i := n - 1; i >= 0; i-- {
					bands := k
					if i+bands >= n {
						bands = n - i - 1
					}
					atmp := a[i*lda+1:]
					xtmp := x[i+1 : i+bands+1]
					var sum float64
					for j, v := range xtmp {
						sum += v * atmp[j]
					}
					x[i] -= sum
					if nonUnit {
						x[i] /= a[i*lda]
					}
				}
				return
			}
			ix := kx + (n-1)*incX
			for i := n - 1; i >= 0; i-- {
				max := k + 1
				if i+max > n {
					max = n - i
				}
				atmp := a[i*lda:]
				var (
					jx  int
					sum float64
				)
				for j := 1; j < max; j++ {
					jx += incX
					sum += x[ix+jx] * atmp[j]
				}
				x[ix] -= sum
				if nonUnit {
					x[ix] /= atmp[0]
				}
				ix -= incX
			}
			return
		}
		if incX == 1 {
			for i := 0; i < n; i++ {
				bands := k
				if i-k < 0 {
					bands = i
				}
				atmp := a[i*lda+k-bands:]
				xtmp := x[i-bands : i]
				var sum float64
				for j, v := range xtmp {
					sum += v * atmp[j]
				}
				x[i] -= sum
				if nonUnit {
					x[i] /= atmp[bands]
				}
			}
			return
		}
		ix := kx
		for i := 0; i < n; i++ {
			bands := k
			if i-k < 0 {
				bands = i
			}
			atmp := a[i*lda+k-bands:]
			var (
				sum float64
				jx  int
			)
			for j := 0; j < bands; j++ {
				sum += x[ix-bands*incX+jx] * atmp[j]
				jx += incX
			}
			x[ix] -= sum
			if nonUnit {
				x[ix] /= atmp[bands]
			}
			ix += incX
		}
		return
	}
	// Cases where a is transposed.
	if ul == blas.Upper {
		if incX == 1 {
			for i := 0; i < n; i++ {
				bands := k
				if i-k < 0 {
					bands = i
				}
				var sum float64
				for j := 0; j < bands; j++ {
					sum += x[i-bands+j] * a[(i-bands+j)*lda+bands-j]
				}
				x[i] -= sum
				if nonUnit {
					x[i] /= a[i*lda]
				}
			}
			return
		}
		ix := kx
		for i := 0; i < n; i++ {
			bands := k
			if i-k < 0 {
				bands = i
			}
			var (
				sum float64
				jx  int
			)
			for j := 0; j < bands; j++ {
				sum += x[ix-bands*incX+jx] * a[(i-bands+j)*lda+bands-j]
				jx += incX
			}
			x[ix] -= sum
			if nonUnit {
				x[ix] /= a[i*lda]
			}
			ix += incX
		}
		return
	}
	if incX == 1 {
		for i := n - 1; i >= 0; i-- {
			bands := k
			if i+bands >= n {
				bands = n - i - 1
			}
			var sum float64
			xtmp := x[i+1 : i+1+bands]
			for j, v := range xtmp {
				sum += v * a[(i+j+1)*lda+k-j-1]
			}
			x[i] -= sum
			if nonUnit {
				x[i] /= a[i*lda+k]
			}
		}
		return
	}
	ix := kx + (n-1)*incX
	for i := n - 1; i >= 0; i-- {
		bands := k
		if i+bands >= n {
			bands = n - i - 1
		}
		var (
			sum float64
			jx  int
		)
		for j := 0; j < bands; j++ {
			sum += x[ix+jx+incX] * a[(i+j+1)*lda+k-j-1]
			jx += incX
		}
		x[ix] -= sum
		if nonUnit {
			x[ix] /= a[i*lda+k]
		}
		ix -= incX
	}
}

// Dsbmv performs
//  y = alpha * A * x + beta * y
// where A is an n×n symmetric banded matrix, x and y are vectors, and alpha
// and beta are scalars.
func (Implementation) Dsbmv(ul blas.Uplo, n, k int, alpha float64, a []float64, lda int, x []float64, incX int, beta float64, y []float64, incY int) {
	if ul != blas.Lower && ul != blas.Upper {
		panic(badUplo)
	}
	if n < 0 {
		panic(nLT0)
	}

	if incX == 0 {
		panic(zeroIncX)
	}
	if incY == 0 {
		panic(zeroIncY)
	}
	if (incX > 0 && (n-1)*incX >= len(x)) || (incX < 0 && (1-n)*incX >= len(x)) {
		panic(badX)
	}
	if (incY > 0 && (n-1)*incY >= len(y)) || (incY < 0 && (1-n)*incY >= len(y)) {
		panic(badY)
	}
	if lda*(n-1)+k+1 > len(a) || lda < k+1 {
		panic(badLdA)
	}

	// Quick return if possible
	if n == 0 || (alpha == 0 && beta == 1) {
		return
	}

	// Set up indexes
	lenX := n
	lenY := n
	var kx, ky int
	if incX > 0 {
		kx = 0
	} else {
		kx = -(lenX - 1) * incX
	}
	if incY > 0 {
		ky = 0
	} else {
		ky = -(lenY - 1) * incY
	}

	// First form y := beta * y
	if incY > 0 {
		Implementation{}.Dscal(lenY, beta, y, incY)
	} else {
		Implementation{}.Dscal(lenY, beta, y, -incY)
	}

	if alpha == 0 {
		return
	}

	if ul == blas.Upper {
		if incX == 1 {
			iy := ky
			for i := 0; i < n; i++ {
				atmp := a[i*lda:]
				tmp := alpha * x[i]
				sum := tmp * atmp[0]
				u := min(k, n-i-1)
				jy := incY
				for j := 1; j <= u; j++ {
					v := atmp[j]
					sum += alpha * x[i+j] * v
					y[iy+jy] += tmp * v
					jy += incY
				}
				y[iy] += sum
				iy += incY
			}
			return
		}
		ix := kx
		iy := ky
		for i := 0; i < n; i++ {
			atmp := a[i*lda:]
			tmp := alpha * x[ix]
			sum := tmp * atmp[0]
			u := min(k, n-i-1)
			jx := incX
			jy := incY
			for j := 1; j <= u; j++ {
				v := atmp[j]
				sum += alpha * x[ix+jx] * v
				y[iy+jy] += tmp * v
				jx += incX
				jy += incY
			}
			y[iy] += sum
			ix += incX
			iy += incY
		}
		return
	}

	// Casses where a has bands below the diagonal.
	if incX == 1 {
		iy := ky
		for i := 0; i < n; i++ {
			l := max(0, k-i)
			tmp := alpha * x[i]
			jy := l * incY
			atmp := a[i*lda:]
			for j := l; j < k; j++ {
				v := atmp[j]
				y[iy] += alpha * v * x[i-k+j]
				y[iy-k*incY+jy] += tmp * v
				jy += incY
			}
			y[iy] += tmp * atmp[k]
			iy += incY
		}
		return
	}
	ix := kx
	iy := ky
	for i := 0; i < n; i++ {
		l := max(0, k-i)
		tmp := alpha * x[ix]
		jx := l * incX
		jy := l * incY
		atmp := a[i*lda:]
		for j := l; j < k; j++ {
			v := atmp[j]
			y[iy] += alpha * v * x[ix-k*incX+jx]
			y[iy-k*incY+jy] += tmp * v
			jx += incX
			jy += incY
		}
		y[iy] += tmp * atmp[k]
		ix += incX
		iy += incY
	}
}

// Dsyr performs the rank-one update
//  a += alpha * x * x^T
// where a is an n×n symmetric matrix, and x is a vector.
func (Implementation) Dsyr(ul blas.Uplo, n int, alpha float64, x []float64, incX int, a []float64, lda int) {
	if ul != blas.Lower && ul != blas.Upper {
		panic(badUplo)
	}
	if n < 0 {
		panic(nLT0)
	}
	if incX == 0 {
		panic(zeroIncX)
	}
	if (incX > 0 && (n-1)*incX >= len(x)) || (incX < 0 && (1-n)*incX >= len(x)) {
		panic(badX)
	}
	if lda*(n-1)+n > len(a) || lda < max(1, n) {
		panic(badLdA)
	}
	if alpha == 0 || n == 0 {
		return
	}

	lenX := n
	var kx int
	if incX > 0 {
		kx = 0
	} else {
		kx = -(lenX - 1) * incX
	}
	if ul == blas.Upper {
		if incX == 1 {
			for i := 0; i < n; i++ {
				tmp := x[i] * alpha
				if tmp != 0 {
					atmp := a[i*lda+i : i*lda+n]
					xtmp := x[i:n]
					for j, v := range xtmp {
						atmp[j] += v * tmp
					}
				}
			}
			return
		}
		ix := kx
		for i := 0; i < n; i++ {
			tmp := x[ix] * alpha
			if tmp != 0 {
				jx := ix
				atmp := a[i*lda:]
				for j := i; j < n; j++ {
					atmp[j] += x[jx] * tmp
					jx += incX
				}
			}
			ix += incX
		}
		return
	}
	// Cases where a is lower triangular.
	if incX == 1 {
		for i := 0; i < n; i++ {
			tmp := x[i] * alpha
			if tmp != 0 {
				atmp := a[i*lda:]
				xtmp := x[:i+1]
				for j, v := range xtmp {
					atmp[j] += tmp * v
				}
			}
		}
		return
	}
	ix := kx
	for i := 0; i < n; i++ {
		tmp := x[ix] * alpha
		if tmp != 0 {
			atmp := a[i*lda:]
			jx := kx
			for j := 0; j < i+1; j++ {
				atmp[j] += tmp * x[jx]
				jx += incX
			}
		}
		ix += incX
	}
}

// Dsyr2 performs the symmetric rank-two update
//  A += alpha * x * y^T + alpha * y * x^T
// where A is a symmetric n×n matrix, x and y are vectors, and alpha is a scalar.
func (Implementation) Dsyr2(ul blas.Uplo, n int, alpha float64, x []float64, incX int, y []float64, incY int, a []float64, lda int) {
	if ul != blas.Lower && ul != blas.Upper {
		panic(badUplo)
	}
	if n < 0 {
		panic(nLT0)
	}
	if incX == 0 {
		panic(zeroIncX)
	}
	if incY == 0 {
		panic(zeroIncY)
	}
	if (incX > 0 && (n-1)*incX >= len(x)) || (incX < 0 && (1-n)*incX >= len(x)) {
		panic(badX)
	}
	if (incY > 0 && (n-1)*incY >= len(y)) || (incY < 0 && (1-n)*incY >= len(y)) {
		panic(badY)
	}
	if lda*(n-1)+n > len(a) || lda < max(1, n) {
		panic(badLdA)
	}
	if alpha == 0 {
		return
	}

	var ky, kx int
	if incY > 0 {
		ky = 0
	} else {
		ky = -(n - 1) * incY
	}
	if incX > 0 {
		kx = 0
	} else {
		kx = -(n - 1) * incX
	}
	if ul == blas.Upper {
		if incX == 1 && incY == 1 {
			for i := 0; i < n; i++ {
				xi := x[i]
				yi := y[i]
				atmp := a[i*lda:]
				for j := i; j < n; j++ {
					atmp[j] += alpha * (xi*y[j] + x[j]*yi)
				}
			}
			return
		}
		ix := kx
		iy := ky
		for i := 0; i < n; i++ {
			jx := kx + i*incX
			jy := ky + i*incY
			xi := x[ix]
			yi := y[iy]
			atmp := a[i*lda:]
			for j := i; j < n; j++ {
				atmp[j] += alpha * (xi*y[jy] + x[jx]*yi)
				jx += incX
				jy += incY
			}
			ix += incX
			iy += incY
		}
		return
	}
	if incX == 1 && incY == 1 {
		for i := 0; i < n; i++ {
			xi := x[i]
			yi := y[i]
			atmp := a[i*lda:]
			for j := 0; j <= i; j++ {
				atmp[j] += alpha * (xi*y[j] + x[j]*yi)
			}
		}
		return
	}
	ix := kx
	iy := ky
	for i := 0; i < n; i++ {
		jx := kx
		jy := ky
		xi := x[ix]
		yi := y[iy]
		atmp := a[i*lda:]
		for j := 0; j <= i; j++ {
			atmp[j] += alpha * (xi*y[jy] + x[jx]*yi)
			jx += incX
			jy += incY
		}
		ix += incX
		iy += incY
	}
}

// Dtpsv solves
//  A * x = b if tA == blas.NoTrans
//  A^T * x = b if tA == blas.Trans or blas.ConjTrans
// where A is an n×n triangular matrix in packed format and x is a vector.
// At entry to the function, x contains the values of b, and the result 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.
func (Implementation) Dtpsv(ul blas.Uplo, tA blas.Transpose, d blas.Diag, n int, ap []float64, x []float64, incX int) {
	// Verify inputs
	if ul != blas.Lower && ul != blas.Upper {
		panic(badUplo)
	}
	if tA != blas.NoTrans && tA != blas.Trans && tA != blas.ConjTrans {
		panic(badTranspose)
	}
	if d != blas.NonUnit && d != blas.Unit {
		panic(badDiag)
	}
	if n < 0 {
		panic(nLT0)
	}
	if len(ap) < (n*(n+1))/2 {
		panic(badLdA)
	}
	if incX == 0 {
		panic(zeroIncX)
	}
	if (incX > 0 && (n-1)*incX >= len(x)) || (incX < 0 && (1-n)*incX >= len(x)) {
		panic(badX)
	}
	if n == 0 {
		return
	}
	var kx int
	if incX <= 0 {
		kx = -(n - 1) * incX
	}

	nonUnit := d == blas.NonUnit
	var offset int // Offset is the index of (i,i)
	if tA == blas.NoTrans {
		if ul == blas.Upper {
			offset = n*(n+1)/2 - 1
			if incX == 1 {
				for i := n - 1; i >= 0; i-- {
					atmp := ap[offset+1 : offset+n-i]
					xtmp := x[i+1:]
					var sum float64
					for j, v := range atmp {
						sum += v * xtmp[j]
					}
					x[i] -= sum
					if nonUnit {
						x[i] /= ap[offset]
					}
					offset -= n - i + 1
				}
				return
			}
			ix := kx + (n-1)*incX
			for i := n - 1; i >= 0; i-- {
				atmp := ap[offset+1 : offset+n-i]
				jx := kx + (i+1)*incX
				var sum float64
				for _, v := range atmp {
					sum += v * x[jx]
					jx += incX
				}
				x[ix] -= sum
				if nonUnit {
					x[ix] /= ap[offset]
				}
				ix -= incX
				offset -= n - i + 1
			}
			return
		}
		if incX == 1 {
			for i := 0; i < n; i++ {
				atmp := ap[offset-i : offset]
				var sum float64
				for j, v := range atmp {
					sum += v * x[j]
				}
				x[i] -= sum
				if nonUnit {
					x[i] /= ap[offset]
				}
				offset += i + 2
			}
			return
		}
		ix := kx
		for i := 0; i < n; i++ {
			jx := kx
			atmp := ap[offset-i : offset]
			var sum float64
			for _, v := range atmp {
				sum += v * x[jx]
				jx += incX
			}
			x[ix] -= sum
			if nonUnit {
				x[ix] /= ap[offset]
			}
			ix += incX
			offset += i + 2
		}
		return
	}
	// Cases where ap is transposed.
	if ul == blas.Upper {
		if incX == 1 {
			for i := 0; i < n; i++ {
				if nonUnit {
					x[i] /= ap[offset]
				}
				xi := x[i]
				atmp := ap[offset+1 : offset+n-i]
				xtmp := x[i+1:]
				for j, v := range atmp {
					xtmp[j] -= v * xi
				}
				offset += n - i
			}
			return
		}
		ix := kx
		for i := 0; i < n; i++ {
			if nonUnit {
				x[ix] /= ap[offset]
			}
			xix := x[ix]
			atmp := ap[offset+1 : offset+n-i]
			jx := kx + (i+1)*incX
			for _, v := range atmp {
				x[jx] -= v * xix
				jx += incX
			}
			ix += incX
			offset += n - i
		}
		return
	}
	if incX == 1 {
		offset = n*(n+1)/2 - 1
		for i := n - 1; i >= 0; i-- {
			if nonUnit {
				x[i] /= ap[offset]
			}
			xi := x[i]
			atmp := ap[offset-i : offset]
			for j, v := range atmp {
				x[j] -= v * xi
			}
			offset -= i + 1
		}
		return
	}
	ix := kx + (n-1)*incX
	offset = n*(n+1)/2 - 1
	for i := n - 1; i >= 0; i-- {
		if nonUnit {
			x[ix] /= ap[offset]
		}
		xix := x[ix]
		atmp := ap[offset-i : offset]
		jx := kx
		for _, v := range atmp {
			x[jx] -= v * xix
			jx += incX
		}
		ix -= incX
		offset -= i + 1
	}
}

// Dspmv performs
//    y = alpha * A * x + beta * y,
// where A is an n×n symmetric matrix in packed format, x and y are vectors
// and alpha and beta are scalars.
func (Implementation) Dspmv(ul blas.Uplo, n int, alpha float64, a []float64, x []float64, incX int, beta float64, y []float64, incY int) {
	// Verify inputs
	if ul != blas.Lower && ul != blas.Upper {
		panic(badUplo)
	}
	if n < 0 {
		panic(nLT0)
	}
	if len(a) < (n*(n+1))/2 {
		panic(badLdA)
	}
	if incX == 0 {
		panic(zeroIncX)
	}
	if incY == 0 {
		panic(zeroIncY)
	}
	if (incX > 0 && (n-1)*incX >= len(x)) || (incX < 0 && (1-n)*incX >= len(x)) {
		panic(badX)
	}
	if (incY > 0 && (n-1)*incY >= len(y)) || (incY < 0 && (1-n)*incY >= len(y)) {
		panic(badY)
	}
	// Quick return if possible
	if n == 0 || (alpha == 0 && beta == 1) {
		return
	}

	// Set up start points
	var kx, ky int
	if incX > 0 {
		kx = 0
	} else {
		kx = -(n - 1) * incX
	}
	if incY > 0 {
		ky = 0
	} else {
		ky = -(n - 1) * incY
	}

	// Form y = beta * y
	if beta != 1 {
		if incY > 0 {
			Implementation{}.Dscal(n, beta, y, incY)
		} else {
			Implementation{}.Dscal(n, beta, y, -incY)
		}
	}

	if alpha == 0 {
		return
	}

	if n == 1 {
		y[0] += alpha * a[0] * x[0]
		return
	}
	var offset int // Offset is the index of (i,i).
	if ul == blas.Upper {
		if incX == 1 {
			iy := ky
			for i := 0; i < n; i++ {
				xv := x[i] * alpha
				sum := a[offset] * x[i]
				atmp := a[offset+1 : offset+n-i]
				xtmp := x[i+1:]
				jy := ky + (i+1)*incY
				for j, v := range atmp {
					sum += v * xtmp[j]
					y[jy] += v * xv
					jy += incY
				}
				y[iy] += alpha * sum
				iy += incY
				offset += n - i
			}
			return
		}
		ix := kx
		iy := ky
		for i := 0; i < n; i++ {
			xv := x[ix] * alpha
			sum := a[offset] * x[ix]
			atmp := a[offset+1 : offset+n-i]
			jx := kx + (i+1)*incX
			jy := ky + (i+1)*incY
			for _, v := range atmp {
				sum += v * x[jx]
				y[jy] += v * xv
				jx += incX
				jy += incY
			}
			y[iy] += alpha * sum
			ix += incX
			iy += incY
			offset += n - i
		}
		return
	}
	if incX == 1 {
		iy := ky
		for i := 0; i < n; i++ {
			xv := x[i] * alpha
			atmp := a[offset-i : offset]
			jy := ky
			var sum float64
			for j, v := range atmp {
				sum += v * x[j]
				y[jy] += v * xv
				jy += incY
			}
			sum += a[offset] * x[i]
			y[iy] += alpha * sum
			iy += incY
			offset += i + 2
		}
		return
	}
	ix := kx
	iy := ky
	for i := 0; i < n; i++ {
		xv := x[ix] * alpha
		atmp := a[offset-i : offset]
		jx := kx
		jy := ky
		var sum float64
		for _, v := range atmp {
			sum += v * x[jx]
			y[jy] += v * xv
			jx += incX
			jy += incY
		}

		sum += a[offset] * x[ix]
		y[iy] += alpha * sum
		ix += incX
		iy += incY
		offset += i + 2
	}
}

// Dspr computes the rank-one operation
//  a += alpha * x * x^T
// where a is an n×n symmetric matrix in packed format, x is a vector, and
// alpha is a scalar.
func (Implementation) Dspr(ul blas.Uplo, n int, alpha float64, x []float64, incX int, a []float64) {
	if ul != blas.Lower && ul != blas.Upper {
		panic(badUplo)
	}
	if n < 0 {
		panic(nLT0)
	}
	if incX == 0 {
		panic(zeroIncX)
	}
	if (incX > 0 && (n-1)*incX >= len(x)) || (incX < 0 && (1-n)*incX >= len(x)) {
		panic(badX)
	}
	if len(a) < (n*(n+1))/2 {
		panic(badLdA)
	}
	if alpha == 0 || n == 0 {
		return
	}
	lenX := n
	var kx int
	if incX > 0 {
		kx = 0
	} else {
		kx = -(lenX - 1) * incX
	}
	var offset int // Offset is the index of (i,i).
	if ul == blas.Upper {
		if incX == 1 {
			for i := 0; i < n; i++ {
				atmp := a[offset:]
				xv := alpha * x[i]
				xtmp := x[i:n]
				for j, v := range xtmp {
					atmp[j] += xv * v
				}
				offset += n - i
			}
			return
		}
		ix := kx
		for i := 0; i < n; i++ {
			jx := kx + i*incX
			atmp := a[offset:]
			xv := alpha * x[ix]
			for j := 0; j < n-i; j++ {
				atmp[j] += xv * x[jx]
				jx += incX
			}
			ix += incX
			offset += n - i
		}
		return
	}
	if incX == 1 {
		for i := 0; i < n; i++ {
			atmp := a[offset-i:]
			xv := alpha * x[i]
			xtmp := x[:i+1]
			for j, v := range xtmp {
				atmp[j] += xv * v
			}
			offset += i + 2
		}
		return
	}
	ix := kx
	for i := 0; i < n; i++ {
		jx := kx
		atmp := a[offset-i:]
		xv := alpha * x[ix]
		for j := 0; j <= i; j++ {
			atmp[j] += xv * x[jx]
			jx += incX
		}
		ix += incX
		offset += i + 2
	}
}

// Dspr2 performs the symmetric rank-2 update
//  A += alpha * x * y^T + alpha * y * x^T,
// where A is an n×n symmetric matrix in packed format, x and y are vectors,
// and alpha is a scalar.
func (Implementation) Dspr2(ul blas.Uplo, n int, alpha float64, x []float64, incX int, y []float64, incY int, ap []float64) {
	if ul != blas.Lower && ul != blas.Upper {
		panic(badUplo)
	}
	if n < 0 {
		panic(nLT0)
	}
	if incX == 0 {
		panic(zeroIncX)
	}
	if incY == 0 {
		panic(zeroIncY)
	}
	if (incX > 0 && (n-1)*incX >= len(x)) || (incX < 0 && (1-n)*incX >= len(x)) {
		panic(badX)
	}
	if (incY > 0 && (n-1)*incY >= len(y)) || (incY < 0 && (1-n)*incY >= len(y)) {
		panic(badY)
	}
	if len(ap) < (n*(n+1))/2 {
		panic(badLdA)
	}
	if alpha == 0 {
		return
	}
	var ky, kx int
	if incY > 0 {
		ky = 0
	} else {
		ky = -(n - 1) * incY
	}
	if incX > 0 {
		kx = 0
	} else {
		kx = -(n - 1) * incX
	}
	var offset int // Offset is the index of (i,i).
	if ul == blas.Upper {
		if incX == 1 && incY == 1 {
			for i := 0; i < n; i++ {
				atmp := ap[offset:]
				xi := x[i]
				yi := y[i]
				xtmp := x[i:n]
				ytmp := y[i:n]
				for j, v := range xtmp {
					atmp[j] += alpha * (xi*ytmp[j] + v*yi)
				}
				offset += n - i
			}
			return
		}
		ix := kx
		iy := ky
		for i := 0; i < n; i++ {
			jx := kx + i*incX
			jy := ky + i*incY
			atmp := ap[offset:]
			xi := x[ix]
			yi := y[iy]
			for j := 0; j < n-i; j++ {
				atmp[j] += alpha * (xi*y[jy] + x[jx]*yi)
				jx += incX
				jy += incY
			}
			ix += incX
			iy += incY
			offset += n - i
		}
		return
	}
	if incX == 1 && incY == 1 {
		for i := 0; i < n; i++ {
			atmp := ap[offset-i:]
			xi := x[i]
			yi := y[i]
			xtmp := x[:i+1]
			for j, v := range xtmp {
				atmp[j] += alpha * (xi*y[j] + v*yi)
			}
			offset += i + 2
		}
		return
	}
	ix := kx
	iy := ky
	for i := 0; i < n; i++ {
		jx := kx
		jy := ky
		atmp := ap[offset-i:]
		for j := 0; j <= i; j++ {
			atmp[j] += alpha * (x[ix]*y[jy] + x[jx]*y[iy])
			jx += incX
			jy += incY
		}
		ix += incX
		iy += incY
		offset += i + 2
	}
}