--- /dev/null
+// 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 (
+ "math/cmplx"
+
+ "gonum.org/v1/gonum/blas"
+ "gonum.org/v1/gonum/internal/asm/c128"
+)
+
+// Zgemv 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.
+func (Implementation) Zgemv(trans blas.Transpose, m, n int, alpha complex128, a []complex128, lda int, x []complex128, incX int, beta complex128, y []complex128, incY int) {
+ checkZMatrix('A', m, n, a, lda)
+ switch trans {
+ default:
+ panic(badTranspose)
+ case blas.NoTrans:
+ checkZVector('x', n, x, incX)
+ checkZVector('y', m, y, incY)
+ case blas.Trans, blas.ConjTrans:
+ checkZVector('x', m, x, incX)
+ checkZVector('y', n, y, incY)
+ }
+
+ if m == 0 || n == 0 || (alpha == 0 && beta == 1) {
+ return
+ }
+
+ var lenX, lenY int
+ if trans == blas.NoTrans {
+ lenX = n
+ lenY = m
+ } else {
+ lenX = m
+ lenY = n
+ }
+ 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 {
+ c128.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 {
+ c128.ScalInc(beta, y, uintptr(lenY), uintptr(incY))
+ } else {
+ c128.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 * c128.DotuUnitary(a[i*lda:i*lda+n], x[:n])
+ iy += incY
+ }
+ return
+ }
+ for i := 0; i < m; i++ {
+ y[iy] += alpha * c128.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++ {
+ c128.AxpyUnitary(alpha*x[ix], a[i*lda:i*lda+n], y[:n])
+ ix += incX
+ }
+ return
+ }
+ for i := 0; i < m; i++ {
+ c128.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
+ }
+}
+
+// Zgerc 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.
+func (Implementation) Zgerc(m, n int, alpha complex128, x []complex128, incX int, y []complex128, incY int, a []complex128, lda int) {
+ checkZMatrix('A', m, n, a, lda)
+ checkZVector('x', m, x, incX)
+ checkZVector('y', n, y, incY)
+
+ if m == 0 || n == 0 || 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])
+ c128.AxpyInc(tmp, x, a[j:], uintptr(m), uintptr(incX), uintptr(lda), uintptr(kx), 0)
+ }
+ jy += incY
+ }
+}
+
+// Zgeru 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.
+func (Implementation) Zgeru(m, n int, alpha complex128, x []complex128, incX int, y []complex128, incY int, a []complex128, lda int) {
+ checkZMatrix('A', m, n, a, lda)
+ checkZVector('x', m, x, incX)
+ checkZVector('y', n, y, incY)
+
+ if m == 0 || n == 0 || 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]
+ c128.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]
+ c128.AxpyInc(tmp, y, a[i*lda:i*lda+n], uintptr(n), uintptr(incY), 1, uintptr(jy), 0)
+ }
+ kx += incX
+ }
+}
+
+// Zhemv 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.
+func (Implementation) Zhemv(uplo blas.Uplo, n int, alpha complex128, a []complex128, lda int, x []complex128, incX int, beta complex128, y []complex128, incY int) {
+ if uplo != blas.Upper && uplo != blas.Lower {
+ panic(badUplo)
+ }
+ checkZMatrix('A', n, n, a, lda)
+ checkZVector('x', n, x, incX)
+ checkZVector('y', n, y, incY)
+
+ if n == 0 || (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 complex128
+ 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 complex128
+ 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 complex128
+ 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 complex128
+ 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
+ }
+ }
+}
+
+// Zher 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.
+func (Implementation) Zher(uplo blas.Uplo, n int, alpha float64, x []complex128, incX int, a []complex128, lda int) {
+ if uplo != blas.Upper && uplo != blas.Lower {
+ panic(badUplo)
+ }
+ checkZMatrix('A', n, n, a, lda)
+ checkZVector('x', n, x, incX)
+
+ if n == 0 || 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
+ }
+}
+
+// Zher2 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.
+func (Implementation) Zher2(uplo blas.Uplo, n int, alpha complex128, x []complex128, incX int, y []complex128, incY int, a []complex128, lda int) {
+ if uplo != blas.Upper && uplo != blas.Lower {
+ panic(badUplo)
+ }
+ checkZMatrix('A', n, n, a, lda)
+ checkZVector('x', n, x, incX)
+ checkZVector('y', n, y, incY)
+
+ if n == 0 || 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
+ }
+}
+
+// Zhpmv 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.
+func (Implementation) Zhpmv(uplo blas.Uplo, n int, alpha complex128, ap []complex128, x []complex128, incX int, beta complex128, y []complex128, incY int) {
+ if uplo != blas.Upper && uplo != blas.Lower {
+ panic(badUplo)
+ }
+ checkZVector('x', n, x, incX)
+ checkZVector('y', n, y, incY)
+ if len(ap) < n*(n+1)/2 {
+ panic("blas: insufficient A packed matrix slice length")
+ }
+
+ if n == 0 || (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 complex128
+ 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 complex128
+ 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 complex128
+ 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 complex128
+ 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
+ }
+ }
+}
+
+// Zhpr 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.
+func (Implementation) Zhpr(uplo blas.Uplo, n int, alpha float64, x []complex128, incX int, ap []complex128) {
+ if uplo != blas.Upper && uplo != blas.Lower {
+ panic(badUplo)
+ }
+ if n < 0 {
+ panic(nLT0)
+ }
+ checkZVector('x', n, x, incX)
+ if len(ap) < n*(n+1)/2 {
+ panic("blas: insufficient A packed matrix slice length")
+ }
+
+ if n == 0 || 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
+ }
+ }
+}
+
+// Zhpr2 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.
+func (Implementation) Zhpr2(uplo blas.Uplo, n int, alpha complex128, x []complex128, incX int, y []complex128, incY int, ap []complex128) {
+ if uplo != blas.Upper && uplo != blas.Lower {
+ panic(badUplo)
+ }
+ if n < 0 {
+ panic(nLT0)
+ }
+ checkZVector('x', n, x, incX)
+ checkZVector('y', n, y, incY)
+ if len(ap) < n*(n+1)/2 {
+ panic("blas: insufficient A packed matrix slice length")
+ }
+
+ if n == 0 || 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
+ }
+ }
+}
+
+// Ztpmv 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.
+func (Implementation) Ztpmv(uplo blas.Uplo, trans blas.Transpose, diag blas.Diag, n int, ap []complex128, x []complex128, incX int) {
+ if uplo != blas.Upper && uplo != blas.Lower {
+ panic(badUplo)
+ }
+ if trans != blas.NoTrans && trans != blas.Trans && trans != blas.ConjTrans {
+ panic(badTranspose)
+ }
+ if diag != blas.Unit && diag != blas.NonUnit {
+ panic(badDiag)
+ }
+ checkZVector('x', n, x, incX)
+ if len(ap) < n*(n+1)/2 {
+ panic("blas: insufficient A packed matrix slice length")
+ }
+
+ if n == 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 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] += c128.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] += c128.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] += c128.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] += c128.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 {
+ c128.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 {
+ c128.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 {
+ c128.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 {
+ c128.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
+ }
+ }
+ }
+}
+
+// Ztrmv 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.
+func (Implementation) Ztrmv(uplo blas.Uplo, trans blas.Transpose, diag blas.Diag, n int, a []complex128, lda int, x []complex128, incX int) {
+ if uplo != blas.Upper && uplo != blas.Lower {
+ panic(badUplo)
+ }
+ if trans != blas.NoTrans && trans != blas.Trans && trans != blas.ConjTrans {
+ panic(badTranspose)
+ }
+ if diag != blas.Unit && diag != blas.NonUnit {
+ panic(badDiag)
+ }
+ checkZMatrix('A', n, n, a, lda)
+ checkZVector('x', n, x, incX)
+
+ if n == 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 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] += c128.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] += c128.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] += c128.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] += c128.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 {
+ c128.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 {
+ c128.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 {
+ c128.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 {
+ c128.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
+ }
+ }
+ }
+}
+
+// Ztrsv 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.
+func (Implementation) Ztrsv(uplo blas.Uplo, trans blas.Transpose, diag blas.Diag, n int, a []complex128, lda int, x []complex128, incX int) {
+ if uplo != blas.Upper && uplo != blas.Lower {
+ panic(badUplo)
+ }
+ if trans != blas.NoTrans && trans != blas.Trans && trans != blas.ConjTrans {
+ panic(badTranspose)
+ }
+ if diag != blas.Unit && diag != blas.NonUnit {
+ panic(badDiag)
+ }
+ checkZMatrix('A', n, n, a, lda)
+ checkZVector('x', n, x, incX)
+
+ if n == 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 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] -= c128.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] -= c128.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] -= c128.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] -= c128.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 {
+ c128.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 {
+ c128.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 {
+ c128.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 {
+ c128.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
+ }
+ }
+ }
+}
OSDN Git Service
