vendor/gonum.org/v1/gonum/blas/gonum/level2cmplx128.go

   1 // Copyright ©2017 The Gonum Authors. All rights reserved.
   2 // Use of this source code is governed by a BSD-style
   3 // license that can be found in the LICENSE file.
   4
   5 package gonum
   6
   7 import (
   8         "math/cmplx"
   9
  10         "gonum.org/v1/gonum/blas"
  11         "gonum.org/v1/gonum/internal/asm/c128"
  12 )
  13
  14 // Zgemv performs one of the matrix-vector operations
  15 //  y = alpha * A * x + beta * y    if trans = blas.NoTrans
  16 //  y = alpha * A^T * x + beta * y  if trans = blas.Trans
  17 //  y = alpha * A^H * x + beta * y  if trans = blas.ConjTrans
  18 // where alpha and beta are scalars, x and y are vectors, and A is an m×n dense matrix.
  19 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) {
  20         checkZMatrix('A', m, n, a, lda)
  21         switch trans {
  22         default:
  23                 panic(badTranspose)
  24         case blas.NoTrans:
  25                 checkZVector('x', n, x, incX)
  26                 checkZVector('y', m, y, incY)
  27         case blas.Trans, blas.ConjTrans:
  28                 checkZVector('x', m, x, incX)
  29                 checkZVector('y', n, y, incY)
  30         }
  31
  32         if m == 0 || n == 0 || (alpha == 0 && beta == 1) {
  33                 return
  34         }
  35
  36         var lenX, lenY int
  37         if trans == blas.NoTrans {
  38                 lenX = n
  39                 lenY = m
  40         } else {
  41                 lenX = m
  42                 lenY = n
  43         }
  44         var kx int
  45         if incX < 0 {
  46                 kx = (1 - lenX) * incX
  47         }
  48         var ky int
  49         if incY < 0 {
  50                 ky = (1 - lenY) * incY
  51         }
  52
  53         // Form y = beta*y.
  54         if beta != 1 {
  55                 if incY == 1 {
  56                         if beta == 0 {
  57                                 for i := range y[:lenY] {
  58                                         y[i] = 0
  59                                 }
  60                         } else {
  61                                 c128.ScalUnitary(beta, y[:lenY])
  62                         }
  63                 } else {
  64                         iy := ky
  65                         if beta == 0 {
  66                                 for i := 0; i < lenY; i++ {
  67                                         y[iy] = 0
  68                                         iy += incY
  69                                 }
  70                         } else {
  71                                 if incY > 0 {
  72                                         c128.ScalInc(beta, y, uintptr(lenY), uintptr(incY))
  73                                 } else {
  74                                         c128.ScalInc(beta, y, uintptr(lenY), uintptr(-incY))
  75                                 }
  76                         }
  77                 }
  78         }
  79
  80         if alpha == 0 {
  81                 return
  82         }
  83
  84         switch trans {
  85         default:
  86                 // Form y = alpha*A*x + y.
  87                 iy := ky
  88                 if incX == 1 {
  89                         for i := 0; i < m; i++ {
  90                                 y[iy] += alpha * c128.DotuUnitary(a[i*lda:i*lda+n], x[:n])
  91                                 iy += incY
  92                         }
  93                         return
  94                 }
  95                 for i := 0; i < m; i++ {
  96                         y[iy] += alpha * c128.DotuInc(a[i*lda:i*lda+n], x, uintptr(n), 1, uintptr(incX), 0, uintptr(kx))
  97                         iy += incY
  98                 }
  99                 return
 100
 101         case blas.Trans:
 102                 // Form y = alpha*A^T*x + y.
 103                 ix := kx
 104                 if incY == 1 {
 105                         for i := 0; i < m; i++ {
 106                                 c128.AxpyUnitary(alpha*x[ix], a[i*lda:i*lda+n], y[:n])
 107                                 ix += incX
 108                         }
 109                         return
 110                 }
 111                 for i := 0; i < m; i++ {
 112                         c128.AxpyInc(alpha*x[ix], a[i*lda:i*lda+n], y, uintptr(n), 1, uintptr(incY), 0, uintptr(ky))
 113                         ix += incX
 114                 }
 115                 return
 116
 117         case blas.ConjTrans:
 118                 // Form y = alpha*A^H*x + y.
 119                 ix := kx
 120                 if incY == 1 {
 121                         for i := 0; i < m; i++ {
 122                                 tmp := alpha * x[ix]
 123                                 for j := 0; j < n; j++ {
 124                                         y[j] += tmp * cmplx.Conj(a[i*lda+j])
 125                                 }
 126                                 ix += incX
 127                         }
 128                         return
 129                 }
 130                 for i := 0; i < m; i++ {
 131                         tmp := alpha * x[ix]
 132                         jy := ky
 133                         for j := 0; j < n; j++ {
 134                                 y[jy] += tmp * cmplx.Conj(a[i*lda+j])
 135                                 jy += incY
 136                         }
 137                         ix += incX
 138                 }
 139                 return
 140         }
 141 }
 142
 143 // Zgerc performs the rank-one operation
 144 //  A += alpha * x * y^H
 145 // where A is an m×n dense matrix, alpha is a scalar, x is an m element vector,
 146 // and y is an n element vector.
 147 func (Implementation) Zgerc(m, n int, alpha complex128, x []complex128, incX int, y []complex128, incY int, a []complex128, lda int) {
 148         checkZMatrix('A', m, n, a, lda)
 149         checkZVector('x', m, x, incX)
 150         checkZVector('y', n, y, incY)
 151
 152         if m == 0 || n == 0 || alpha == 0 {
 153                 return
 154         }
 155
 156         var kx, jy int
 157         if incX < 0 {
 158                 kx = (1 - m) * incX
 159         }
 160         if incY < 0 {
 161                 jy = (1 - n) * incY
 162         }
 163         for j := 0; j < n; j++ {
 164                 if y[jy] != 0 {
 165                         tmp := alpha * cmplx.Conj(y[jy])
 166                         c128.AxpyInc(tmp, x, a[j:], uintptr(m), uintptr(incX), uintptr(lda), uintptr(kx), 0)
 167                 }
 168                 jy += incY
 169         }
 170 }
 171
 172 // Zgeru performs the rank-one operation
 173 //  A += alpha * x * y^T
 174 // where A is an m×n dense matrix, alpha is a scalar, x is an m element vector,
 175 // and y is an n element vector.
 176 func (Implementation) Zgeru(m, n int, alpha complex128, x []complex128, incX int, y []complex128, incY int, a []complex128, lda int) {
 177         checkZMatrix('A', m, n, a, lda)
 178         checkZVector('x', m, x, incX)
 179         checkZVector('y', n, y, incY)
 180
 181         if m == 0 || n == 0 || alpha == 0 {
 182                 return
 183         }
 184
 185         var kx int
 186         if incX < 0 {
 187                 kx = (1 - m) * incX
 188         }
 189         if incY == 1 {
 190                 for i := 0; i < m; i++ {
 191                         if x[kx] != 0 {
 192                                 tmp := alpha * x[kx]
 193                                 c128.AxpyUnitary(tmp, y[:n], a[i*lda:i*lda+n])
 194                         }
 195                         kx += incX
 196                 }
 197                 return
 198         }
 199         var jy int
 200         if incY < 0 {
 201                 jy = (1 - n) * incY
 202         }
 203         for i := 0; i < m; i++ {
 204                 if x[kx] != 0 {
 205                         tmp := alpha * x[kx]
 206                         c128.AxpyInc(tmp, y, a[i*lda:i*lda+n], uintptr(n), uintptr(incY), 1, uintptr(jy), 0)
 207                 }
 208                 kx += incX
 209         }
 210 }
 211
 212 // Zhemv performs the matrix-vector operation
 213 //  y = alpha * A * x + beta * y
 214 // where alpha and beta are scalars, x and y are vectors, and A is an n×n
 215 // Hermitian matrix. The imaginary parts of the diagonal elements of A are
 216 // ignored and assumed to be zero.
 217 func (Implementation) Zhemv(uplo blas.Uplo, n int, alpha complex128, a []complex128, lda int, x []complex128, incX int, beta complex128, y []complex128, incY int) {
 218         if uplo != blas.Upper && uplo != blas.Lower {
 219                 panic(badUplo)
 220         }
 221         checkZMatrix('A', n, n, a, lda)
 222         checkZVector('x', n, x, incX)
 223         checkZVector('y', n, y, incY)
 224
 225         if n == 0 || (alpha == 0 && beta == 1) {
 226                 return
 227         }
 228
 229         // Set up the start indices in X and Y.
 230         var kx int
 231         if incX < 0 {
 232                 kx = (1 - n) * incX
 233         }
 234         var ky int
 235         if incY < 0 {
 236                 ky = (1 - n) * incY
 237         }
 238
 239         // Form y = beta*y.
 240         if beta != 1 {
 241                 if incY == 1 {
 242                         if beta == 0 {
 243                                 for i := range y[:n] {
 244                                         y[i] = 0
 245                                 }
 246                         } else {
 247                                 for i, v := range y[:n] {
 248                                         y[i] = beta * v
 249                                 }
 250                         }
 251                 } else {
 252                         iy := ky
 253                         if beta == 0 {
 254                                 for i := 0; i < n; i++ {
 255                                         y[iy] = 0
 256                                         iy += incY
 257                                 }
 258                         } else {
 259                                 for i := 0; i < n; i++ {
 260                                         y[iy] = beta * y[iy]
 261                                         iy += incY
 262                                 }
 263                         }
 264                 }
 265         }
 266
 267         if alpha == 0 {
 268                 return
 269         }
 270
 271         // The elements of A are accessed sequentially with one pass through
 272         // the triangular part of A.
 273
 274         if uplo == blas.Upper {
 275                 // Form y when A is stored in upper triangle.
 276                 if incX == 1 && incY == 1 {
 277                         for i := 0; i < n; i++ {
 278                                 tmp1 := alpha * x[i]
 279                                 var tmp2 complex128
 280                                 for j := i + 1; j < n; j++ {
 281                                         y[j] += tmp1 * cmplx.Conj(a[i*lda+j])
 282                                         tmp2 += a[i*lda+j] * x[j]
 283                                 }
 284                                 aii := complex(real(a[i*lda+i]), 0)
 285                                 y[i] += tmp1*aii + alpha*tmp2
 286                         }
 287                 } else {
 288                         ix := kx
 289                         iy := ky
 290                         for i := 0; i < n; i++ {
 291                                 tmp1 := alpha * x[ix]
 292                                 var tmp2 complex128
 293                                 jx := ix
 294                                 jy := iy
 295                                 for j := i + 1; j < n; j++ {
 296                                         jx += incX
 297                                         jy += incY
 298                                         y[jy] += tmp1 * cmplx.Conj(a[i*lda+j])
 299                                         tmp2 += a[i*lda+j] * x[jx]
 300                                 }
 301                                 aii := complex(real(a[i*lda+i]), 0)
 302                                 y[iy] += tmp1*aii + alpha*tmp2
 303                                 ix += incX
 304                                 iy += incY
 305                         }
 306                 }
 307                 return
 308         }
 309
 310         // Form y when A is stored in lower triangle.
 311         if incX == 1 && incY == 1 {
 312                 for i := 0; i < n; i++ {
 313                         tmp1 := alpha * x[i]
 314                         var tmp2 complex128
 315                         for j := 0; j < i; j++ {
 316                                 y[j] += tmp1 * cmplx.Conj(a[i*lda+j])
 317                                 tmp2 += a[i*lda+j] * x[j]
 318                         }
 319                         aii := complex(real(a[i*lda+i]), 0)
 320                         y[i] += tmp1*aii + alpha*tmp2
 321                 }
 322         } else {
 323                 ix := kx
 324                 iy := ky
 325                 for i := 0; i < n; i++ {
 326                         tmp1 := alpha * x[ix]
 327                         var tmp2 complex128
 328                         jx := kx
 329                         jy := ky
 330                         for j := 0; j < i; j++ {
 331                                 y[jy] += tmp1 * cmplx.Conj(a[i*lda+j])
 332                                 tmp2 += a[i*lda+j] * x[jx]
 333                                 jx += incX
 334                                 jy += incY
 335                         }
 336                         aii := complex(real(a[i*lda+i]), 0)
 337                         y[iy] += tmp1*aii + alpha*tmp2
 338                         ix += incX
 339                         iy += incY
 340                 }
 341         }
 342 }
 343
 344 // Zher performs the Hermitian rank-one operation
 345 //  A += alpha * x * x^H
 346 // where A is an n×n Hermitian matrix, alpha is a real scalar, and x is an n
 347 // element vector. On entry, the imaginary parts of the diagonal elements of A
 348 // are ignored and assumed to be zero, on return they will be set to zero.
 349 func (Implementation) Zher(uplo blas.Uplo, n int, alpha float64, x []complex128, incX int, a []complex128, lda int) {
 350         if uplo != blas.Upper && uplo != blas.Lower {
 351                 panic(badUplo)
 352         }
 353         checkZMatrix('A', n, n, a, lda)
 354         checkZVector('x', n, x, incX)
 355
 356         if n == 0 || alpha == 0 {
 357                 return
 358         }
 359
 360         var kx int
 361         if incX < 0 {
 362                 kx = (1 - n) * incX
 363         }
 364         if uplo == blas.Upper {
 365                 if incX == 1 {
 366                         for i := 0; i < n; i++ {
 367                                 if x[i] != 0 {
 368                                         tmp := complex(alpha*real(x[i]), alpha*imag(x[i]))
 369                                         aii := real(a[i*lda+i])
 370                                         xtmp := real(tmp * cmplx.Conj(x[i]))
 371                                         a[i*lda+i] = complex(aii+xtmp, 0)
 372                                         for j := i + 1; j < n; j++ {
 373                                                 a[i*lda+j] += tmp * cmplx.Conj(x[j])
 374                                         }
 375                                 } else {
 376                                         aii := real(a[i*lda+i])
 377                                         a[i*lda+i] = complex(aii, 0)
 378                                 }
 379                         }
 380                         return
 381                 }
 382
 383                 ix := kx
 384                 for i := 0; i < n; i++ {
 385                         if x[ix] != 0 {
 386                                 tmp := complex(alpha*real(x[ix]), alpha*imag(x[ix]))
 387                                 aii := real(a[i*lda+i])
 388                                 xtmp := real(tmp * cmplx.Conj(x[ix]))
 389                                 a[i*lda+i] = complex(aii+xtmp, 0)
 390                                 jx := ix + incX
 391                                 for j := i + 1; j < n; j++ {
 392                                         a[i*lda+j] += tmp * cmplx.Conj(x[jx])
 393                                         jx += incX
 394                                 }
 395                         } else {
 396                                 aii := real(a[i*lda+i])
 397                                 a[i*lda+i] = complex(aii, 0)
 398                         }
 399                         ix += incX
 400                 }
 401                 return
 402         }
 403
 404         if incX == 1 {
 405                 for i := 0; i < n; i++ {
 406                         if x[i] != 0 {
 407                                 tmp := complex(alpha*real(x[i]), alpha*imag(x[i]))
 408                                 for j := 0; j < i; j++ {
 409                                         a[i*lda+j] += tmp * cmplx.Conj(x[j])
 410                                 }
 411                                 aii := real(a[i*lda+i])
 412                                 xtmp := real(tmp * cmplx.Conj(x[i]))
 413                                 a[i*lda+i] = complex(aii+xtmp, 0)
 414                         } else {
 415                                 aii := real(a[i*lda+i])
 416                                 a[i*lda+i] = complex(aii, 0)
 417                         }
 418                 }
 419                 return
 420         }
 421
 422         ix := kx
 423         for i := 0; i < n; i++ {
 424                 if x[ix] != 0 {
 425                         tmp := complex(alpha*real(x[ix]), alpha*imag(x[ix]))
 426                         jx := kx
 427                         for j := 0; j < i; j++ {
 428                                 a[i*lda+j] += tmp * cmplx.Conj(x[jx])
 429                                 jx += incX
 430                         }
 431                         aii := real(a[i*lda+i])
 432                         xtmp := real(tmp * cmplx.Conj(x[ix]))
 433                         a[i*lda+i] = complex(aii+xtmp, 0)
 434
 435                 } else {
 436                         aii := real(a[i*lda+i])
 437                         a[i*lda+i] = complex(aii, 0)
 438                 }
 439                 ix += incX
 440         }
 441 }
 442
 443 // Zher2 performs the Hermitian rank-two operation
 444 //  A += alpha*x*y^H + conj(alpha)*y*x^H
 445 // where alpha is a scalar, x and y are n element vectors and A is an n×n
 446 // Hermitian matrix. On entry, the imaginary parts of the diagonal elements are
 447 // ignored and assumed to be zero. On return they will be set to zero.
 448 func (Implementation) Zher2(uplo blas.Uplo, n int, alpha complex128, x []complex128, incX int, y []complex128, incY int, a []complex128, lda int) {
 449         if uplo != blas.Upper && uplo != blas.Lower {
 450                 panic(badUplo)
 451         }
 452         checkZMatrix('A', n, n, a, lda)
 453         checkZVector('x', n, x, incX)
 454         checkZVector('y', n, y, incY)
 455
 456         if n == 0 || alpha == 0 {
 457                 return
 458         }
 459
 460         var kx, ky int
 461         var ix, iy int
 462         if incX != 1 || incY != 1 {
 463                 if incX < 0 {
 464                         kx = (1 - n) * incX
 465                 }
 466                 if incY < 0 {
 467                         ky = (1 - n) * incY
 468                 }
 469                 ix = kx
 470                 iy = ky
 471         }
 472         if uplo == blas.Upper {
 473                 if incX == 1 && incY == 1 {
 474                         for i := 0; i < n; i++ {
 475                                 if x[i] != 0 || y[i] != 0 {
 476                                         tmp1 := alpha * x[i]
 477                                         tmp2 := cmplx.Conj(alpha) * y[i]
 478                                         aii := real(a[i*lda+i]) + real(tmp1*cmplx.Conj(y[i])) + real(tmp2*cmplx.Conj(x[i]))
 479                                         a[i*lda+i] = complex(aii, 0)
 480                                         for j := i + 1; j < n; j++ {
 481                                                 a[i*lda+j] += tmp1*cmplx.Conj(y[j]) + tmp2*cmplx.Conj(x[j])
 482                                         }
 483                                 } else {
 484                                         aii := real(a[i*lda+i])
 485                                         a[i*lda+i] = complex(aii, 0)
 486                                 }
 487                         }
 488                         return
 489                 }
 490                 for i := 0; i < n; i++ {
 491                         if x[ix] != 0 || y[iy] != 0 {
 492                                 tmp1 := alpha * x[ix]
 493                                 tmp2 := cmplx.Conj(alpha) * y[iy]
 494                                 aii := real(a[i*lda+i]) + real(tmp1*cmplx.Conj(y[iy])) + real(tmp2*cmplx.Conj(x[ix]))
 495                                 a[i*lda+i] = complex(aii, 0)
 496                                 jx := ix + incX
 497                                 jy := iy + incY
 498                                 for j := i + 1; j < n; j++ {
 499                                         a[i*lda+j] += tmp1*cmplx.Conj(y[jy]) + tmp2*cmplx.Conj(x[jx])
 500                                         jx += incX
 501                                         jy += incY
 502                                 }
 503                         } else {
 504                                 aii := real(a[i*lda+i])
 505                                 a[i*lda+i] = complex(aii, 0)
 506                         }
 507                         ix += incX
 508                         iy += incY
 509                 }
 510                 return
 511         }
 512
 513         if incX == 1 && incY == 1 {
 514                 for i := 0; i < n; i++ {
 515                         if x[i] != 0 || y[i] != 0 {
 516                                 tmp1 := alpha * x[i]
 517                                 tmp2 := cmplx.Conj(alpha) * y[i]
 518                                 for j := 0; j < i; j++ {
 519                                         a[i*lda+j] += tmp1*cmplx.Conj(y[j]) + tmp2*cmplx.Conj(x[j])
 520                                 }
 521                                 aii := real(a[i*lda+i]) + real(tmp1*cmplx.Conj(y[i])) + real(tmp2*cmplx.Conj(x[i]))
 522                                 a[i*lda+i] = complex(aii, 0)
 523                         } else {
 524                                 aii := real(a[i*lda+i])
 525                                 a[i*lda+i] = complex(aii, 0)
 526                         }
 527                 }
 528                 return
 529         }
 530         for i := 0; i < n; i++ {
 531                 if x[ix] != 0 || y[iy] != 0 {
 532                         tmp1 := alpha * x[ix]
 533                         tmp2 := cmplx.Conj(alpha) * y[iy]
 534                         jx := kx
 535                         jy := ky
 536                         for j := 0; j < i; j++ {
 537                                 a[i*lda+j] += tmp1*cmplx.Conj(y[jy]) + tmp2*cmplx.Conj(x[jx])
 538                                 jx += incX
 539                                 jy += incY
 540                         }
 541                         aii := real(a[i*lda+i]) + real(tmp1*cmplx.Conj(y[iy])) + real(tmp2*cmplx.Conj(x[ix]))
 542                         a[i*lda+i] = complex(aii, 0)
 543                 } else {
 544                         aii := real(a[i*lda+i])
 545                         a[i*lda+i] = complex(aii, 0)
 546                 }
 547                 ix += incX
 548                 iy += incY
 549         }
 550 }
 551
 552 // Zhpmv performs the matrix-vector operation
 553 //  y = alpha * A * x + beta * y
 554 // where alpha and beta are scalars, x and y are vectors, and A is an n×n
 555 // Hermitian matrix in packed form. The imaginary parts of the diagonal
 556 // elements of A are ignored and assumed to be zero.
 557 func (Implementation) Zhpmv(uplo blas.Uplo, n int, alpha complex128, ap []complex128, x []complex128, incX int, beta complex128, y []complex128, incY int) {
 558         if uplo != blas.Upper && uplo != blas.Lower {
 559                 panic(badUplo)
 560         }
 561         checkZVector('x', n, x, incX)
 562         checkZVector('y', n, y, incY)
 563         if len(ap) < n*(n+1)/2 {
 564                 panic("blas: insufficient A packed matrix slice length")
 565         }
 566
 567         if n == 0 || (alpha == 0 && beta == 1) {
 568                 return
 569         }
 570
 571         // Set up the start indices in X and Y.
 572         var kx int
 573         if incX < 0 {
 574                 kx = (1 - n) * incX
 575         }
 576         var ky int
 577         if incY < 0 {
 578                 ky = (1 - n) * incY
 579         }
 580
 581         // Form y = beta*y.
 582         if beta != 1 {
 583                 if incY == 1 {
 584                         if beta == 0 {
 585                                 for i := range y[:n] {
 586                                         y[i] = 0
 587                                 }
 588                         } else {
 589                                 for i, v := range y[:n] {
 590                                         y[i] = beta * v
 591                                 }
 592                         }
 593                 } else {
 594                         iy := ky
 595                         if beta == 0 {
 596                                 for i := 0; i < n; i++ {
 597                                         y[iy] = 0
 598                                         iy += incY
 599                                 }
 600                         } else {
 601                                 for i := 0; i < n; i++ {
 602                                         y[iy] *= beta
 603                                         iy += incY
 604                                 }
 605                         }
 606                 }
 607         }
 608
 609         if alpha == 0 {
 610                 return
 611         }
 612
 613         // The elements of A are accessed sequentially with one pass through ap.
 614
 615         var kk int
 616         if uplo == blas.Upper {
 617                 // Form y when ap contains the upper triangle.
 618                 // Here, kk points to the current diagonal element in ap.
 619                 if incX == 1 && incY == 1 {
 620                         for i := 0; i < n; i++ {
 621                                 tmp1 := alpha * x[i]
 622                                 y[i] += tmp1 * complex(real(ap[kk]), 0)
 623                                 var tmp2 complex128
 624                                 k := kk + 1
 625                                 for j := i + 1; j < n; j++ {
 626                                         y[j] += tmp1 * cmplx.Conj(ap[k])
 627                                         tmp2 += ap[k] * x[j]
 628                                         k++
 629                                 }
 630                                 y[i] += alpha * tmp2
 631                                 kk += n - i
 632                         }
 633                 } else {
 634                         ix := kx
 635                         iy := ky
 636                         for i := 0; i < n; i++ {
 637                                 tmp1 := alpha * x[ix]
 638                                 y[iy] += tmp1 * complex(real(ap[kk]), 0)
 639                                 var tmp2 complex128
 640                                 jx := ix
 641                                 jy := iy
 642                                 for k := kk + 1; k < kk+n-i; k++ {
 643                                         jx += incX
 644                                         jy += incY
 645                                         y[jy] += tmp1 * cmplx.Conj(ap[k])
 646                                         tmp2 += ap[k] * x[jx]
 647                                 }
 648                                 y[iy] += alpha * tmp2
 649                                 ix += incX
 650                                 iy += incY
 651                                 kk += n - i
 652                         }
 653                 }
 654                 return
 655         }
 656
 657         // Form y when ap contains the lower triangle.
 658         // Here, kk points to the beginning of current row in ap.
 659         if incX == 1 && incY == 1 {
 660                 for i := 0; i < n; i++ {
 661                         tmp1 := alpha * x[i]
 662                         var tmp2 complex128
 663                         k := kk
 664                         for j := 0; j < i; j++ {
 665                                 y[j] += tmp1 * cmplx.Conj(ap[k])
 666                                 tmp2 += ap[k] * x[j]
 667                                 k++
 668                         }
 669                         aii := complex(real(ap[kk+i]), 0)
 670                         y[i] += tmp1*aii + alpha*tmp2
 671                         kk += i + 1
 672                 }
 673         } else {
 674                 ix := kx
 675                 iy := ky
 676                 for i := 0; i < n; i++ {
 677                         tmp1 := alpha * x[ix]
 678                         var tmp2 complex128
 679                         jx := kx
 680                         jy := ky
 681                         for k := kk; k < kk+i; k++ {
 682                                 y[jy] += tmp1 * cmplx.Conj(ap[k])
 683                                 tmp2 += ap[k] * x[jx]
 684                                 jx += incX
 685                                 jy += incY
 686                         }
 687                         aii := complex(real(ap[kk+i]), 0)
 688                         y[iy] += tmp1*aii + alpha*tmp2
 689                         ix += incX
 690                         iy += incY
 691                         kk += i + 1
 692                 }
 693         }
 694 }
 695
 696 // Zhpr performs the Hermitian rank-1 operation
 697 //  A += alpha * x * x^H,
 698 // where alpha is a real scalar, x is a vector, and A is an n×n hermitian matrix
 699 // in packed form. On entry, the imaginary parts of the diagonal elements are
 700 // assumed to be zero, and on return they are set to zero.
 701 func (Implementation) Zhpr(uplo blas.Uplo, n int, alpha float64, x []complex128, incX int, ap []complex128) {
 702         if uplo != blas.Upper && uplo != blas.Lower {
 703                 panic(badUplo)
 704         }
 705         if n < 0 {
 706                 panic(nLT0)
 707         }
 708         checkZVector('x', n, x, incX)
 709         if len(ap) < n*(n+1)/2 {
 710                 panic("blas: insufficient A packed matrix slice length")
 711         }
 712
 713         if n == 0 || alpha == 0 {
 714                 return
 715         }
 716
 717         // Set up start index in X.
 718         var kx int
 719         if incX < 0 {
 720                 kx = (1 - n) * incX
 721         }
 722
 723         // The elements of A are accessed sequentially with one pass through ap.
 724
 725         var kk int
 726         if uplo == blas.Upper {
 727                 // Form A when upper triangle is stored in AP.
 728                 // Here, kk points to the current diagonal element in ap.
 729                 if incX == 1 {
 730                         for i := 0; i < n; i++ {
 731                                 xi := x[i]
 732                                 if xi != 0 {
 733                                         aii := real(ap[kk]) + alpha*real(cmplx.Conj(xi)*xi)
 734                                         ap[kk] = complex(aii, 0)
 735
 736                                         tmp := complex(alpha, 0) * xi
 737                                         a := ap[kk+1 : kk+n-i]
 738                                         x := x[i+1 : n]
 739                                         for j, v := range x {
 740                                                 a[j] += tmp * cmplx.Conj(v)
 741                                         }
 742                                 } else {
 743                                         ap[kk] = complex(real(ap[kk]), 0)
 744                                 }
 745                                 kk += n - i
 746                         }
 747                 } else {
 748                         ix := kx
 749                         for i := 0; i < n; i++ {
 750                                 xi := x[ix]
 751                                 if xi != 0 {
 752                                         aii := real(ap[kk]) + alpha*real(cmplx.Conj(xi)*xi)
 753                                         ap[kk] = complex(aii, 0)
 754
 755                                         tmp := complex(alpha, 0) * xi
 756                                         jx := ix + incX
 757                                         a := ap[kk+1 : kk+n-i]
 758                                         for k := range a {
 759                                                 a[k] += tmp * cmplx.Conj(x[jx])
 760                                                 jx += incX
 761                                         }
 762                                 } else {
 763                                         ap[kk] = complex(real(ap[kk]), 0)
 764                                 }
 765                                 ix += incX
 766                                 kk += n - i
 767                         }
 768                 }
 769                 return
 770         }
 771
 772         // Form A when lower triangle is stored in AP.
 773         // Here, kk points to the beginning of current row in ap.
 774         if incX == 1 {
 775                 for i := 0; i < n; i++ {
 776                         xi := x[i]
 777                         if xi != 0 {
 778                                 tmp := complex(alpha, 0) * xi
 779                                 a := ap[kk : kk+i]
 780                                 for j, v := range x[:i] {
 781                                         a[j] += tmp * cmplx.Conj(v)
 782                                 }
 783
 784                                 aii := real(ap[kk+i]) + alpha*real(cmplx.Conj(xi)*xi)
 785                                 ap[kk+i] = complex(aii, 0)
 786                         } else {
 787                                 ap[kk+i] = complex(real(ap[kk+i]), 0)
 788                         }
 789                         kk += i + 1
 790                 }
 791         } else {
 792                 ix := kx
 793                 for i := 0; i < n; i++ {
 794                         xi := x[ix]
 795                         if xi != 0 {
 796                                 tmp := complex(alpha, 0) * xi
 797                                 a := ap[kk : kk+i]
 798                                 jx := kx
 799                                 for k := range a {
 800                                         a[k] += tmp * cmplx.Conj(x[jx])
 801                                         jx += incX
 802                                 }
 803
 804                                 aii := real(ap[kk+i]) + alpha*real(cmplx.Conj(xi)*xi)
 805                                 ap[kk+i] = complex(aii, 0)
 806                         } else {
 807                                 ap[kk+i] = complex(real(ap[kk+i]), 0)
 808                         }
 809                         ix += incX
 810                         kk += i + 1
 811                 }
 812         }
 813 }
 814
 815 // Zhpr2 performs the Hermitian rank-2 operation
 816 //  A += alpha*x*y^H + conj(alpha)*y*x^H,
 817 // where alpha is a complex scalar, x and y are n element vectors, and A is an
 818 // n×n Hermitian matrix, supplied in packed form. On entry, the imaginary parts
 819 // of the diagonal elements are assumed to be zero, and on return they are set to zero.
 820 func (Implementation) Zhpr2(uplo blas.Uplo, n int, alpha complex128, x []complex128, incX int, y []complex128, incY int, ap []complex128) {
 821         if uplo != blas.Upper && uplo != blas.Lower {
 822                 panic(badUplo)
 823         }
 824         if n < 0 {
 825                 panic(nLT0)
 826         }
 827         checkZVector('x', n, x, incX)
 828         checkZVector('y', n, y, incY)
 829         if len(ap) < n*(n+1)/2 {
 830                 panic("blas: insufficient A packed matrix slice length")
 831         }
 832
 833         if n == 0 || alpha == 0 {
 834                 return
 835         }
 836
 837         // Set up start indices in X and Y.
 838         var kx int
 839         if incX < 0 {
 840                 kx = (1 - n) * incX
 841         }
 842         var ky int
 843         if incY < 0 {
 844                 ky = (1 - n) * incY
 845         }
 846
 847         // The elements of A are accessed sequentially with one pass through ap.
 848
 849         var kk int
 850         if uplo == blas.Upper {
 851                 // Form A when upper triangle is stored in AP.
 852                 // Here, kk points to the current diagonal element in ap.
 853                 if incX == 1 && incY == 1 {
 854                         for i := 0; i < n; i++ {
 855                                 if x[i] != 0 || y[i] != 0 {
 856                                         tmp1 := alpha * x[i]
 857                                         tmp2 := cmplx.Conj(alpha) * y[i]
 858                                         aii := real(ap[kk]) + real(tmp1*cmplx.Conj(y[i])) + real(tmp2*cmplx.Conj(x[i]))
 859                                         ap[kk] = complex(aii, 0)
 860                                         k := kk + 1
 861                                         for j := i + 1; j < n; j++ {
 862                                                 ap[k] += tmp1*cmplx.Conj(y[j]) + tmp2*cmplx.Conj(x[j])
 863                                                 k++
 864                                         }
 865                                 } else {
 866                                         ap[kk] = complex(real(ap[kk]), 0)
 867                                 }
 868                                 kk += n - i
 869                         }
 870                 } else {
 871                         ix := kx
 872                         iy := ky
 873                         for i := 0; i < n; i++ {
 874                                 if x[ix] != 0 || y[iy] != 0 {
 875                                         tmp1 := alpha * x[ix]
 876                                         tmp2 := cmplx.Conj(alpha) * y[iy]
 877                                         aii := real(ap[kk]) + real(tmp1*cmplx.Conj(y[iy])) + real(tmp2*cmplx.Conj(x[ix]))
 878                                         ap[kk] = complex(aii, 0)
 879                                         jx := ix + incX
 880                                         jy := iy + incY
 881                                         for k := kk + 1; k < kk+n-i; k++ {
 882                                                 ap[k] += tmp1*cmplx.Conj(y[jy]) + tmp2*cmplx.Conj(x[jx])
 883                                                 jx += incX
 884                                                 jy += incY
 885                                         }
 886                                 } else {
 887                                         ap[kk] = complex(real(ap[kk]), 0)
 888                                 }
 889                                 ix += incX
 890                                 iy += incY
 891                                 kk += n - i
 892                         }
 893                 }
 894                 return
 895         }
 896
 897         // Form A when lower triangle is stored in AP.
 898         // Here, kk points to the beginning of current row in ap.
 899         if incX == 1 && incY == 1 {
 900                 for i := 0; i < n; i++ {
 901                         if x[i] != 0 || y[i] != 0 {
 902                                 tmp1 := alpha * x[i]
 903                                 tmp2 := cmplx.Conj(alpha) * y[i]
 904                                 k := kk
 905                                 for j := 0; j < i; j++ {
 906                                         ap[k] += tmp1*cmplx.Conj(y[j]) + tmp2*cmplx.Conj(x[j])
 907                                         k++
 908                                 }
 909                                 aii := real(ap[kk+i]) + real(tmp1*cmplx.Conj(y[i])) + real(tmp2*cmplx.Conj(x[i]))
 910                                 ap[kk+i] = complex(aii, 0)
 911                         } else {
 912                                 ap[kk+i] = complex(real(ap[kk+i]), 0)
 913                         }
 914                         kk += i + 1
 915                 }
 916         } else {
 917                 ix := kx
 918                 iy := ky
 919                 for i := 0; i < n; i++ {
 920                         if x[ix] != 0 || y[iy] != 0 {
 921                                 tmp1 := alpha * x[ix]
 922                                 tmp2 := cmplx.Conj(alpha) * y[iy]
 923                                 jx := kx
 924                                 jy := ky
 925                                 for k := kk; k < kk+i; k++ {
 926                                         ap[k] += tmp1*cmplx.Conj(y[jy]) + tmp2*cmplx.Conj(x[jx])
 927                                         jx += incX
 928                                         jy += incY
 929                                 }
 930                                 aii := real(ap[kk+i]) + real(tmp1*cmplx.Conj(y[iy])) + real(tmp2*cmplx.Conj(x[ix]))
 931                                 ap[kk+i] = complex(aii, 0)
 932                         } else {
 933                                 ap[kk+i] = complex(real(ap[kk+i]), 0)
 934                         }
 935                         ix += incX
 936                         iy += incY
 937                         kk += i + 1
 938                 }
 939         }
 940 }
 941
 942 // Ztpmv performs one of the matrix-vector operations
 943 //  x = A * x    if trans = blas.NoTrans
 944 //  x = A^T * x  if trans = blas.Trans
 945 //  x = A^H * x  if trans = blas.ConjTrans
 946 // where x is an n element vector and A is an n×n triangular matrix, supplied in
 947 // packed form.
 948 func (Implementation) Ztpmv(uplo blas.Uplo, trans blas.Transpose, diag blas.Diag, n int, ap []complex128, x []complex128, incX int) {
 949         if uplo != blas.Upper && uplo != blas.Lower {
 950                 panic(badUplo)
 951         }
 952         if trans != blas.NoTrans && trans != blas.Trans && trans != blas.ConjTrans {
 953                 panic(badTranspose)
 954         }
 955         if diag != blas.Unit && diag != blas.NonUnit {
 956                 panic(badDiag)
 957         }
 958         checkZVector('x', n, x, incX)
 959         if len(ap) < n*(n+1)/2 {
 960                 panic("blas: insufficient A packed matrix slice length")
 961         }
 962
 963         if n == 0 {
 964                 return
 965         }
 966
 967         // Set up start index in X.
 968         var kx int
 969         if incX < 0 {
 970                 kx = (1 - n) * incX
 971         }
 972
 973         // The elements of A are accessed sequentially with one pass through A.
 974
 975         if trans == blas.NoTrans {
 976                 // Form x = A*x.
 977                 if uplo == blas.Upper {
 978                         // kk points to the current diagonal element in ap.
 979                         kk := 0
 980                         if incX == 1 {
 981                                 x = x[:n]
 982                                 for i := range x {
 983                                         if diag == blas.NonUnit {
 984                                                 x[i] *= ap[kk]
 985                                         }
 986                                         if n-i-1 > 0 {
 987                                                 x[i] += c128.DotuUnitary(ap[kk+1:kk+n-i], x[i+1:])
 988                                         }
 989                                         kk += n - i
 990                                 }
 991                         } else {
 992                                 ix := kx
 993                                 for i := 0; i < n; i++ {
 994                                         if diag == blas.NonUnit {
 995                                                 x[ix] *= ap[kk]
 996                                         }
 997                                         if n-i-1 > 0 {
 998                                                 x[ix] += c128.DotuInc(ap[kk+1:kk+n-i], x, uintptr(n-i-1), 1, uintptr(incX), 0, uintptr(ix+incX))
 999                                         }
1000                                         ix += incX
1001                                         kk += n - i
1002                                 }
1003                         }
1004                 } else {
1005                         // kk points to the beginning of current row in ap.
1006                         kk := n*(n+1)/2 - n
1007                         if incX == 1 {
1008                                 for i := n - 1; i >= 0; i-- {
1009                                         if diag == blas.NonUnit {
1010                                                 x[i] *= ap[kk+i]
1011                                         }
1012                                         if i > 0 {
1013                                                 x[i] += c128.DotuUnitary(ap[kk:kk+i], x[:i])
1014                                         }
1015                                         kk -= i
1016                                 }
1017                         } else {
1018                                 ix := kx + (n-1)*incX
1019                                 for i := n - 1; i >= 0; i-- {
1020                                         if diag == blas.NonUnit {
1021                                                 x[ix] *= ap[kk+i]
1022                                         }
1023                                         if i > 0 {
1024                                                 x[ix] += c128.DotuInc(ap[kk:kk+i], x, uintptr(i), 1, uintptr(incX), 0, uintptr(kx))
1025                                         }
1026                                         ix -= incX
1027                                         kk -= i
1028                                 }
1029                         }
1030                 }
1031                 return
1032         }
1033
1034         if trans == blas.Trans {
1035                 // Form x = A^T*x.
1036                 if uplo == blas.Upper {
1037                         // kk points to the current diagonal element in ap.
1038                         kk := n*(n+1)/2 - 1
1039                         if incX == 1 {
1040                                 for i := n - 1; i >= 0; i-- {
1041                                         xi := x[i]
1042                                         if diag == blas.NonUnit {
1043                                                 x[i] *= ap[kk]
1044                                         }
1045                                         if n-i-1 > 0 {
1046                                                 c128.AxpyUnitary(xi, ap[kk+1:kk+n-i], x[i+1:n])
1047                                         }
1048                                         kk -= n - i + 1
1049                                 }
1050                         } else {
1051                                 ix := kx + (n-1)*incX
1052                                 for i := n - 1; i >= 0; i-- {
1053                                         xi := x[ix]
1054                                         if diag == blas.NonUnit {
1055                                                 x[ix] *= ap[kk]
1056                                         }
1057                                         if n-i-1 > 0 {
1058                                                 c128.AxpyInc(xi, ap[kk+1:kk+n-i], x, uintptr(n-i-1), 1, uintptr(incX), 0, uintptr(ix+incX))
1059                                         }
1060                                         ix -= incX
1061                                         kk -= n - i + 1
1062                                 }
1063                         }
1064                 } else {
1065                         // kk points to the beginning of current row in ap.
1066                         kk := 0
1067                         if incX == 1 {
1068                                 x = x[:n]
1069                                 for i := range x {
1070                                         if i > 0 {
1071                                                 c128.AxpyUnitary(x[i], ap[kk:kk+i], x[:i])
1072                                         }
1073                                         if diag == blas.NonUnit {
1074                                                 x[i] *= ap[kk+i]
1075                                         }
1076                                         kk += i + 1
1077                                 }
1078                         } else {
1079                                 ix := kx
1080                                 for i := 0; i < n; i++ {
1081                                         if i > 0 {
1082                                                 c128.AxpyInc(x[ix], ap[kk:kk+i], x, uintptr(i), 1, uintptr(incX), 0, uintptr(kx))
1083                                         }
1084                                         if diag == blas.NonUnit {
1085                                                 x[ix] *= ap[kk+i]
1086                                         }
1087                                         ix += incX
1088                                         kk += i + 1
1089                                 }
1090                         }
1091                 }
1092                 return
1093         }
1094
1095         // Form x = A^H*x.
1096         if uplo == blas.Upper {
1097                 // kk points to the current diagonal element in ap.
1098                 kk := n*(n+1)/2 - 1
1099                 if incX == 1 {
1100                         for i := n - 1; i >= 0; i-- {
1101                                 xi := x[i]
1102                                 if diag == blas.NonUnit {
1103                                         x[i] *= cmplx.Conj(ap[kk])
1104                                 }
1105                                 k := kk + 1
1106                                 for j := i + 1; j < n; j++ {
1107                                         x[j] += xi * cmplx.Conj(ap[k])
1108                                         k++
1109                                 }
1110                                 kk -= n - i + 1
1111                         }
1112                 } else {
1113                         ix := kx + (n-1)*incX
1114                         for i := n - 1; i >= 0; i-- {
1115                                 xi := x[ix]
1116                                 if diag == blas.NonUnit {
1117                                         x[ix] *= cmplx.Conj(ap[kk])
1118                                 }
1119                                 jx := ix + incX
1120                                 k := kk + 1
1121                                 for j := i + 1; j < n; j++ {
1122                                         x[jx] += xi * cmplx.Conj(ap[k])
1123                                         jx += incX
1124                                         k++
1125                                 }
1126                                 ix -= incX
1127                                 kk -= n - i + 1
1128                         }
1129                 }
1130         } else {
1131                 // kk points to the beginning of current row in ap.
1132                 kk := 0
1133                 if incX == 1 {
1134                         x = x[:n]
1135                         for i, xi := range x {
1136                                 for j := 0; j < i; j++ {
1137                                         x[j] += xi * cmplx.Conj(ap[kk+j])
1138                                 }
1139                                 if diag == blas.NonUnit {
1140                                         x[i] *= cmplx.Conj(ap[kk+i])
1141                                 }
1142                                 kk += i + 1
1143                         }
1144                 } else {
1145                         ix := kx
1146                         for i := 0; i < n; i++ {
1147                                 xi := x[ix]
1148                                 jx := kx
1149                                 for j := 0; j < i; j++ {
1150                                         x[jx] += xi * cmplx.Conj(ap[kk+j])
1151                                         jx += incX
1152                                 }
1153                                 if diag == blas.NonUnit {
1154                                         x[ix] *= cmplx.Conj(ap[kk+i])
1155                                 }
1156                                 ix += incX
1157                                 kk += i + 1
1158                         }
1159                 }
1160         }
1161 }
1162
1163 // Ztrmv performs one of the matrix-vector operations
1164 //  x = A * x    if trans = blas.NoTrans
1165 //  x = A^T * x  if trans = blas.Trans
1166 //  x = A^H * x  if trans = blas.ConjTrans
1167 // where x is a vector, and A is an n×n triangular matrix.
1168 func (Implementation) Ztrmv(uplo blas.Uplo, trans blas.Transpose, diag blas.Diag, n int, a []complex128, lda int, x []complex128, incX int) {
1169         if uplo != blas.Upper && uplo != blas.Lower {
1170                 panic(badUplo)
1171         }
1172         if trans != blas.NoTrans && trans != blas.Trans && trans != blas.ConjTrans {
1173                 panic(badTranspose)
1174         }
1175         if diag != blas.Unit && diag != blas.NonUnit {
1176                 panic(badDiag)
1177         }
1178         checkZMatrix('A', n, n, a, lda)
1179         checkZVector('x', n, x, incX)
1180
1181         if n == 0 {
1182                 return
1183         }
1184
1185         // Set up start index in X.
1186         var kx int
1187         if incX < 0 {
1188                 kx = (1 - n) * incX
1189         }
1190
1191         // The elements of A are accessed sequentially with one pass through A.
1192
1193         if trans == blas.NoTrans {
1194                 // Form x = A*x.
1195                 if uplo == blas.Upper {
1196                         if incX == 1 {
1197                                 for i := 0; i < n; i++ {
1198                                         if diag == blas.NonUnit {
1199                                                 x[i] *= a[i*lda+i]
1200                                         }
1201                                         if n-i-1 > 0 {
1202                                                 x[i] += c128.DotuUnitary(a[i*lda+i+1:i*lda+n], x[i+1:n])
1203                                         }
1204                                 }
1205                         } else {
1206                                 ix := kx
1207                                 for i := 0; i < n; i++ {
1208                                         if diag == blas.NonUnit {
1209                                                 x[ix] *= a[i*lda+i]
1210                                         }
1211                                         if n-i-1 > 0 {
1212                                                 x[ix] += c128.DotuInc(a[i*lda+i+1:i*lda+n], x, uintptr(n-i-1), 1, uintptr(incX), 0, uintptr(ix+incX))
1213                                         }
1214                                         ix += incX
1215                                 }
1216                         }
1217                 } else {
1218                         if incX == 1 {
1219                                 for i := n - 1; i >= 0; i-- {
1220                                         if diag == blas.NonUnit {
1221                                                 x[i] *= a[i*lda+i]
1222                                         }
1223                                         if i > 0 {
1224                                                 x[i] += c128.DotuUnitary(a[i*lda:i*lda+i], x[:i])
1225                                         }
1226                                 }
1227                         } else {
1228                                 ix := kx + (n-1)*incX
1229                                 for i := n - 1; i >= 0; i-- {
1230                                         if diag == blas.NonUnit {
1231                                                 x[ix] *= a[i*lda+i]
1232                                         }
1233                                         if i > 0 {
1234                                                 x[ix] += c128.DotuInc(a[i*lda:i*lda+i], x, uintptr(i), 1, uintptr(incX), 0, uintptr(kx))
1235                                         }
1236                                         ix -= incX
1237                                 }
1238                         }
1239                 }
1240                 return
1241         }
1242
1243         if trans == blas.Trans {
1244                 // Form x = A^T*x.
1245                 if uplo == blas.Upper {
1246                         if incX == 1 {
1247                                 for i := n - 1; i >= 0; i-- {
1248                                         xi := x[i]
1249                                         if diag == blas.NonUnit {
1250                                                 x[i] *= a[i*lda+i]
1251                                         }
1252                                         if n-i-1 > 0 {
1253                                                 c128.AxpyUnitary(xi, a[i*lda+i+1:i*lda+n], x[i+1:n])
1254                                         }
1255                                 }
1256                         } else {
1257                                 ix := kx + (n-1)*incX
1258                                 for i := n - 1; i >= 0; i-- {
1259                                         xi := x[ix]
1260                                         if diag == blas.NonUnit {
1261                                                 x[ix] *= a[i*lda+i]
1262                                         }
1263                                         if n-i-1 > 0 {
1264                                                 c128.AxpyInc(xi, a[i*lda+i+1:i*lda+n], x, uintptr(n-i-1), 1, uintptr(incX), 0, uintptr(ix+incX))
1265                                         }
1266                                         ix -= incX
1267                                 }
1268                         }
1269                 } else {
1270                         if incX == 1 {
1271                                 for i := 0; i < n; i++ {
1272                                         if i > 0 {
1273                                                 c128.AxpyUnitary(x[i], a[i*lda:i*lda+i], x[:i])
1274                                         }
1275                                         if diag == blas.NonUnit {
1276                                                 x[i] *= a[i*lda+i]
1277                                         }
1278                                 }
1279                         } else {
1280                                 ix := kx
1281                                 for i := 0; i < n; i++ {
1282                                         if i > 0 {
1283                                                 c128.AxpyInc(x[ix], a[i*lda:i*lda+i], x, uintptr(i), 1, uintptr(incX), 0, uintptr(kx))
1284                                         }
1285                                         if diag == blas.NonUnit {
1286                                                 x[ix] *= a[i*lda+i]
1287                                         }
1288                                         ix += incX
1289                                 }
1290                         }
1291                 }
1292                 return
1293         }
1294
1295         // Form x = A^H*x.
1296         if uplo == blas.Upper {
1297                 if incX == 1 {
1298                         for i := n - 1; i >= 0; i-- {
1299                                 xi := x[i]
1300                                 if diag == blas.NonUnit {
1301                                         x[i] *= cmplx.Conj(a[i*lda+i])
1302                                 }
1303                                 for j := i + 1; j < n; j++ {
1304                                         x[j] += xi * cmplx.Conj(a[i*lda+j])
1305                                 }
1306                         }
1307                 } else {
1308                         ix := kx + (n-1)*incX
1309                         for i := n - 1; i >= 0; i-- {
1310                                 xi := x[ix]
1311                                 if diag == blas.NonUnit {
1312                                         x[ix] *= cmplx.Conj(a[i*lda+i])
1313                                 }
1314                                 jx := ix + incX
1315                                 for j := i + 1; j < n; j++ {
1316                                         x[jx] += xi * cmplx.Conj(a[i*lda+j])
1317                                         jx += incX
1318                                 }
1319                                 ix -= incX
1320                         }
1321                 }
1322         } else {
1323                 if incX == 1 {
1324                         for i := 0; i < n; i++ {
1325                                 for j := 0; j < i; j++ {
1326                                         x[j] += x[i] * cmplx.Conj(a[i*lda+j])
1327                                 }
1328                                 if diag == blas.NonUnit {
1329                                         x[i] *= cmplx.Conj(a[i*lda+i])
1330                                 }
1331                         }
1332                 } else {
1333                         ix := kx
1334                         for i := 0; i < n; i++ {
1335                                 jx := kx
1336                                 for j := 0; j < i; j++ {
1337                                         x[jx] += x[ix] * cmplx.Conj(a[i*lda+j])
1338                                         jx += incX
1339                                 }
1340                                 if diag == blas.NonUnit {
1341                                         x[ix] *= cmplx.Conj(a[i*lda+i])
1342                                 }
1343                                 ix += incX
1344                         }
1345                 }
1346         }
1347 }
1348
1349 // Ztrsv solves one of the systems of equations
1350 //  A*x = b     if trans == blas.NoTrans,
1351 //  A^T*x = b,  if trans == blas.Trans,
1352 //  A^H*x = b,  if trans == blas.ConjTrans,
1353 // where b and x are n element vectors and A is an n×n triangular matrix.
1354 //
1355 // On entry, x contains the values of b, and the solution is
1356 // stored in-place into x.
1357 //
1358 // No test for singularity or near-singularity is included in this
1359 // routine. Such tests must be performed before calling this routine.
1360 func (Implementation) Ztrsv(uplo blas.Uplo, trans blas.Transpose, diag blas.Diag, n int, a []complex128, lda int, x []complex128, incX int) {
1361         if uplo != blas.Upper && uplo != blas.Lower {
1362                 panic(badUplo)
1363         }
1364         if trans != blas.NoTrans && trans != blas.Trans && trans != blas.ConjTrans {
1365                 panic(badTranspose)
1366         }
1367         if diag != blas.Unit && diag != blas.NonUnit {
1368                 panic(badDiag)
1369         }
1370         checkZMatrix('A', n, n, a, lda)
1371         checkZVector('x', n, x, incX)
1372
1373         if n == 0 {
1374                 return
1375         }
1376
1377         // Set up start index in X.
1378         var kx int
1379         if incX < 0 {
1380                 kx = (1 - n) * incX
1381         }
1382
1383         // The elements of A are accessed sequentially with one pass through A.
1384
1385         if trans == blas.NoTrans {
1386                 // Form x = inv(A)*x.
1387                 if uplo == blas.Upper {
1388                         if incX == 1 {
1389                                 for i := n - 1; i >= 0; i-- {
1390                                         aii := a[i*lda+i]
1391                                         if n-i-1 > 0 {
1392                                                 x[i] -= c128.DotuUnitary(x[i+1:n], a[i*lda+i+1:i*lda+n])
1393                                         }
1394                                         if diag == blas.NonUnit {
1395                                                 x[i] /= aii
1396                                         }
1397                                 }
1398                         } else {
1399                                 ix := kx + (n-1)*incX
1400                                 for i := n - 1; i >= 0; i-- {
1401                                         aii := a[i*lda+i]
1402                                         if n-i-1 > 0 {
1403                                                 x[ix] -= c128.DotuInc(x, a[i*lda+i+1:i*lda+n], uintptr(n-i-1), uintptr(incX), 1, uintptr(ix+incX), 0)
1404                                         }
1405                                         if diag == blas.NonUnit {
1406                                                 x[ix] /= aii
1407                                         }
1408                                         ix -= incX
1409                                 }
1410                         }
1411                 } else {
1412                         if incX == 1 {
1413                                 for i := 0; i < n; i++ {
1414                                         if i > 0 {
1415                                                 x[i] -= c128.DotuUnitary(x[:i], a[i*lda:i*lda+i])
1416                                         }
1417                                         if diag == blas.NonUnit {
1418                                                 x[i] /= a[i*lda+i]
1419                                         }
1420                                 }
1421                         } else {
1422                                 ix := kx
1423                                 for i := 0; i < n; i++ {
1424                                         if i > 0 {
1425                                                 x[ix] -= c128.DotuInc(x, a[i*lda:i*lda+i], uintptr(i), uintptr(incX), 1, uintptr(kx), 0)
1426                                         }
1427                                         if diag == blas.NonUnit {
1428                                                 x[ix] /= a[i*lda+i]
1429                                         }
1430                                         ix += incX
1431                                 }
1432                         }
1433                 }
1434                 return
1435         }
1436
1437         if trans == blas.Trans {
1438                 // Form x = inv(A^T)*x.
1439                 if uplo == blas.Upper {
1440                         if incX == 1 {
1441                                 for j := 0; j < n; j++ {
1442                                         if diag == blas.NonUnit {
1443                                                 x[j] /= a[j*lda+j]
1444                                         }
1445                                         if n-j-1 > 0 {
1446                                                 c128.AxpyUnitary(-x[j], a[j*lda+j+1:j*lda+n], x[j+1:n])
1447                                         }
1448                                 }
1449                         } else {
1450                                 jx := kx
1451                                 for j := 0; j < n; j++ {
1452                                         if diag == blas.NonUnit {
1453                                                 x[jx] /= a[j*lda+j]
1454                                         }
1455                                         if n-j-1 > 0 {
1456                                                 c128.AxpyInc(-x[jx], a[j*lda+j+1:j*lda+n], x, uintptr(n-j-1), 1, uintptr(incX), 0, uintptr(jx+incX))
1457                                         }
1458                                         jx += incX
1459                                 }
1460                         }
1461                 } else {
1462                         if incX == 1 {
1463                                 for j := n - 1; j >= 0; j-- {
1464                                         if diag == blas.NonUnit {
1465                                                 x[j] /= a[j*lda+j]
1466                                         }
1467                                         xj := x[j]
1468                                         if j > 0 {
1469                                                 c128.AxpyUnitary(-xj, a[j*lda:j*lda+j], x[:j])
1470                                         }
1471                                 }
1472                         } else {
1473                                 jx := kx + (n-1)*incX
1474                                 for j := n - 1; j >= 0; j-- {
1475                                         if diag == blas.NonUnit {
1476                                                 x[jx] /= a[j*lda+j]
1477                                         }
1478                                         if j > 0 {
1479                                                 c128.AxpyInc(-x[jx], a[j*lda:j*lda+j], x, uintptr(j), 1, uintptr(incX), 0, uintptr(kx))
1480                                         }
1481                                         jx -= incX
1482                                 }
1483                         }
1484                 }
1485                 return
1486         }
1487
1488         // Form x = inv(A^H)*x.
1489         if uplo == blas.Upper {
1490                 if incX == 1 {
1491                         for j := 0; j < n; j++ {
1492                                 if diag == blas.NonUnit {
1493                                         x[j] /= cmplx.Conj(a[j*lda+j])
1494                                 }
1495                                 xj := x[j]
1496                                 for i := j + 1; i < n; i++ {
1497                                         x[i] -= xj * cmplx.Conj(a[j*lda+i])
1498                                 }
1499                         }
1500                 } else {
1501                         jx := kx
1502                         for j := 0; j < n; j++ {
1503                                 if diag == blas.NonUnit {
1504                                         x[jx] /= cmplx.Conj(a[j*lda+j])
1505                                 }
1506                                 xj := x[jx]
1507                                 ix := jx + incX
1508                                 for i := j + 1; i < n; i++ {
1509                                         x[ix] -= xj * cmplx.Conj(a[j*lda+i])
1510                                         ix += incX
1511                                 }
1512                                 jx += incX
1513                         }
1514                 }
1515         } else {
1516                 if incX == 1 {
1517                         for j := n - 1; j >= 0; j-- {
1518                                 if diag == blas.NonUnit {
1519                                         x[j] /= cmplx.Conj(a[j*lda+j])
1520                                 }
1521                                 xj := x[j]
1522                                 for i := 0; i < j; i++ {
1523                                         x[i] -= xj * cmplx.Conj(a[j*lda+i])
1524                                 }
1525                         }
1526                 } else {
1527                         jx := kx + (n-1)*incX
1528                         for j := n - 1; j >= 0; j-- {
1529                                 if diag == blas.NonUnit {
1530                                         x[jx] /= cmplx.Conj(a[j*lda+j])
1531                                 }
1532                                 xj := x[jx]
1533                                 ix := kx
1534                                 for i := 0; i < j; i++ {
1535                                         x[ix] -= xj * cmplx.Conj(a[j*lda+i])
1536                                         ix += incX
1537                                 }
1538                                 jx -= incX
1539                         }
1540                 }
1541         }
1542 }