+++ /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
