github.com/jingcheng-WU/gonum@v0.9.1-0.20210323123734-f1a2a11a8f7b/blas/gonum/level2cmplx128.go (about) 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 "github.com/jingcheng-WU/gonum/blas" 11 "github.com/jingcheng-WU/gonum/internal/asm/c128" 12 ) 13 14 var _ blas.Complex128Level2 = Implementation{} 15 16 // Zgbmv performs one of the matrix-vector operations 17 // y = alpha * A * x + beta * y if trans = blas.NoTrans 18 // y = alpha * Aᵀ * x + beta * y if trans = blas.Trans 19 // y = alpha * Aᴴ * x + beta * y if trans = blas.ConjTrans 20 // where alpha and beta are scalars, x and y are vectors, and A is an m×n band matrix 21 // with kL sub-diagonals and kU super-diagonals. 22 func (Implementation) Zgbmv(trans blas.Transpose, m, n, kL, kU int, alpha complex128, a []complex128, lda int, x []complex128, incX int, beta complex128, y []complex128, incY int) { 23 switch trans { 24 default: 25 panic(badTranspose) 26 case blas.NoTrans, blas.Trans, blas.ConjTrans: 27 } 28 if m < 0 { 29 panic(mLT0) 30 } 31 if n < 0 { 32 panic(nLT0) 33 } 34 if kL < 0 { 35 panic(kLLT0) 36 } 37 if kU < 0 { 38 panic(kULT0) 39 } 40 if lda < kL+kU+1 { 41 panic(badLdA) 42 } 43 if incX == 0 { 44 panic(zeroIncX) 45 } 46 if incY == 0 { 47 panic(zeroIncY) 48 } 49 50 // Quick return if possible. 51 if m == 0 || n == 0 { 52 return 53 } 54 55 // For zero matrix size the following slice length checks are trivially satisfied. 56 if len(a) < lda*(min(m, n+kL)-1)+kL+kU+1 { 57 panic(shortA) 58 } 59 var lenX, lenY int 60 if trans == blas.NoTrans { 61 lenX, lenY = n, m 62 } else { 63 lenX, lenY = m, n 64 } 65 if (incX > 0 && len(x) <= (lenX-1)*incX) || (incX < 0 && len(x) <= (1-lenX)*incX) { 66 panic(shortX) 67 } 68 if (incY > 0 && len(y) <= (lenY-1)*incY) || (incY < 0 && len(y) <= (1-lenY)*incY) { 69 panic(shortY) 70 } 71 72 // Quick return if possible. 73 if alpha == 0 && beta == 1 { 74 return 75 } 76 77 var kx int 78 if incX < 0 { 79 kx = (1 - lenX) * incX 80 } 81 var ky int 82 if incY < 0 { 83 ky = (1 - lenY) * incY 84 } 85 86 // Form y = beta*y. 87 if beta != 1 { 88 if incY == 1 { 89 if beta == 0 { 90 for i := range y[:lenY] { 91 y[i] = 0 92 } 93 } else { 94 c128.ScalUnitary(beta, y[:lenY]) 95 } 96 } else { 97 iy := ky 98 if beta == 0 { 99 for i := 0; i < lenY; i++ { 100 y[iy] = 0 101 iy += incY 102 } 103 } else { 104 if incY > 0 { 105 c128.ScalInc(beta, y, uintptr(lenY), uintptr(incY)) 106 } else { 107 c128.ScalInc(beta, y, uintptr(lenY), uintptr(-incY)) 108 } 109 } 110 } 111 } 112 113 nRow := min(m, n+kL) 114 nCol := kL + 1 + kU 115 switch trans { 116 case blas.NoTrans: 117 iy := ky 118 if incX == 1 { 119 for i := 0; i < nRow; i++ { 120 l := max(0, kL-i) 121 u := min(nCol, n+kL-i) 122 aRow := a[i*lda+l : i*lda+u] 123 off := max(0, i-kL) 124 xtmp := x[off : off+u-l] 125 var sum complex128 126 for j, v := range aRow { 127 sum += xtmp[j] * v 128 } 129 y[iy] += alpha * sum 130 iy += incY 131 } 132 } else { 133 for i := 0; i < nRow; i++ { 134 l := max(0, kL-i) 135 u := min(nCol, n+kL-i) 136 aRow := a[i*lda+l : i*lda+u] 137 off := max(0, i-kL) * incX 138 jx := kx 139 var sum complex128 140 for _, v := range aRow { 141 sum += x[off+jx] * v 142 jx += incX 143 } 144 y[iy] += alpha * sum 145 iy += incY 146 } 147 } 148 case blas.Trans: 149 if incX == 1 { 150 for i := 0; i < nRow; i++ { 151 l := max(0, kL-i) 152 u := min(nCol, n+kL-i) 153 aRow := a[i*lda+l : i*lda+u] 154 off := max(0, i-kL) * incY 155 alphaxi := alpha * x[i] 156 jy := ky 157 for _, v := range aRow { 158 y[off+jy] += alphaxi * v 159 jy += incY 160 } 161 } 162 } else { 163 ix := kx 164 for i := 0; i < nRow; i++ { 165 l := max(0, kL-i) 166 u := min(nCol, n+kL-i) 167 aRow := a[i*lda+l : i*lda+u] 168 off := max(0, i-kL) * incY 169 alphaxi := alpha * x[ix] 170 jy := ky 171 for _, v := range aRow { 172 y[off+jy] += alphaxi * v 173 jy += incY 174 } 175 ix += incX 176 } 177 } 178 case blas.ConjTrans: 179 if incX == 1 { 180 for i := 0; i < nRow; i++ { 181 l := max(0, kL-i) 182 u := min(nCol, n+kL-i) 183 aRow := a[i*lda+l : i*lda+u] 184 off := max(0, i-kL) * incY 185 alphaxi := alpha * x[i] 186 jy := ky 187 for _, v := range aRow { 188 y[off+jy] += alphaxi * cmplx.Conj(v) 189 jy += incY 190 } 191 } 192 } else { 193 ix := kx 194 for i := 0; i < nRow; i++ { 195 l := max(0, kL-i) 196 u := min(nCol, n+kL-i) 197 aRow := a[i*lda+l : i*lda+u] 198 off := max(0, i-kL) * incY 199 alphaxi := alpha * x[ix] 200 jy := ky 201 for _, v := range aRow { 202 y[off+jy] += alphaxi * cmplx.Conj(v) 203 jy += incY 204 } 205 ix += incX 206 } 207 } 208 } 209 } 210 211 // Zgemv performs one of the matrix-vector operations 212 // y = alpha * A * x + beta * y if trans = blas.NoTrans 213 // y = alpha * Aᵀ * x + beta * y if trans = blas.Trans 214 // y = alpha * Aᴴ * x + beta * y if trans = blas.ConjTrans 215 // where alpha and beta are scalars, x and y are vectors, and A is an m×n dense matrix. 216 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) { 217 switch trans { 218 default: 219 panic(badTranspose) 220 case blas.NoTrans, blas.Trans, blas.ConjTrans: 221 } 222 if m < 0 { 223 panic(mLT0) 224 } 225 if n < 0 { 226 panic(nLT0) 227 } 228 if lda < max(1, n) { 229 panic(badLdA) 230 } 231 if incX == 0 { 232 panic(zeroIncX) 233 } 234 if incY == 0 { 235 panic(zeroIncY) 236 } 237 238 // Quick return if possible. 239 if m == 0 || n == 0 { 240 return 241 } 242 243 // For zero matrix size the following slice length checks are trivially satisfied. 244 var lenX, lenY int 245 if trans == blas.NoTrans { 246 lenX = n 247 lenY = m 248 } else { 249 lenX = m 250 lenY = n 251 } 252 if len(a) < lda*(m-1)+n { 253 panic(shortA) 254 } 255 if (incX > 0 && len(x) <= (lenX-1)*incX) || (incX < 0 && len(x) <= (1-lenX)*incX) { 256 panic(shortX) 257 } 258 if (incY > 0 && len(y) <= (lenY-1)*incY) || (incY < 0 && len(y) <= (1-lenY)*incY) { 259 panic(shortY) 260 } 261 262 // Quick return if possible. 263 if alpha == 0 && beta == 1 { 264 return 265 } 266 267 var kx int 268 if incX < 0 { 269 kx = (1 - lenX) * incX 270 } 271 var ky int 272 if incY < 0 { 273 ky = (1 - lenY) * incY 274 } 275 276 // Form y = beta*y. 277 if beta != 1 { 278 if incY == 1 { 279 if beta == 0 { 280 for i := range y[:lenY] { 281 y[i] = 0 282 } 283 } else { 284 c128.ScalUnitary(beta, y[:lenY]) 285 } 286 } else { 287 iy := ky 288 if beta == 0 { 289 for i := 0; i < lenY; i++ { 290 y[iy] = 0 291 iy += incY 292 } 293 } else { 294 if incY > 0 { 295 c128.ScalInc(beta, y, uintptr(lenY), uintptr(incY)) 296 } else { 297 c128.ScalInc(beta, y, uintptr(lenY), uintptr(-incY)) 298 } 299 } 300 } 301 } 302 303 if alpha == 0 { 304 return 305 } 306 307 switch trans { 308 default: 309 // Form y = alpha*A*x + y. 310 iy := ky 311 if incX == 1 { 312 for i := 0; i < m; i++ { 313 y[iy] += alpha * c128.DotuUnitary(a[i*lda:i*lda+n], x[:n]) 314 iy += incY 315 } 316 return 317 } 318 for i := 0; i < m; i++ { 319 y[iy] += alpha * c128.DotuInc(a[i*lda:i*lda+n], x, uintptr(n), 1, uintptr(incX), 0, uintptr(kx)) 320 iy += incY 321 } 322 return 323 324 case blas.Trans: 325 // Form y = alpha*Aᵀ*x + y. 326 ix := kx 327 if incY == 1 { 328 for i := 0; i < m; i++ { 329 c128.AxpyUnitary(alpha*x[ix], a[i*lda:i*lda+n], y[:n]) 330 ix += incX 331 } 332 return 333 } 334 for i := 0; i < m; i++ { 335 c128.AxpyInc(alpha*x[ix], a[i*lda:i*lda+n], y, uintptr(n), 1, uintptr(incY), 0, uintptr(ky)) 336 ix += incX 337 } 338 return 339 340 case blas.ConjTrans: 341 // Form y = alpha*Aᴴ*x + y. 342 ix := kx 343 if incY == 1 { 344 for i := 0; i < m; i++ { 345 tmp := alpha * x[ix] 346 for j := 0; j < n; j++ { 347 y[j] += tmp * cmplx.Conj(a[i*lda+j]) 348 } 349 ix += incX 350 } 351 return 352 } 353 for i := 0; i < m; i++ { 354 tmp := alpha * x[ix] 355 jy := ky 356 for j := 0; j < n; j++ { 357 y[jy] += tmp * cmplx.Conj(a[i*lda+j]) 358 jy += incY 359 } 360 ix += incX 361 } 362 return 363 } 364 } 365 366 // Zgerc performs the rank-one operation 367 // A += alpha * x * yᴴ 368 // where A is an m×n dense matrix, alpha is a scalar, x is an m element vector, 369 // and y is an n element vector. 370 func (Implementation) Zgerc(m, n int, alpha complex128, x []complex128, incX int, y []complex128, incY int, a []complex128, lda int) { 371 if m < 0 { 372 panic(mLT0) 373 } 374 if n < 0 { 375 panic(nLT0) 376 } 377 if lda < max(1, n) { 378 panic(badLdA) 379 } 380 if incX == 0 { 381 panic(zeroIncX) 382 } 383 if incY == 0 { 384 panic(zeroIncY) 385 } 386 387 // Quick return if possible. 388 if m == 0 || n == 0 { 389 return 390 } 391 392 // For zero matrix size the following slice length checks are trivially satisfied. 393 if (incX > 0 && len(x) <= (m-1)*incX) || (incX < 0 && len(x) <= (1-m)*incX) { 394 panic(shortX) 395 } 396 if (incY > 0 && len(y) <= (n-1)*incY) || (incY < 0 && len(y) <= (1-n)*incY) { 397 panic(shortY) 398 } 399 if len(a) < lda*(m-1)+n { 400 panic(shortA) 401 } 402 403 // Quick return if possible. 404 if alpha == 0 { 405 return 406 } 407 408 var kx, jy int 409 if incX < 0 { 410 kx = (1 - m) * incX 411 } 412 if incY < 0 { 413 jy = (1 - n) * incY 414 } 415 for j := 0; j < n; j++ { 416 if y[jy] != 0 { 417 tmp := alpha * cmplx.Conj(y[jy]) 418 c128.AxpyInc(tmp, x, a[j:], uintptr(m), uintptr(incX), uintptr(lda), uintptr(kx), 0) 419 } 420 jy += incY 421 } 422 } 423 424 // Zgeru performs the rank-one operation 425 // A += alpha * x * yᵀ 426 // where A is an m×n dense matrix, alpha is a scalar, x is an m element vector, 427 // and y is an n element vector. 428 func (Implementation) Zgeru(m, n int, alpha complex128, x []complex128, incX int, y []complex128, incY int, a []complex128, lda int) { 429 if m < 0 { 430 panic(mLT0) 431 } 432 if n < 0 { 433 panic(nLT0) 434 } 435 if lda < max(1, n) { 436 panic(badLdA) 437 } 438 if incX == 0 { 439 panic(zeroIncX) 440 } 441 if incY == 0 { 442 panic(zeroIncY) 443 } 444 445 // Quick return if possible. 446 if m == 0 || n == 0 { 447 return 448 } 449 450 // For zero matrix size the following slice length checks are trivially satisfied. 451 if (incX > 0 && len(x) <= (m-1)*incX) || (incX < 0 && len(x) <= (1-m)*incX) { 452 panic(shortX) 453 } 454 if (incY > 0 && len(y) <= (n-1)*incY) || (incY < 0 && len(y) <= (1-n)*incY) { 455 panic(shortY) 456 } 457 if len(a) < lda*(m-1)+n { 458 panic(shortA) 459 } 460 461 // Quick return if possible. 462 if alpha == 0 { 463 return 464 } 465 466 var kx int 467 if incX < 0 { 468 kx = (1 - m) * incX 469 } 470 if incY == 1 { 471 for i := 0; i < m; i++ { 472 if x[kx] != 0 { 473 tmp := alpha * x[kx] 474 c128.AxpyUnitary(tmp, y[:n], a[i*lda:i*lda+n]) 475 } 476 kx += incX 477 } 478 return 479 } 480 var jy int 481 if incY < 0 { 482 jy = (1 - n) * incY 483 } 484 for i := 0; i < m; i++ { 485 if x[kx] != 0 { 486 tmp := alpha * x[kx] 487 c128.AxpyInc(tmp, y, a[i*lda:i*lda+n], uintptr(n), uintptr(incY), 1, uintptr(jy), 0) 488 } 489 kx += incX 490 } 491 } 492 493 // Zhbmv performs the matrix-vector operation 494 // y = alpha * A * x + beta * y 495 // where alpha and beta are scalars, x and y are vectors, and A is an n×n 496 // Hermitian band matrix with k super-diagonals. The imaginary parts of 497 // the diagonal elements of A are ignored and assumed to be zero. 498 func (Implementation) Zhbmv(uplo blas.Uplo, n, k int, alpha complex128, a []complex128, lda int, x []complex128, incX int, beta complex128, y []complex128, incY int) { 499 switch uplo { 500 default: 501 panic(badUplo) 502 case blas.Upper, blas.Lower: 503 } 504 if n < 0 { 505 panic(nLT0) 506 } 507 if k < 0 { 508 panic(kLT0) 509 } 510 if lda < k+1 { 511 panic(badLdA) 512 } 513 if incX == 0 { 514 panic(zeroIncX) 515 } 516 if incY == 0 { 517 panic(zeroIncY) 518 } 519 520 // Quick return if possible. 521 if n == 0 { 522 return 523 } 524 525 // For zero matrix size the following slice length checks are trivially satisfied. 526 if len(a) < lda*(n-1)+k+1 { 527 panic(shortA) 528 } 529 if (incX > 0 && len(x) <= (n-1)*incX) || (incX < 0 && len(x) <= (1-n)*incX) { 530 panic(shortX) 531 } 532 if (incY > 0 && len(y) <= (n-1)*incY) || (incY < 0 && len(y) <= (1-n)*incY) { 533 panic(shortY) 534 } 535 536 // Quick return if possible. 537 if alpha == 0 && beta == 1 { 538 return 539 } 540 541 // Set up the start indices in X and Y. 542 var kx int 543 if incX < 0 { 544 kx = (1 - n) * incX 545 } 546 var ky int 547 if incY < 0 { 548 ky = (1 - n) * incY 549 } 550 551 // Form y = beta*y. 552 if beta != 1 { 553 if incY == 1 { 554 if beta == 0 { 555 for i := range y[:n] { 556 y[i] = 0 557 } 558 } else { 559 for i, v := range y[:n] { 560 y[i] = beta * v 561 } 562 } 563 } else { 564 iy := ky 565 if beta == 0 { 566 for i := 0; i < n; i++ { 567 y[iy] = 0 568 iy += incY 569 } 570 } else { 571 for i := 0; i < n; i++ { 572 y[iy] = beta * y[iy] 573 iy += incY 574 } 575 } 576 } 577 } 578 579 if alpha == 0 { 580 return 581 } 582 583 // The elements of A are accessed sequentially with one pass through a. 584 switch uplo { 585 case blas.Upper: 586 iy := ky 587 if incX == 1 { 588 for i := 0; i < n; i++ { 589 aRow := a[i*lda:] 590 alphaxi := alpha * x[i] 591 sum := alphaxi * complex(real(aRow[0]), 0) 592 u := min(k+1, n-i) 593 jy := incY 594 for j := 1; j < u; j++ { 595 v := aRow[j] 596 sum += alpha * x[i+j] * v 597 y[iy+jy] += alphaxi * cmplx.Conj(v) 598 jy += incY 599 } 600 y[iy] += sum 601 iy += incY 602 } 603 } else { 604 ix := kx 605 for i := 0; i < n; i++ { 606 aRow := a[i*lda:] 607 alphaxi := alpha * x[ix] 608 sum := alphaxi * complex(real(aRow[0]), 0) 609 u := min(k+1, n-i) 610 jx := incX 611 jy := incY 612 for j := 1; j < u; j++ { 613 v := aRow[j] 614 sum += alpha * x[ix+jx] * v 615 y[iy+jy] += alphaxi * cmplx.Conj(v) 616 jx += incX 617 jy += incY 618 } 619 y[iy] += sum 620 ix += incX 621 iy += incY 622 } 623 } 624 case blas.Lower: 625 iy := ky 626 if incX == 1 { 627 for i := 0; i < n; i++ { 628 l := max(0, k-i) 629 alphaxi := alpha * x[i] 630 jy := l * incY 631 aRow := a[i*lda:] 632 for j := l; j < k; j++ { 633 v := aRow[j] 634 y[iy] += alpha * v * x[i-k+j] 635 y[iy-k*incY+jy] += alphaxi * cmplx.Conj(v) 636 jy += incY 637 } 638 y[iy] += alphaxi * complex(real(aRow[k]), 0) 639 iy += incY 640 } 641 } else { 642 ix := kx 643 for i := 0; i < n; i++ { 644 l := max(0, k-i) 645 alphaxi := alpha * x[ix] 646 jx := l * incX 647 jy := l * incY 648 aRow := a[i*lda:] 649 for j := l; j < k; j++ { 650 v := aRow[j] 651 y[iy] += alpha * v * x[ix-k*incX+jx] 652 y[iy-k*incY+jy] += alphaxi * cmplx.Conj(v) 653 jx += incX 654 jy += incY 655 } 656 y[iy] += alphaxi * complex(real(aRow[k]), 0) 657 ix += incX 658 iy += incY 659 } 660 } 661 } 662 } 663 664 // Zhemv performs the matrix-vector operation 665 // y = alpha * A * x + beta * y 666 // where alpha and beta are scalars, x and y are vectors, and A is an n×n 667 // Hermitian matrix. The imaginary parts of the diagonal elements of A are 668 // ignored and assumed to be zero. 669 func (Implementation) Zhemv(uplo blas.Uplo, n int, alpha complex128, a []complex128, lda int, x []complex128, incX int, beta complex128, y []complex128, incY int) { 670 switch uplo { 671 default: 672 panic(badUplo) 673 case blas.Upper, blas.Lower: 674 } 675 if n < 0 { 676 panic(nLT0) 677 } 678 if lda < max(1, n) { 679 panic(badLdA) 680 } 681 if incX == 0 { 682 panic(zeroIncX) 683 } 684 if incY == 0 { 685 panic(zeroIncY) 686 } 687 688 // Quick return if possible. 689 if n == 0 { 690 return 691 } 692 693 // For zero matrix size the following slice length checks are trivially satisfied. 694 if len(a) < lda*(n-1)+n { 695 panic(shortA) 696 } 697 if (incX > 0 && len(x) <= (n-1)*incX) || (incX < 0 && len(x) <= (1-n)*incX) { 698 panic(shortX) 699 } 700 if (incY > 0 && len(y) <= (n-1)*incY) || (incY < 0 && len(y) <= (1-n)*incY) { 701 panic(shortY) 702 } 703 704 // Quick return if possible. 705 if alpha == 0 && beta == 1 { 706 return 707 } 708 709 // Set up the start indices in X and Y. 710 var kx int 711 if incX < 0 { 712 kx = (1 - n) * incX 713 } 714 var ky int 715 if incY < 0 { 716 ky = (1 - n) * incY 717 } 718 719 // Form y = beta*y. 720 if beta != 1 { 721 if incY == 1 { 722 if beta == 0 { 723 for i := range y[:n] { 724 y[i] = 0 725 } 726 } else { 727 for i, v := range y[:n] { 728 y[i] = beta * v 729 } 730 } 731 } else { 732 iy := ky 733 if beta == 0 { 734 for i := 0; i < n; i++ { 735 y[iy] = 0 736 iy += incY 737 } 738 } else { 739 for i := 0; i < n; i++ { 740 y[iy] = beta * y[iy] 741 iy += incY 742 } 743 } 744 } 745 } 746 747 if alpha == 0 { 748 return 749 } 750 751 // The elements of A are accessed sequentially with one pass through 752 // the triangular part of A. 753 754 if uplo == blas.Upper { 755 // Form y when A is stored in upper triangle. 756 if incX == 1 && incY == 1 { 757 for i := 0; i < n; i++ { 758 tmp1 := alpha * x[i] 759 var tmp2 complex128 760 for j := i + 1; j < n; j++ { 761 y[j] += tmp1 * cmplx.Conj(a[i*lda+j]) 762 tmp2 += a[i*lda+j] * x[j] 763 } 764 aii := complex(real(a[i*lda+i]), 0) 765 y[i] += tmp1*aii + alpha*tmp2 766 } 767 } else { 768 ix := kx 769 iy := ky 770 for i := 0; i < n; i++ { 771 tmp1 := alpha * x[ix] 772 var tmp2 complex128 773 jx := ix 774 jy := iy 775 for j := i + 1; j < n; j++ { 776 jx += incX 777 jy += incY 778 y[jy] += tmp1 * cmplx.Conj(a[i*lda+j]) 779 tmp2 += a[i*lda+j] * x[jx] 780 } 781 aii := complex(real(a[i*lda+i]), 0) 782 y[iy] += tmp1*aii + alpha*tmp2 783 ix += incX 784 iy += incY 785 } 786 } 787 return 788 } 789 790 // Form y when A is stored in lower triangle. 791 if incX == 1 && incY == 1 { 792 for i := 0; i < n; i++ { 793 tmp1 := alpha * x[i] 794 var tmp2 complex128 795 for j := 0; j < i; j++ { 796 y[j] += tmp1 * cmplx.Conj(a[i*lda+j]) 797 tmp2 += a[i*lda+j] * x[j] 798 } 799 aii := complex(real(a[i*lda+i]), 0) 800 y[i] += tmp1*aii + alpha*tmp2 801 } 802 } else { 803 ix := kx 804 iy := ky 805 for i := 0; i < n; i++ { 806 tmp1 := alpha * x[ix] 807 var tmp2 complex128 808 jx := kx 809 jy := ky 810 for j := 0; j < i; j++ { 811 y[jy] += tmp1 * cmplx.Conj(a[i*lda+j]) 812 tmp2 += a[i*lda+j] * x[jx] 813 jx += incX 814 jy += incY 815 } 816 aii := complex(real(a[i*lda+i]), 0) 817 y[iy] += tmp1*aii + alpha*tmp2 818 ix += incX 819 iy += incY 820 } 821 } 822 } 823 824 // Zher performs the Hermitian rank-one operation 825 // A += alpha * x * xᴴ 826 // where A is an n×n Hermitian matrix, alpha is a real scalar, and x is an n 827 // element vector. On entry, the imaginary parts of the diagonal elements of A 828 // are ignored and assumed to be zero, on return they will be set to zero. 829 func (Implementation) Zher(uplo blas.Uplo, n int, alpha float64, x []complex128, incX int, a []complex128, lda int) { 830 switch uplo { 831 default: 832 panic(badUplo) 833 case blas.Upper, blas.Lower: 834 } 835 if n < 0 { 836 panic(nLT0) 837 } 838 if lda < max(1, n) { 839 panic(badLdA) 840 } 841 if incX == 0 { 842 panic(zeroIncX) 843 } 844 845 // Quick return if possible. 846 if n == 0 { 847 return 848 } 849 850 // For zero matrix size the following slice length checks are trivially satisfied. 851 if (incX > 0 && len(x) <= (n-1)*incX) || (incX < 0 && len(x) <= (1-n)*incX) { 852 panic(shortX) 853 } 854 if len(a) < lda*(n-1)+n { 855 panic(shortA) 856 } 857 858 // Quick return if possible. 859 if alpha == 0 { 860 return 861 } 862 863 var kx int 864 if incX < 0 { 865 kx = (1 - n) * incX 866 } 867 if uplo == blas.Upper { 868 if incX == 1 { 869 for i := 0; i < n; i++ { 870 if x[i] != 0 { 871 tmp := complex(alpha*real(x[i]), alpha*imag(x[i])) 872 aii := real(a[i*lda+i]) 873 xtmp := real(tmp * cmplx.Conj(x[i])) 874 a[i*lda+i] = complex(aii+xtmp, 0) 875 for j := i + 1; j < n; j++ { 876 a[i*lda+j] += tmp * cmplx.Conj(x[j]) 877 } 878 } else { 879 aii := real(a[i*lda+i]) 880 a[i*lda+i] = complex(aii, 0) 881 } 882 } 883 return 884 } 885 886 ix := kx 887 for i := 0; i < n; i++ { 888 if x[ix] != 0 { 889 tmp := complex(alpha*real(x[ix]), alpha*imag(x[ix])) 890 aii := real(a[i*lda+i]) 891 xtmp := real(tmp * cmplx.Conj(x[ix])) 892 a[i*lda+i] = complex(aii+xtmp, 0) 893 jx := ix + incX 894 for j := i + 1; j < n; j++ { 895 a[i*lda+j] += tmp * cmplx.Conj(x[jx]) 896 jx += incX 897 } 898 } else { 899 aii := real(a[i*lda+i]) 900 a[i*lda+i] = complex(aii, 0) 901 } 902 ix += incX 903 } 904 return 905 } 906 907 if incX == 1 { 908 for i := 0; i < n; i++ { 909 if x[i] != 0 { 910 tmp := complex(alpha*real(x[i]), alpha*imag(x[i])) 911 for j := 0; j < i; j++ { 912 a[i*lda+j] += tmp * cmplx.Conj(x[j]) 913 } 914 aii := real(a[i*lda+i]) 915 xtmp := real(tmp * cmplx.Conj(x[i])) 916 a[i*lda+i] = complex(aii+xtmp, 0) 917 } else { 918 aii := real(a[i*lda+i]) 919 a[i*lda+i] = complex(aii, 0) 920 } 921 } 922 return 923 } 924 925 ix := kx 926 for i := 0; i < n; i++ { 927 if x[ix] != 0 { 928 tmp := complex(alpha*real(x[ix]), alpha*imag(x[ix])) 929 jx := kx 930 for j := 0; j < i; j++ { 931 a[i*lda+j] += tmp * cmplx.Conj(x[jx]) 932 jx += incX 933 } 934 aii := real(a[i*lda+i]) 935 xtmp := real(tmp * cmplx.Conj(x[ix])) 936 a[i*lda+i] = complex(aii+xtmp, 0) 937 938 } else { 939 aii := real(a[i*lda+i]) 940 a[i*lda+i] = complex(aii, 0) 941 } 942 ix += incX 943 } 944 } 945 946 // Zher2 performs the Hermitian rank-two operation 947 // A += alpha * x * yᴴ + conj(alpha) * y * xᴴ 948 // where alpha is a scalar, x and y are n element vectors and A is an n×n 949 // Hermitian matrix. On entry, the imaginary parts of the diagonal elements are 950 // ignored and assumed to be zero. On return they will be set to zero. 951 func (Implementation) Zher2(uplo blas.Uplo, n int, alpha complex128, x []complex128, incX int, y []complex128, incY int, a []complex128, lda int) { 952 switch uplo { 953 default: 954 panic(badUplo) 955 case blas.Upper, blas.Lower: 956 } 957 if n < 0 { 958 panic(nLT0) 959 } 960 if lda < max(1, n) { 961 panic(badLdA) 962 } 963 if incX == 0 { 964 panic(zeroIncX) 965 } 966 if incY == 0 { 967 panic(zeroIncY) 968 } 969 970 // Quick return if possible. 971 if n == 0 { 972 return 973 } 974 975 // For zero matrix size the following slice length checks are trivially satisfied. 976 if (incX > 0 && len(x) <= (n-1)*incX) || (incX < 0 && len(x) <= (1-n)*incX) { 977 panic(shortX) 978 } 979 if (incY > 0 && len(y) <= (n-1)*incY) || (incY < 0 && len(y) <= (1-n)*incY) { 980 panic(shortY) 981 } 982 if len(a) < lda*(n-1)+n { 983 panic(shortA) 984 } 985 986 // Quick return if possible. 987 if alpha == 0 { 988 return 989 } 990 991 var kx, ky int 992 var ix, iy int 993 if incX != 1 || incY != 1 { 994 if incX < 0 { 995 kx = (1 - n) * incX 996 } 997 if incY < 0 { 998 ky = (1 - n) * incY 999 } 1000 ix = kx 1001 iy = ky 1002 } 1003 if uplo == blas.Upper { 1004 if incX == 1 && incY == 1 { 1005 for i := 0; i < n; i++ { 1006 if x[i] != 0 || y[i] != 0 { 1007 tmp1 := alpha * x[i] 1008 tmp2 := cmplx.Conj(alpha) * y[i] 1009 aii := real(a[i*lda+i]) + real(tmp1*cmplx.Conj(y[i])) + real(tmp2*cmplx.Conj(x[i])) 1010 a[i*lda+i] = complex(aii, 0) 1011 for j := i + 1; j < n; j++ { 1012 a[i*lda+j] += tmp1*cmplx.Conj(y[j]) + tmp2*cmplx.Conj(x[j]) 1013 } 1014 } else { 1015 aii := real(a[i*lda+i]) 1016 a[i*lda+i] = complex(aii, 0) 1017 } 1018 } 1019 return 1020 } 1021 for i := 0; i < n; i++ { 1022 if x[ix] != 0 || y[iy] != 0 { 1023 tmp1 := alpha * x[ix] 1024 tmp2 := cmplx.Conj(alpha) * y[iy] 1025 aii := real(a[i*lda+i]) + real(tmp1*cmplx.Conj(y[iy])) + real(tmp2*cmplx.Conj(x[ix])) 1026 a[i*lda+i] = complex(aii, 0) 1027 jx := ix + incX 1028 jy := iy + incY 1029 for j := i + 1; j < n; j++ { 1030 a[i*lda+j] += tmp1*cmplx.Conj(y[jy]) + tmp2*cmplx.Conj(x[jx]) 1031 jx += incX 1032 jy += incY 1033 } 1034 } else { 1035 aii := real(a[i*lda+i]) 1036 a[i*lda+i] = complex(aii, 0) 1037 } 1038 ix += incX 1039 iy += incY 1040 } 1041 return 1042 } 1043 1044 if incX == 1 && incY == 1 { 1045 for i := 0; i < n; i++ { 1046 if x[i] != 0 || y[i] != 0 { 1047 tmp1 := alpha * x[i] 1048 tmp2 := cmplx.Conj(alpha) * y[i] 1049 for j := 0; j < i; j++ { 1050 a[i*lda+j] += tmp1*cmplx.Conj(y[j]) + tmp2*cmplx.Conj(x[j]) 1051 } 1052 aii := real(a[i*lda+i]) + real(tmp1*cmplx.Conj(y[i])) + real(tmp2*cmplx.Conj(x[i])) 1053 a[i*lda+i] = complex(aii, 0) 1054 } else { 1055 aii := real(a[i*lda+i]) 1056 a[i*lda+i] = complex(aii, 0) 1057 } 1058 } 1059 return 1060 } 1061 for i := 0; i < n; i++ { 1062 if x[ix] != 0 || y[iy] != 0 { 1063 tmp1 := alpha * x[ix] 1064 tmp2 := cmplx.Conj(alpha) * y[iy] 1065 jx := kx 1066 jy := ky 1067 for j := 0; j < i; j++ { 1068 a[i*lda+j] += tmp1*cmplx.Conj(y[jy]) + tmp2*cmplx.Conj(x[jx]) 1069 jx += incX 1070 jy += incY 1071 } 1072 aii := real(a[i*lda+i]) + real(tmp1*cmplx.Conj(y[iy])) + real(tmp2*cmplx.Conj(x[ix])) 1073 a[i*lda+i] = complex(aii, 0) 1074 } else { 1075 aii := real(a[i*lda+i]) 1076 a[i*lda+i] = complex(aii, 0) 1077 } 1078 ix += incX 1079 iy += incY 1080 } 1081 } 1082 1083 // Zhpmv performs the matrix-vector operation 1084 // y = alpha * A * x + beta * y 1085 // where alpha and beta are scalars, x and y are vectors, and A is an n×n 1086 // Hermitian matrix in packed form. The imaginary parts of the diagonal 1087 // elements of A are ignored and assumed to be zero. 1088 func (Implementation) Zhpmv(uplo blas.Uplo, n int, alpha complex128, ap []complex128, x []complex128, incX int, beta complex128, y []complex128, incY int) { 1089 switch uplo { 1090 default: 1091 panic(badUplo) 1092 case blas.Upper, blas.Lower: 1093 } 1094 if n < 0 { 1095 panic(nLT0) 1096 } 1097 if incX == 0 { 1098 panic(zeroIncX) 1099 } 1100 if incY == 0 { 1101 panic(zeroIncY) 1102 } 1103 1104 // Quick return if possible. 1105 if n == 0 { 1106 return 1107 } 1108 1109 // For zero matrix size the following slice length checks are trivially satisfied. 1110 if len(ap) < n*(n+1)/2 { 1111 panic(shortAP) 1112 } 1113 if (incX > 0 && len(x) <= (n-1)*incX) || (incX < 0 && len(x) <= (1-n)*incX) { 1114 panic(shortX) 1115 } 1116 if (incY > 0 && len(y) <= (n-1)*incY) || (incY < 0 && len(y) <= (1-n)*incY) { 1117 panic(shortY) 1118 } 1119 1120 // Quick return if possible. 1121 if alpha == 0 && beta == 1 { 1122 return 1123 } 1124 1125 // Set up the start indices in X and Y. 1126 var kx int 1127 if incX < 0 { 1128 kx = (1 - n) * incX 1129 } 1130 var ky int 1131 if incY < 0 { 1132 ky = (1 - n) * incY 1133 } 1134 1135 // Form y = beta*y. 1136 if beta != 1 { 1137 if incY == 1 { 1138 if beta == 0 { 1139 for i := range y[:n] { 1140 y[i] = 0 1141 } 1142 } else { 1143 for i, v := range y[:n] { 1144 y[i] = beta * v 1145 } 1146 } 1147 } else { 1148 iy := ky 1149 if beta == 0 { 1150 for i := 0; i < n; i++ { 1151 y[iy] = 0 1152 iy += incY 1153 } 1154 } else { 1155 for i := 0; i < n; i++ { 1156 y[iy] *= beta 1157 iy += incY 1158 } 1159 } 1160 } 1161 } 1162 1163 if alpha == 0 { 1164 return 1165 } 1166 1167 // The elements of A are accessed sequentially with one pass through ap. 1168 1169 var kk int 1170 if uplo == blas.Upper { 1171 // Form y when ap contains the upper triangle. 1172 // Here, kk points to the current diagonal element in ap. 1173 if incX == 1 && incY == 1 { 1174 for i := 0; i < n; i++ { 1175 tmp1 := alpha * x[i] 1176 y[i] += tmp1 * complex(real(ap[kk]), 0) 1177 var tmp2 complex128 1178 k := kk + 1 1179 for j := i + 1; j < n; j++ { 1180 y[j] += tmp1 * cmplx.Conj(ap[k]) 1181 tmp2 += ap[k] * x[j] 1182 k++ 1183 } 1184 y[i] += alpha * tmp2 1185 kk += n - i 1186 } 1187 } else { 1188 ix := kx 1189 iy := ky 1190 for i := 0; i < n; i++ { 1191 tmp1 := alpha * x[ix] 1192 y[iy] += tmp1 * complex(real(ap[kk]), 0) 1193 var tmp2 complex128 1194 jx := ix 1195 jy := iy 1196 for k := kk + 1; k < kk+n-i; k++ { 1197 jx += incX 1198 jy += incY 1199 y[jy] += tmp1 * cmplx.Conj(ap[k]) 1200 tmp2 += ap[k] * x[jx] 1201 } 1202 y[iy] += alpha * tmp2 1203 ix += incX 1204 iy += incY 1205 kk += n - i 1206 } 1207 } 1208 return 1209 } 1210 1211 // Form y when ap contains the lower triangle. 1212 // Here, kk points to the beginning of current row in ap. 1213 if incX == 1 && incY == 1 { 1214 for i := 0; i < n; i++ { 1215 tmp1 := alpha * x[i] 1216 var tmp2 complex128 1217 k := kk 1218 for j := 0; j < i; j++ { 1219 y[j] += tmp1 * cmplx.Conj(ap[k]) 1220 tmp2 += ap[k] * x[j] 1221 k++ 1222 } 1223 aii := complex(real(ap[kk+i]), 0) 1224 y[i] += tmp1*aii + alpha*tmp2 1225 kk += i + 1 1226 } 1227 } else { 1228 ix := kx 1229 iy := ky 1230 for i := 0; i < n; i++ { 1231 tmp1 := alpha * x[ix] 1232 var tmp2 complex128 1233 jx := kx 1234 jy := ky 1235 for k := kk; k < kk+i; k++ { 1236 y[jy] += tmp1 * cmplx.Conj(ap[k]) 1237 tmp2 += ap[k] * x[jx] 1238 jx += incX 1239 jy += incY 1240 } 1241 aii := complex(real(ap[kk+i]), 0) 1242 y[iy] += tmp1*aii + alpha*tmp2 1243 ix += incX 1244 iy += incY 1245 kk += i + 1 1246 } 1247 } 1248 } 1249 1250 // Zhpr performs the Hermitian rank-1 operation 1251 // A += alpha * x * xᴴ 1252 // where alpha is a real scalar, x is a vector, and A is an n×n hermitian matrix 1253 // in packed form. On entry, the imaginary parts of the diagonal elements are 1254 // assumed to be zero, and on return they are set to zero. 1255 func (Implementation) Zhpr(uplo blas.Uplo, n int, alpha float64, x []complex128, incX int, ap []complex128) { 1256 switch uplo { 1257 default: 1258 panic(badUplo) 1259 case blas.Upper, blas.Lower: 1260 } 1261 if n < 0 { 1262 panic(nLT0) 1263 } 1264 if incX == 0 { 1265 panic(zeroIncX) 1266 } 1267 1268 // Quick return if possible. 1269 if n == 0 { 1270 return 1271 } 1272 1273 // For zero matrix size the following slice length checks are trivially satisfied. 1274 if (incX > 0 && len(x) <= (n-1)*incX) || (incX < 0 && len(x) <= (1-n)*incX) { 1275 panic(shortX) 1276 } 1277 if len(ap) < n*(n+1)/2 { 1278 panic(shortAP) 1279 } 1280 1281 // Quick return if possible. 1282 if alpha == 0 { 1283 return 1284 } 1285 1286 // Set up start index in X. 1287 var kx int 1288 if incX < 0 { 1289 kx = (1 - n) * incX 1290 } 1291 1292 // The elements of A are accessed sequentially with one pass through ap. 1293 1294 var kk int 1295 if uplo == blas.Upper { 1296 // Form A when upper triangle is stored in AP. 1297 // Here, kk points to the current diagonal element in ap. 1298 if incX == 1 { 1299 for i := 0; i < n; i++ { 1300 xi := x[i] 1301 if xi != 0 { 1302 aii := real(ap[kk]) + alpha*real(cmplx.Conj(xi)*xi) 1303 ap[kk] = complex(aii, 0) 1304 1305 tmp := complex(alpha, 0) * xi 1306 a := ap[kk+1 : kk+n-i] 1307 x := x[i+1 : n] 1308 for j, v := range x { 1309 a[j] += tmp * cmplx.Conj(v) 1310 } 1311 } else { 1312 ap[kk] = complex(real(ap[kk]), 0) 1313 } 1314 kk += n - i 1315 } 1316 } else { 1317 ix := kx 1318 for i := 0; i < n; i++ { 1319 xi := x[ix] 1320 if xi != 0 { 1321 aii := real(ap[kk]) + alpha*real(cmplx.Conj(xi)*xi) 1322 ap[kk] = complex(aii, 0) 1323 1324 tmp := complex(alpha, 0) * xi 1325 jx := ix + incX 1326 a := ap[kk+1 : kk+n-i] 1327 for k := range a { 1328 a[k] += tmp * cmplx.Conj(x[jx]) 1329 jx += incX 1330 } 1331 } else { 1332 ap[kk] = complex(real(ap[kk]), 0) 1333 } 1334 ix += incX 1335 kk += n - i 1336 } 1337 } 1338 return 1339 } 1340 1341 // Form A when lower triangle is stored in AP. 1342 // Here, kk points to the beginning of current row in ap. 1343 if incX == 1 { 1344 for i := 0; i < n; i++ { 1345 xi := x[i] 1346 if xi != 0 { 1347 tmp := complex(alpha, 0) * xi 1348 a := ap[kk : kk+i] 1349 for j, v := range x[:i] { 1350 a[j] += tmp * cmplx.Conj(v) 1351 } 1352 1353 aii := real(ap[kk+i]) + alpha*real(cmplx.Conj(xi)*xi) 1354 ap[kk+i] = complex(aii, 0) 1355 } else { 1356 ap[kk+i] = complex(real(ap[kk+i]), 0) 1357 } 1358 kk += i + 1 1359 } 1360 } else { 1361 ix := kx 1362 for i := 0; i < n; i++ { 1363 xi := x[ix] 1364 if xi != 0 { 1365 tmp := complex(alpha, 0) * xi 1366 a := ap[kk : kk+i] 1367 jx := kx 1368 for k := range a { 1369 a[k] += tmp * cmplx.Conj(x[jx]) 1370 jx += incX 1371 } 1372 1373 aii := real(ap[kk+i]) + alpha*real(cmplx.Conj(xi)*xi) 1374 ap[kk+i] = complex(aii, 0) 1375 } else { 1376 ap[kk+i] = complex(real(ap[kk+i]), 0) 1377 } 1378 ix += incX 1379 kk += i + 1 1380 } 1381 } 1382 } 1383 1384 // Zhpr2 performs the Hermitian rank-2 operation 1385 // A += alpha * x * yᴴ + conj(alpha) * y * xᴴ 1386 // where alpha is a complex scalar, x and y are n element vectors, and A is an 1387 // n×n Hermitian matrix, supplied in packed form. On entry, the imaginary parts 1388 // of the diagonal elements are assumed to be zero, and on return they are set to zero. 1389 func (Implementation) Zhpr2(uplo blas.Uplo, n int, alpha complex128, x []complex128, incX int, y []complex128, incY int, ap []complex128) { 1390 switch uplo { 1391 default: 1392 panic(badUplo) 1393 case blas.Upper, blas.Lower: 1394 } 1395 if n < 0 { 1396 panic(nLT0) 1397 } 1398 if incX == 0 { 1399 panic(zeroIncX) 1400 } 1401 if incY == 0 { 1402 panic(zeroIncY) 1403 } 1404 1405 // Quick return if possible. 1406 if n == 0 { 1407 return 1408 } 1409 1410 // For zero matrix size the following slice length checks are trivially satisfied. 1411 if (incX > 0 && len(x) <= (n-1)*incX) || (incX < 0 && len(x) <= (1-n)*incX) { 1412 panic(shortX) 1413 } 1414 if (incY > 0 && len(y) <= (n-1)*incY) || (incY < 0 && len(y) <= (1-n)*incY) { 1415 panic(shortY) 1416 } 1417 if len(ap) < n*(n+1)/2 { 1418 panic(shortAP) 1419 } 1420 1421 // Quick return if possible. 1422 if alpha == 0 { 1423 return 1424 } 1425 1426 // Set up start indices in X and Y. 1427 var kx int 1428 if incX < 0 { 1429 kx = (1 - n) * incX 1430 } 1431 var ky int 1432 if incY < 0 { 1433 ky = (1 - n) * incY 1434 } 1435 1436 // The elements of A are accessed sequentially with one pass through ap. 1437 1438 var kk int 1439 if uplo == blas.Upper { 1440 // Form A when upper triangle is stored in AP. 1441 // Here, kk points to the current diagonal element in ap. 1442 if incX == 1 && incY == 1 { 1443 for i := 0; i < n; i++ { 1444 if x[i] != 0 || y[i] != 0 { 1445 tmp1 := alpha * x[i] 1446 tmp2 := cmplx.Conj(alpha) * y[i] 1447 aii := real(ap[kk]) + real(tmp1*cmplx.Conj(y[i])) + real(tmp2*cmplx.Conj(x[i])) 1448 ap[kk] = complex(aii, 0) 1449 k := kk + 1 1450 for j := i + 1; j < n; j++ { 1451 ap[k] += tmp1*cmplx.Conj(y[j]) + tmp2*cmplx.Conj(x[j]) 1452 k++ 1453 } 1454 } else { 1455 ap[kk] = complex(real(ap[kk]), 0) 1456 } 1457 kk += n - i 1458 } 1459 } else { 1460 ix := kx 1461 iy := ky 1462 for i := 0; i < n; i++ { 1463 if x[ix] != 0 || y[iy] != 0 { 1464 tmp1 := alpha * x[ix] 1465 tmp2 := cmplx.Conj(alpha) * y[iy] 1466 aii := real(ap[kk]) + real(tmp1*cmplx.Conj(y[iy])) + real(tmp2*cmplx.Conj(x[ix])) 1467 ap[kk] = complex(aii, 0) 1468 jx := ix + incX 1469 jy := iy + incY 1470 for k := kk + 1; k < kk+n-i; k++ { 1471 ap[k] += tmp1*cmplx.Conj(y[jy]) + tmp2*cmplx.Conj(x[jx]) 1472 jx += incX 1473 jy += incY 1474 } 1475 } else { 1476 ap[kk] = complex(real(ap[kk]), 0) 1477 } 1478 ix += incX 1479 iy += incY 1480 kk += n - i 1481 } 1482 } 1483 return 1484 } 1485 1486 // Form A when lower triangle is stored in AP. 1487 // Here, kk points to the beginning of current row in ap. 1488 if incX == 1 && incY == 1 { 1489 for i := 0; i < n; i++ { 1490 if x[i] != 0 || y[i] != 0 { 1491 tmp1 := alpha * x[i] 1492 tmp2 := cmplx.Conj(alpha) * y[i] 1493 k := kk 1494 for j := 0; j < i; j++ { 1495 ap[k] += tmp1*cmplx.Conj(y[j]) + tmp2*cmplx.Conj(x[j]) 1496 k++ 1497 } 1498 aii := real(ap[kk+i]) + real(tmp1*cmplx.Conj(y[i])) + real(tmp2*cmplx.Conj(x[i])) 1499 ap[kk+i] = complex(aii, 0) 1500 } else { 1501 ap[kk+i] = complex(real(ap[kk+i]), 0) 1502 } 1503 kk += i + 1 1504 } 1505 } else { 1506 ix := kx 1507 iy := ky 1508 for i := 0; i < n; i++ { 1509 if x[ix] != 0 || y[iy] != 0 { 1510 tmp1 := alpha * x[ix] 1511 tmp2 := cmplx.Conj(alpha) * y[iy] 1512 jx := kx 1513 jy := ky 1514 for k := kk; k < kk+i; k++ { 1515 ap[k] += tmp1*cmplx.Conj(y[jy]) + tmp2*cmplx.Conj(x[jx]) 1516 jx += incX 1517 jy += incY 1518 } 1519 aii := real(ap[kk+i]) + real(tmp1*cmplx.Conj(y[iy])) + real(tmp2*cmplx.Conj(x[ix])) 1520 ap[kk+i] = complex(aii, 0) 1521 } else { 1522 ap[kk+i] = complex(real(ap[kk+i]), 0) 1523 } 1524 ix += incX 1525 iy += incY 1526 kk += i + 1 1527 } 1528 } 1529 } 1530 1531 // Ztbmv performs one of the matrix-vector operations 1532 // x = A * x if trans = blas.NoTrans 1533 // x = Aᵀ * x if trans = blas.Trans 1534 // x = Aᴴ * x if trans = blas.ConjTrans 1535 // where x is an n element vector and A is an n×n triangular band matrix, with 1536 // (k+1) diagonals. 1537 func (Implementation) Ztbmv(uplo blas.Uplo, trans blas.Transpose, diag blas.Diag, n, k int, a []complex128, lda int, x []complex128, incX int) { 1538 switch trans { 1539 default: 1540 panic(badTranspose) 1541 case blas.NoTrans, blas.Trans, blas.ConjTrans: 1542 } 1543 switch uplo { 1544 default: 1545 panic(badUplo) 1546 case blas.Upper, blas.Lower: 1547 } 1548 switch diag { 1549 default: 1550 panic(badDiag) 1551 case blas.NonUnit, blas.Unit: 1552 } 1553 if n < 0 { 1554 panic(nLT0) 1555 } 1556 if k < 0 { 1557 panic(kLT0) 1558 } 1559 if lda < k+1 { 1560 panic(badLdA) 1561 } 1562 if incX == 0 { 1563 panic(zeroIncX) 1564 } 1565 1566 // Quick return if possible. 1567 if n == 0 { 1568 return 1569 } 1570 1571 // For zero matrix size the following slice length checks are trivially satisfied. 1572 if len(a) < lda*(n-1)+k+1 { 1573 panic(shortA) 1574 } 1575 if (incX > 0 && len(x) <= (n-1)*incX) || (incX < 0 && len(x) <= (1-n)*incX) { 1576 panic(shortX) 1577 } 1578 1579 // Set up start index in X. 1580 var kx int 1581 if incX < 0 { 1582 kx = (1 - n) * incX 1583 } 1584 1585 switch trans { 1586 case blas.NoTrans: 1587 if uplo == blas.Upper { 1588 if incX == 1 { 1589 for i := 0; i < n; i++ { 1590 xi := x[i] 1591 if diag == blas.NonUnit { 1592 xi *= a[i*lda] 1593 } 1594 kk := min(k, n-i-1) 1595 for j, aij := range a[i*lda+1 : i*lda+kk+1] { 1596 xi += x[i+j+1] * aij 1597 } 1598 x[i] = xi 1599 } 1600 } else { 1601 ix := kx 1602 for i := 0; i < n; i++ { 1603 xi := x[ix] 1604 if diag == blas.NonUnit { 1605 xi *= a[i*lda] 1606 } 1607 kk := min(k, n-i-1) 1608 jx := ix + incX 1609 for _, aij := range a[i*lda+1 : i*lda+kk+1] { 1610 xi += x[jx] * aij 1611 jx += incX 1612 } 1613 x[ix] = xi 1614 ix += incX 1615 } 1616 } 1617 } else { 1618 if incX == 1 { 1619 for i := n - 1; i >= 0; i-- { 1620 xi := x[i] 1621 if diag == blas.NonUnit { 1622 xi *= a[i*lda+k] 1623 } 1624 kk := min(k, i) 1625 for j, aij := range a[i*lda+k-kk : i*lda+k] { 1626 xi += x[i-kk+j] * aij 1627 } 1628 x[i] = xi 1629 } 1630 } else { 1631 ix := kx + (n-1)*incX 1632 for i := n - 1; i >= 0; i-- { 1633 xi := x[ix] 1634 if diag == blas.NonUnit { 1635 xi *= a[i*lda+k] 1636 } 1637 kk := min(k, i) 1638 jx := ix - kk*incX 1639 for _, aij := range a[i*lda+k-kk : i*lda+k] { 1640 xi += x[jx] * aij 1641 jx += incX 1642 } 1643 x[ix] = xi 1644 ix -= incX 1645 } 1646 } 1647 } 1648 case blas.Trans: 1649 if uplo == blas.Upper { 1650 if incX == 1 { 1651 for i := n - 1; i >= 0; i-- { 1652 kk := min(k, n-i-1) 1653 xi := x[i] 1654 for j, aij := range a[i*lda+1 : i*lda+kk+1] { 1655 x[i+j+1] += xi * aij 1656 } 1657 if diag == blas.NonUnit { 1658 x[i] *= a[i*lda] 1659 } 1660 } 1661 } else { 1662 ix := kx + (n-1)*incX 1663 for i := n - 1; i >= 0; i-- { 1664 kk := min(k, n-i-1) 1665 jx := ix + incX 1666 xi := x[ix] 1667 for _, aij := range a[i*lda+1 : i*lda+kk+1] { 1668 x[jx] += xi * aij 1669 jx += incX 1670 } 1671 if diag == blas.NonUnit { 1672 x[ix] *= a[i*lda] 1673 } 1674 ix -= incX 1675 } 1676 } 1677 } else { 1678 if incX == 1 { 1679 for i := 0; i < n; i++ { 1680 kk := min(k, i) 1681 xi := x[i] 1682 for j, aij := range a[i*lda+k-kk : i*lda+k] { 1683 x[i-kk+j] += xi * aij 1684 } 1685 if diag == blas.NonUnit { 1686 x[i] *= a[i*lda+k] 1687 } 1688 } 1689 } else { 1690 ix := kx 1691 for i := 0; i < n; i++ { 1692 kk := min(k, i) 1693 jx := ix - kk*incX 1694 xi := x[ix] 1695 for _, aij := range a[i*lda+k-kk : i*lda+k] { 1696 x[jx] += xi * aij 1697 jx += incX 1698 } 1699 if diag == blas.NonUnit { 1700 x[ix] *= a[i*lda+k] 1701 } 1702 ix += incX 1703 } 1704 } 1705 } 1706 case blas.ConjTrans: 1707 if uplo == blas.Upper { 1708 if incX == 1 { 1709 for i := n - 1; i >= 0; i-- { 1710 kk := min(k, n-i-1) 1711 xi := x[i] 1712 for j, aij := range a[i*lda+1 : i*lda+kk+1] { 1713 x[i+j+1] += xi * cmplx.Conj(aij) 1714 } 1715 if diag == blas.NonUnit { 1716 x[i] *= cmplx.Conj(a[i*lda]) 1717 } 1718 } 1719 } else { 1720 ix := kx + (n-1)*incX 1721 for i := n - 1; i >= 0; i-- { 1722 kk := min(k, n-i-1) 1723 jx := ix + incX 1724 xi := x[ix] 1725 for _, aij := range a[i*lda+1 : i*lda+kk+1] { 1726 x[jx] += xi * cmplx.Conj(aij) 1727 jx += incX 1728 } 1729 if diag == blas.NonUnit { 1730 x[ix] *= cmplx.Conj(a[i*lda]) 1731 } 1732 ix -= incX 1733 } 1734 } 1735 } else { 1736 if incX == 1 { 1737 for i := 0; i < n; i++ { 1738 kk := min(k, i) 1739 xi := x[i] 1740 for j, aij := range a[i*lda+k-kk : i*lda+k] { 1741 x[i-kk+j] += xi * cmplx.Conj(aij) 1742 } 1743 if diag == blas.NonUnit { 1744 x[i] *= cmplx.Conj(a[i*lda+k]) 1745 } 1746 } 1747 } else { 1748 ix := kx 1749 for i := 0; i < n; i++ { 1750 kk := min(k, i) 1751 jx := ix - kk*incX 1752 xi := x[ix] 1753 for _, aij := range a[i*lda+k-kk : i*lda+k] { 1754 x[jx] += xi * cmplx.Conj(aij) 1755 jx += incX 1756 } 1757 if diag == blas.NonUnit { 1758 x[ix] *= cmplx.Conj(a[i*lda+k]) 1759 } 1760 ix += incX 1761 } 1762 } 1763 } 1764 } 1765 } 1766 1767 // Ztbsv solves one of the systems of equations 1768 // A * x = b if trans == blas.NoTrans 1769 // Aᵀ * x = b if trans == blas.Trans 1770 // Aᴴ * x = b if trans == blas.ConjTrans 1771 // where b and x are n element vectors and A is an n×n triangular band matrix 1772 // with (k+1) diagonals. 1773 // 1774 // On entry, x contains the values of b, and the solution is 1775 // stored in-place into x. 1776 // 1777 // No test for singularity or near-singularity is included in this 1778 // routine. Such tests must be performed before calling this routine. 1779 func (Implementation) Ztbsv(uplo blas.Uplo, trans blas.Transpose, diag blas.Diag, n, k int, a []complex128, lda int, x []complex128, incX int) { 1780 switch trans { 1781 default: 1782 panic(badTranspose) 1783 case blas.NoTrans, blas.Trans, blas.ConjTrans: 1784 } 1785 switch uplo { 1786 default: 1787 panic(badUplo) 1788 case blas.Upper, blas.Lower: 1789 } 1790 switch diag { 1791 default: 1792 panic(badDiag) 1793 case blas.NonUnit, blas.Unit: 1794 } 1795 if n < 0 { 1796 panic(nLT0) 1797 } 1798 if k < 0 { 1799 panic(kLT0) 1800 } 1801 if lda < k+1 { 1802 panic(badLdA) 1803 } 1804 if incX == 0 { 1805 panic(zeroIncX) 1806 } 1807 1808 // Quick return if possible. 1809 if n == 0 { 1810 return 1811 } 1812 1813 // For zero matrix size the following slice length checks are trivially satisfied. 1814 if len(a) < lda*(n-1)+k+1 { 1815 panic(shortA) 1816 } 1817 if (incX > 0 && len(x) <= (n-1)*incX) || (incX < 0 && len(x) <= (1-n)*incX) { 1818 panic(shortX) 1819 } 1820 1821 // Set up start index in X. 1822 var kx int 1823 if incX < 0 { 1824 kx = (1 - n) * incX 1825 } 1826 1827 switch trans { 1828 case blas.NoTrans: 1829 if uplo == blas.Upper { 1830 if incX == 1 { 1831 for i := n - 1; i >= 0; i-- { 1832 kk := min(k, n-i-1) 1833 var sum complex128 1834 for j, aij := range a[i*lda+1 : i*lda+kk+1] { 1835 sum += x[i+1+j] * aij 1836 } 1837 x[i] -= sum 1838 if diag == blas.NonUnit { 1839 x[i] /= a[i*lda] 1840 } 1841 } 1842 } else { 1843 ix := kx + (n-1)*incX 1844 for i := n - 1; i >= 0; i-- { 1845 kk := min(k, n-i-1) 1846 var sum complex128 1847 jx := ix + incX 1848 for _, aij := range a[i*lda+1 : i*lda+kk+1] { 1849 sum += x[jx] * aij 1850 jx += incX 1851 } 1852 x[ix] -= sum 1853 if diag == blas.NonUnit { 1854 x[ix] /= a[i*lda] 1855 } 1856 ix -= incX 1857 } 1858 } 1859 } else { 1860 if incX == 1 { 1861 for i := 0; i < n; i++ { 1862 kk := min(k, i) 1863 var sum complex128 1864 for j, aij := range a[i*lda+k-kk : i*lda+k] { 1865 sum += x[i-kk+j] * aij 1866 } 1867 x[i] -= sum 1868 if diag == blas.NonUnit { 1869 x[i] /= a[i*lda+k] 1870 } 1871 } 1872 } else { 1873 ix := kx 1874 for i := 0; i < n; i++ { 1875 kk := min(k, i) 1876 var sum complex128 1877 jx := ix - kk*incX 1878 for _, aij := range a[i*lda+k-kk : i*lda+k] { 1879 sum += x[jx] * aij 1880 jx += incX 1881 } 1882 x[ix] -= sum 1883 if diag == blas.NonUnit { 1884 x[ix] /= a[i*lda+k] 1885 } 1886 ix += incX 1887 } 1888 } 1889 } 1890 case blas.Trans: 1891 if uplo == blas.Upper { 1892 if incX == 1 { 1893 for i := 0; i < n; i++ { 1894 if diag == blas.NonUnit { 1895 x[i] /= a[i*lda] 1896 } 1897 kk := min(k, n-i-1) 1898 xi := x[i] 1899 for j, aij := range a[i*lda+1 : i*lda+kk+1] { 1900 x[i+1+j] -= xi * aij 1901 } 1902 } 1903 } else { 1904 ix := kx 1905 for i := 0; i < n; i++ { 1906 if diag == blas.NonUnit { 1907 x[ix] /= a[i*lda] 1908 } 1909 kk := min(k, n-i-1) 1910 xi := x[ix] 1911 jx := ix + incX 1912 for _, aij := range a[i*lda+1 : i*lda+kk+1] { 1913 x[jx] -= xi * aij 1914 jx += incX 1915 } 1916 ix += incX 1917 } 1918 } 1919 } else { 1920 if incX == 1 { 1921 for i := n - 1; i >= 0; i-- { 1922 if diag == blas.NonUnit { 1923 x[i] /= a[i*lda+k] 1924 } 1925 kk := min(k, i) 1926 xi := x[i] 1927 for j, aij := range a[i*lda+k-kk : i*lda+k] { 1928 x[i-kk+j] -= xi * aij 1929 } 1930 } 1931 } else { 1932 ix := kx + (n-1)*incX 1933 for i := n - 1; i >= 0; i-- { 1934 if diag == blas.NonUnit { 1935 x[ix] /= a[i*lda+k] 1936 } 1937 kk := min(k, i) 1938 xi := x[ix] 1939 jx := ix - kk*incX 1940 for _, aij := range a[i*lda+k-kk : i*lda+k] { 1941 x[jx] -= xi * aij 1942 jx += incX 1943 } 1944 ix -= incX 1945 } 1946 } 1947 } 1948 case blas.ConjTrans: 1949 if uplo == blas.Upper { 1950 if incX == 1 { 1951 for i := 0; i < n; i++ { 1952 if diag == blas.NonUnit { 1953 x[i] /= cmplx.Conj(a[i*lda]) 1954 } 1955 kk := min(k, n-i-1) 1956 xi := x[i] 1957 for j, aij := range a[i*lda+1 : i*lda+kk+1] { 1958 x[i+1+j] -= xi * cmplx.Conj(aij) 1959 } 1960 } 1961 } else { 1962 ix := kx 1963 for i := 0; i < n; i++ { 1964 if diag == blas.NonUnit { 1965 x[ix] /= cmplx.Conj(a[i*lda]) 1966 } 1967 kk := min(k, n-i-1) 1968 xi := x[ix] 1969 jx := ix + incX 1970 for _, aij := range a[i*lda+1 : i*lda+kk+1] { 1971 x[jx] -= xi * cmplx.Conj(aij) 1972 jx += incX 1973 } 1974 ix += incX 1975 } 1976 } 1977 } else { 1978 if incX == 1 { 1979 for i := n - 1; i >= 0; i-- { 1980 if diag == blas.NonUnit { 1981 x[i] /= cmplx.Conj(a[i*lda+k]) 1982 } 1983 kk := min(k, i) 1984 xi := x[i] 1985 for j, aij := range a[i*lda+k-kk : i*lda+k] { 1986 x[i-kk+j] -= xi * cmplx.Conj(aij) 1987 } 1988 } 1989 } else { 1990 ix := kx + (n-1)*incX 1991 for i := n - 1; i >= 0; i-- { 1992 if diag == blas.NonUnit { 1993 x[ix] /= cmplx.Conj(a[i*lda+k]) 1994 } 1995 kk := min(k, i) 1996 xi := x[ix] 1997 jx := ix - kk*incX 1998 for _, aij := range a[i*lda+k-kk : i*lda+k] { 1999 x[jx] -= xi * cmplx.Conj(aij) 2000 jx += incX 2001 } 2002 ix -= incX 2003 } 2004 } 2005 } 2006 } 2007 } 2008 2009 // Ztpmv performs one of the matrix-vector operations 2010 // x = A * x if trans = blas.NoTrans 2011 // x = Aᵀ * x if trans = blas.Trans 2012 // x = Aᴴ * x if trans = blas.ConjTrans 2013 // where x is an n element vector and A is an n×n triangular matrix, supplied in 2014 // packed form. 2015 func (Implementation) Ztpmv(uplo blas.Uplo, trans blas.Transpose, diag blas.Diag, n int, ap []complex128, x []complex128, incX int) { 2016 switch uplo { 2017 default: 2018 panic(badUplo) 2019 case blas.Upper, blas.Lower: 2020 } 2021 switch trans { 2022 default: 2023 panic(badTranspose) 2024 case blas.NoTrans, blas.Trans, blas.ConjTrans: 2025 } 2026 switch diag { 2027 default: 2028 panic(badDiag) 2029 case blas.NonUnit, blas.Unit: 2030 } 2031 if n < 0 { 2032 panic(nLT0) 2033 } 2034 if incX == 0 { 2035 panic(zeroIncX) 2036 } 2037 2038 // Quick return if possible. 2039 if n == 0 { 2040 return 2041 } 2042 2043 // For zero matrix size the following slice length checks are trivially satisfied. 2044 if len(ap) < n*(n+1)/2 { 2045 panic(shortAP) 2046 } 2047 if (incX > 0 && len(x) <= (n-1)*incX) || (incX < 0 && len(x) <= (1-n)*incX) { 2048 panic(shortX) 2049 } 2050 2051 // Set up start index in X. 2052 var kx int 2053 if incX < 0 { 2054 kx = (1 - n) * incX 2055 } 2056 2057 // The elements of A are accessed sequentially with one pass through A. 2058 2059 if trans == blas.NoTrans { 2060 // Form x = A*x. 2061 if uplo == blas.Upper { 2062 // kk points to the current diagonal element in ap. 2063 kk := 0 2064 if incX == 1 { 2065 x = x[:n] 2066 for i := range x { 2067 if diag == blas.NonUnit { 2068 x[i] *= ap[kk] 2069 } 2070 if n-i-1 > 0 { 2071 x[i] += c128.DotuUnitary(ap[kk+1:kk+n-i], x[i+1:]) 2072 } 2073 kk += n - i 2074 } 2075 } else { 2076 ix := kx 2077 for i := 0; i < n; i++ { 2078 if diag == blas.NonUnit { 2079 x[ix] *= ap[kk] 2080 } 2081 if n-i-1 > 0 { 2082 x[ix] += c128.DotuInc(ap[kk+1:kk+n-i], x, uintptr(n-i-1), 1, uintptr(incX), 0, uintptr(ix+incX)) 2083 } 2084 ix += incX 2085 kk += n - i 2086 } 2087 } 2088 } else { 2089 // kk points to the beginning of current row in ap. 2090 kk := n*(n+1)/2 - n 2091 if incX == 1 { 2092 for i := n - 1; i >= 0; i-- { 2093 if diag == blas.NonUnit { 2094 x[i] *= ap[kk+i] 2095 } 2096 if i > 0 { 2097 x[i] += c128.DotuUnitary(ap[kk:kk+i], x[:i]) 2098 } 2099 kk -= i 2100 } 2101 } else { 2102 ix := kx + (n-1)*incX 2103 for i := n - 1; i >= 0; i-- { 2104 if diag == blas.NonUnit { 2105 x[ix] *= ap[kk+i] 2106 } 2107 if i > 0 { 2108 x[ix] += c128.DotuInc(ap[kk:kk+i], x, uintptr(i), 1, uintptr(incX), 0, uintptr(kx)) 2109 } 2110 ix -= incX 2111 kk -= i 2112 } 2113 } 2114 } 2115 return 2116 } 2117 2118 if trans == blas.Trans { 2119 // Form x = Aᵀ*x. 2120 if uplo == blas.Upper { 2121 // kk points to the current diagonal element in ap. 2122 kk := n*(n+1)/2 - 1 2123 if incX == 1 { 2124 for i := n - 1; i >= 0; i-- { 2125 xi := x[i] 2126 if diag == blas.NonUnit { 2127 x[i] *= ap[kk] 2128 } 2129 if n-i-1 > 0 { 2130 c128.AxpyUnitary(xi, ap[kk+1:kk+n-i], x[i+1:n]) 2131 } 2132 kk -= n - i + 1 2133 } 2134 } else { 2135 ix := kx + (n-1)*incX 2136 for i := n - 1; i >= 0; i-- { 2137 xi := x[ix] 2138 if diag == blas.NonUnit { 2139 x[ix] *= ap[kk] 2140 } 2141 if n-i-1 > 0 { 2142 c128.AxpyInc(xi, ap[kk+1:kk+n-i], x, uintptr(n-i-1), 1, uintptr(incX), 0, uintptr(ix+incX)) 2143 } 2144 ix -= incX 2145 kk -= n - i + 1 2146 } 2147 } 2148 } else { 2149 // kk points to the beginning of current row in ap. 2150 kk := 0 2151 if incX == 1 { 2152 x = x[:n] 2153 for i := range x { 2154 if i > 0 { 2155 c128.AxpyUnitary(x[i], ap[kk:kk+i], x[:i]) 2156 } 2157 if diag == blas.NonUnit { 2158 x[i] *= ap[kk+i] 2159 } 2160 kk += i + 1 2161 } 2162 } else { 2163 ix := kx 2164 for i := 0; i < n; i++ { 2165 if i > 0 { 2166 c128.AxpyInc(x[ix], ap[kk:kk+i], x, uintptr(i), 1, uintptr(incX), 0, uintptr(kx)) 2167 } 2168 if diag == blas.NonUnit { 2169 x[ix] *= ap[kk+i] 2170 } 2171 ix += incX 2172 kk += i + 1 2173 } 2174 } 2175 } 2176 return 2177 } 2178 2179 // Form x = Aᴴ*x. 2180 if uplo == blas.Upper { 2181 // kk points to the current diagonal element in ap. 2182 kk := n*(n+1)/2 - 1 2183 if incX == 1 { 2184 for i := n - 1; i >= 0; i-- { 2185 xi := x[i] 2186 if diag == blas.NonUnit { 2187 x[i] *= cmplx.Conj(ap[kk]) 2188 } 2189 k := kk + 1 2190 for j := i + 1; j < n; j++ { 2191 x[j] += xi * cmplx.Conj(ap[k]) 2192 k++ 2193 } 2194 kk -= n - i + 1 2195 } 2196 } else { 2197 ix := kx + (n-1)*incX 2198 for i := n - 1; i >= 0; i-- { 2199 xi := x[ix] 2200 if diag == blas.NonUnit { 2201 x[ix] *= cmplx.Conj(ap[kk]) 2202 } 2203 jx := ix + incX 2204 k := kk + 1 2205 for j := i + 1; j < n; j++ { 2206 x[jx] += xi * cmplx.Conj(ap[k]) 2207 jx += incX 2208 k++ 2209 } 2210 ix -= incX 2211 kk -= n - i + 1 2212 } 2213 } 2214 } else { 2215 // kk points to the beginning of current row in ap. 2216 kk := 0 2217 if incX == 1 { 2218 x = x[:n] 2219 for i, xi := range x { 2220 for j := 0; j < i; j++ { 2221 x[j] += xi * cmplx.Conj(ap[kk+j]) 2222 } 2223 if diag == blas.NonUnit { 2224 x[i] *= cmplx.Conj(ap[kk+i]) 2225 } 2226 kk += i + 1 2227 } 2228 } else { 2229 ix := kx 2230 for i := 0; i < n; i++ { 2231 xi := x[ix] 2232 jx := kx 2233 for j := 0; j < i; j++ { 2234 x[jx] += xi * cmplx.Conj(ap[kk+j]) 2235 jx += incX 2236 } 2237 if diag == blas.NonUnit { 2238 x[ix] *= cmplx.Conj(ap[kk+i]) 2239 } 2240 ix += incX 2241 kk += i + 1 2242 } 2243 } 2244 } 2245 } 2246 2247 // Ztpsv solves one of the systems of equations 2248 // A * x = b if trans == blas.NoTrans 2249 // Aᵀ * x = b if trans == blas.Trans 2250 // Aᴴ * x = b if trans == blas.ConjTrans 2251 // where b and x are n element vectors and A is an n×n triangular matrix in 2252 // packed form. 2253 // 2254 // On entry, x contains the values of b, and the solution is 2255 // stored in-place into x. 2256 // 2257 // No test for singularity or near-singularity is included in this 2258 // routine. Such tests must be performed before calling this routine. 2259 func (Implementation) Ztpsv(uplo blas.Uplo, trans blas.Transpose, diag blas.Diag, n int, ap []complex128, x []complex128, incX int) { 2260 switch uplo { 2261 default: 2262 panic(badUplo) 2263 case blas.Upper, blas.Lower: 2264 } 2265 switch trans { 2266 default: 2267 panic(badTranspose) 2268 case blas.NoTrans, blas.Trans, blas.ConjTrans: 2269 } 2270 switch diag { 2271 default: 2272 panic(badDiag) 2273 case blas.NonUnit, blas.Unit: 2274 } 2275 if n < 0 { 2276 panic(nLT0) 2277 } 2278 if incX == 0 { 2279 panic(zeroIncX) 2280 } 2281 2282 // Quick return if possible. 2283 if n == 0 { 2284 return 2285 } 2286 2287 // For zero matrix size the following slice length checks are trivially satisfied. 2288 if len(ap) < n*(n+1)/2 { 2289 panic(shortAP) 2290 } 2291 if (incX > 0 && len(x) <= (n-1)*incX) || (incX < 0 && len(x) <= (1-n)*incX) { 2292 panic(shortX) 2293 } 2294 2295 // Set up start index in X. 2296 var kx int 2297 if incX < 0 { 2298 kx = (1 - n) * incX 2299 } 2300 2301 // The elements of A are accessed sequentially with one pass through ap. 2302 2303 if trans == blas.NoTrans { 2304 // Form x = inv(A)*x. 2305 if uplo == blas.Upper { 2306 kk := n*(n+1)/2 - 1 2307 if incX == 1 { 2308 for i := n - 1; i >= 0; i-- { 2309 aii := ap[kk] 2310 if n-i-1 > 0 { 2311 x[i] -= c128.DotuUnitary(x[i+1:n], ap[kk+1:kk+n-i]) 2312 } 2313 if diag == blas.NonUnit { 2314 x[i] /= aii 2315 } 2316 kk -= n - i + 1 2317 } 2318 } else { 2319 ix := kx + (n-1)*incX 2320 for i := n - 1; i >= 0; i-- { 2321 aii := ap[kk] 2322 if n-i-1 > 0 { 2323 x[ix] -= c128.DotuInc(x, ap[kk+1:kk+n-i], uintptr(n-i-1), uintptr(incX), 1, uintptr(ix+incX), 0) 2324 } 2325 if diag == blas.NonUnit { 2326 x[ix] /= aii 2327 } 2328 ix -= incX 2329 kk -= n - i + 1 2330 } 2331 } 2332 } else { 2333 kk := 0 2334 if incX == 1 { 2335 for i := 0; i < n; i++ { 2336 if i > 0 { 2337 x[i] -= c128.DotuUnitary(x[:i], ap[kk:kk+i]) 2338 } 2339 if diag == blas.NonUnit { 2340 x[i] /= ap[kk+i] 2341 } 2342 kk += i + 1 2343 } 2344 } else { 2345 ix := kx 2346 for i := 0; i < n; i++ { 2347 if i > 0 { 2348 x[ix] -= c128.DotuInc(x, ap[kk:kk+i], uintptr(i), uintptr(incX), 1, uintptr(kx), 0) 2349 } 2350 if diag == blas.NonUnit { 2351 x[ix] /= ap[kk+i] 2352 } 2353 ix += incX 2354 kk += i + 1 2355 } 2356 } 2357 } 2358 return 2359 } 2360 2361 if trans == blas.Trans { 2362 // Form x = inv(Aᵀ)*x. 2363 if uplo == blas.Upper { 2364 kk := 0 2365 if incX == 1 { 2366 for j := 0; j < n; j++ { 2367 if diag == blas.NonUnit { 2368 x[j] /= ap[kk] 2369 } 2370 if n-j-1 > 0 { 2371 c128.AxpyUnitary(-x[j], ap[kk+1:kk+n-j], x[j+1:n]) 2372 } 2373 kk += n - j 2374 } 2375 } else { 2376 jx := kx 2377 for j := 0; j < n; j++ { 2378 if diag == blas.NonUnit { 2379 x[jx] /= ap[kk] 2380 } 2381 if n-j-1 > 0 { 2382 c128.AxpyInc(-x[jx], ap[kk+1:kk+n-j], x, uintptr(n-j-1), 1, uintptr(incX), 0, uintptr(jx+incX)) 2383 } 2384 jx += incX 2385 kk += n - j 2386 } 2387 } 2388 } else { 2389 kk := n*(n+1)/2 - n 2390 if incX == 1 { 2391 for j := n - 1; j >= 0; j-- { 2392 if diag == blas.NonUnit { 2393 x[j] /= ap[kk+j] 2394 } 2395 if j > 0 { 2396 c128.AxpyUnitary(-x[j], ap[kk:kk+j], x[:j]) 2397 } 2398 kk -= j 2399 } 2400 } else { 2401 jx := kx + (n-1)*incX 2402 for j := n - 1; j >= 0; j-- { 2403 if diag == blas.NonUnit { 2404 x[jx] /= ap[kk+j] 2405 } 2406 if j > 0 { 2407 c128.AxpyInc(-x[jx], ap[kk:kk+j], x, uintptr(j), 1, uintptr(incX), 0, uintptr(kx)) 2408 } 2409 jx -= incX 2410 kk -= j 2411 } 2412 } 2413 } 2414 return 2415 } 2416 2417 // Form x = inv(Aᴴ)*x. 2418 if uplo == blas.Upper { 2419 kk := 0 2420 if incX == 1 { 2421 for j := 0; j < n; j++ { 2422 if diag == blas.NonUnit { 2423 x[j] /= cmplx.Conj(ap[kk]) 2424 } 2425 xj := x[j] 2426 k := kk + 1 2427 for i := j + 1; i < n; i++ { 2428 x[i] -= xj * cmplx.Conj(ap[k]) 2429 k++ 2430 } 2431 kk += n - j 2432 } 2433 } else { 2434 jx := kx 2435 for j := 0; j < n; j++ { 2436 if diag == blas.NonUnit { 2437 x[jx] /= cmplx.Conj(ap[kk]) 2438 } 2439 xj := x[jx] 2440 ix := jx + incX 2441 k := kk + 1 2442 for i := j + 1; i < n; i++ { 2443 x[ix] -= xj * cmplx.Conj(ap[k]) 2444 ix += incX 2445 k++ 2446 } 2447 jx += incX 2448 kk += n - j 2449 } 2450 } 2451 } else { 2452 kk := n*(n+1)/2 - n 2453 if incX == 1 { 2454 for j := n - 1; j >= 0; j-- { 2455 if diag == blas.NonUnit { 2456 x[j] /= cmplx.Conj(ap[kk+j]) 2457 } 2458 xj := x[j] 2459 for i := 0; i < j; i++ { 2460 x[i] -= xj * cmplx.Conj(ap[kk+i]) 2461 } 2462 kk -= j 2463 } 2464 } else { 2465 jx := kx + (n-1)*incX 2466 for j := n - 1; j >= 0; j-- { 2467 if diag == blas.NonUnit { 2468 x[jx] /= cmplx.Conj(ap[kk+j]) 2469 } 2470 xj := x[jx] 2471 ix := kx 2472 for i := 0; i < j; i++ { 2473 x[ix] -= xj * cmplx.Conj(ap[kk+i]) 2474 ix += incX 2475 } 2476 jx -= incX 2477 kk -= j 2478 } 2479 } 2480 } 2481 } 2482 2483 // Ztrmv performs one of the matrix-vector operations 2484 // x = A * x if trans = blas.NoTrans 2485 // x = Aᵀ * x if trans = blas.Trans 2486 // x = Aᴴ * x if trans = blas.ConjTrans 2487 // where x is a vector, and A is an n×n triangular matrix. 2488 func (Implementation) Ztrmv(uplo blas.Uplo, trans blas.Transpose, diag blas.Diag, n int, a []complex128, lda int, x []complex128, incX int) { 2489 switch trans { 2490 default: 2491 panic(badTranspose) 2492 case blas.NoTrans, blas.Trans, blas.ConjTrans: 2493 } 2494 switch uplo { 2495 default: 2496 panic(badUplo) 2497 case blas.Upper, blas.Lower: 2498 } 2499 switch diag { 2500 default: 2501 panic(badDiag) 2502 case blas.NonUnit, blas.Unit: 2503 } 2504 if n < 0 { 2505 panic(nLT0) 2506 } 2507 if lda < max(1, n) { 2508 panic(badLdA) 2509 } 2510 if incX == 0 { 2511 panic(zeroIncX) 2512 } 2513 2514 // Quick return if possible. 2515 if n == 0 { 2516 return 2517 } 2518 2519 // For zero matrix size the following slice length checks are trivially satisfied. 2520 if len(a) < lda*(n-1)+n { 2521 panic(shortA) 2522 } 2523 if (incX > 0 && len(x) <= (n-1)*incX) || (incX < 0 && len(x) <= (1-n)*incX) { 2524 panic(shortX) 2525 } 2526 2527 // Set up start index in X. 2528 var kx int 2529 if incX < 0 { 2530 kx = (1 - n) * incX 2531 } 2532 2533 // The elements of A are accessed sequentially with one pass through A. 2534 2535 if trans == blas.NoTrans { 2536 // Form x = A*x. 2537 if uplo == blas.Upper { 2538 if incX == 1 { 2539 for i := 0; i < n; i++ { 2540 if diag == blas.NonUnit { 2541 x[i] *= a[i*lda+i] 2542 } 2543 if n-i-1 > 0 { 2544 x[i] += c128.DotuUnitary(a[i*lda+i+1:i*lda+n], x[i+1:n]) 2545 } 2546 } 2547 } else { 2548 ix := kx 2549 for i := 0; i < n; i++ { 2550 if diag == blas.NonUnit { 2551 x[ix] *= a[i*lda+i] 2552 } 2553 if n-i-1 > 0 { 2554 x[ix] += c128.DotuInc(a[i*lda+i+1:i*lda+n], x, uintptr(n-i-1), 1, uintptr(incX), 0, uintptr(ix+incX)) 2555 } 2556 ix += incX 2557 } 2558 } 2559 } else { 2560 if incX == 1 { 2561 for i := n - 1; i >= 0; i-- { 2562 if diag == blas.NonUnit { 2563 x[i] *= a[i*lda+i] 2564 } 2565 if i > 0 { 2566 x[i] += c128.DotuUnitary(a[i*lda:i*lda+i], x[:i]) 2567 } 2568 } 2569 } else { 2570 ix := kx + (n-1)*incX 2571 for i := n - 1; i >= 0; i-- { 2572 if diag == blas.NonUnit { 2573 x[ix] *= a[i*lda+i] 2574 } 2575 if i > 0 { 2576 x[ix] += c128.DotuInc(a[i*lda:i*lda+i], x, uintptr(i), 1, uintptr(incX), 0, uintptr(kx)) 2577 } 2578 ix -= incX 2579 } 2580 } 2581 } 2582 return 2583 } 2584 2585 if trans == blas.Trans { 2586 // Form x = Aᵀ*x. 2587 if uplo == blas.Upper { 2588 if incX == 1 { 2589 for i := n - 1; i >= 0; i-- { 2590 xi := x[i] 2591 if diag == blas.NonUnit { 2592 x[i] *= a[i*lda+i] 2593 } 2594 if n-i-1 > 0 { 2595 c128.AxpyUnitary(xi, a[i*lda+i+1:i*lda+n], x[i+1:n]) 2596 } 2597 } 2598 } else { 2599 ix := kx + (n-1)*incX 2600 for i := n - 1; i >= 0; i-- { 2601 xi := x[ix] 2602 if diag == blas.NonUnit { 2603 x[ix] *= a[i*lda+i] 2604 } 2605 if n-i-1 > 0 { 2606 c128.AxpyInc(xi, a[i*lda+i+1:i*lda+n], x, uintptr(n-i-1), 1, uintptr(incX), 0, uintptr(ix+incX)) 2607 } 2608 ix -= incX 2609 } 2610 } 2611 } else { 2612 if incX == 1 { 2613 for i := 0; i < n; i++ { 2614 if i > 0 { 2615 c128.AxpyUnitary(x[i], a[i*lda:i*lda+i], x[:i]) 2616 } 2617 if diag == blas.NonUnit { 2618 x[i] *= a[i*lda+i] 2619 } 2620 } 2621 } else { 2622 ix := kx 2623 for i := 0; i < n; i++ { 2624 if i > 0 { 2625 c128.AxpyInc(x[ix], a[i*lda:i*lda+i], x, uintptr(i), 1, uintptr(incX), 0, uintptr(kx)) 2626 } 2627 if diag == blas.NonUnit { 2628 x[ix] *= a[i*lda+i] 2629 } 2630 ix += incX 2631 } 2632 } 2633 } 2634 return 2635 } 2636 2637 // Form x = Aᴴ*x. 2638 if uplo == blas.Upper { 2639 if incX == 1 { 2640 for i := n - 1; i >= 0; i-- { 2641 xi := x[i] 2642 if diag == blas.NonUnit { 2643 x[i] *= cmplx.Conj(a[i*lda+i]) 2644 } 2645 for j := i + 1; j < n; j++ { 2646 x[j] += xi * cmplx.Conj(a[i*lda+j]) 2647 } 2648 } 2649 } else { 2650 ix := kx + (n-1)*incX 2651 for i := n - 1; i >= 0; i-- { 2652 xi := x[ix] 2653 if diag == blas.NonUnit { 2654 x[ix] *= cmplx.Conj(a[i*lda+i]) 2655 } 2656 jx := ix + incX 2657 for j := i + 1; j < n; j++ { 2658 x[jx] += xi * cmplx.Conj(a[i*lda+j]) 2659 jx += incX 2660 } 2661 ix -= incX 2662 } 2663 } 2664 } else { 2665 if incX == 1 { 2666 for i := 0; i < n; i++ { 2667 for j := 0; j < i; j++ { 2668 x[j] += x[i] * cmplx.Conj(a[i*lda+j]) 2669 } 2670 if diag == blas.NonUnit { 2671 x[i] *= cmplx.Conj(a[i*lda+i]) 2672 } 2673 } 2674 } else { 2675 ix := kx 2676 for i := 0; i < n; i++ { 2677 jx := kx 2678 for j := 0; j < i; j++ { 2679 x[jx] += x[ix] * cmplx.Conj(a[i*lda+j]) 2680 jx += incX 2681 } 2682 if diag == blas.NonUnit { 2683 x[ix] *= cmplx.Conj(a[i*lda+i]) 2684 } 2685 ix += incX 2686 } 2687 } 2688 } 2689 } 2690 2691 // Ztrsv solves one of the systems of equations 2692 // A * x = b if trans == blas.NoTrans 2693 // Aᵀ * x = b if trans == blas.Trans 2694 // Aᴴ * x = b if trans == blas.ConjTrans 2695 // where b and x are n element vectors and A is an n×n triangular matrix. 2696 // 2697 // On entry, x contains the values of b, and the solution is 2698 // stored in-place into x. 2699 // 2700 // No test for singularity or near-singularity is included in this 2701 // routine. Such tests must be performed before calling this routine. 2702 func (Implementation) Ztrsv(uplo blas.Uplo, trans blas.Transpose, diag blas.Diag, n int, a []complex128, lda int, x []complex128, incX int) { 2703 switch trans { 2704 default: 2705 panic(badTranspose) 2706 case blas.NoTrans, blas.Trans, blas.ConjTrans: 2707 } 2708 switch uplo { 2709 default: 2710 panic(badUplo) 2711 case blas.Upper, blas.Lower: 2712 } 2713 switch diag { 2714 default: 2715 panic(badDiag) 2716 case blas.NonUnit, blas.Unit: 2717 } 2718 if n < 0 { 2719 panic(nLT0) 2720 } 2721 if lda < max(1, n) { 2722 panic(badLdA) 2723 } 2724 if incX == 0 { 2725 panic(zeroIncX) 2726 } 2727 2728 // Quick return if possible. 2729 if n == 0 { 2730 return 2731 } 2732 2733 // For zero matrix size the following slice length checks are trivially satisfied. 2734 if len(a) < lda*(n-1)+n { 2735 panic(shortA) 2736 } 2737 if (incX > 0 && len(x) <= (n-1)*incX) || (incX < 0 && len(x) <= (1-n)*incX) { 2738 panic(shortX) 2739 } 2740 2741 // Set up start index in X. 2742 var kx int 2743 if incX < 0 { 2744 kx = (1 - n) * incX 2745 } 2746 2747 // The elements of A are accessed sequentially with one pass through A. 2748 2749 if trans == blas.NoTrans { 2750 // Form x = inv(A)*x. 2751 if uplo == blas.Upper { 2752 if incX == 1 { 2753 for i := n - 1; i >= 0; i-- { 2754 aii := a[i*lda+i] 2755 if n-i-1 > 0 { 2756 x[i] -= c128.DotuUnitary(x[i+1:n], a[i*lda+i+1:i*lda+n]) 2757 } 2758 if diag == blas.NonUnit { 2759 x[i] /= aii 2760 } 2761 } 2762 } else { 2763 ix := kx + (n-1)*incX 2764 for i := n - 1; i >= 0; i-- { 2765 aii := a[i*lda+i] 2766 if n-i-1 > 0 { 2767 x[ix] -= c128.DotuInc(x, a[i*lda+i+1:i*lda+n], uintptr(n-i-1), uintptr(incX), 1, uintptr(ix+incX), 0) 2768 } 2769 if diag == blas.NonUnit { 2770 x[ix] /= aii 2771 } 2772 ix -= incX 2773 } 2774 } 2775 } else { 2776 if incX == 1 { 2777 for i := 0; i < n; i++ { 2778 if i > 0 { 2779 x[i] -= c128.DotuUnitary(x[:i], a[i*lda:i*lda+i]) 2780 } 2781 if diag == blas.NonUnit { 2782 x[i] /= a[i*lda+i] 2783 } 2784 } 2785 } else { 2786 ix := kx 2787 for i := 0; i < n; i++ { 2788 if i > 0 { 2789 x[ix] -= c128.DotuInc(x, a[i*lda:i*lda+i], uintptr(i), uintptr(incX), 1, uintptr(kx), 0) 2790 } 2791 if diag == blas.NonUnit { 2792 x[ix] /= a[i*lda+i] 2793 } 2794 ix += incX 2795 } 2796 } 2797 } 2798 return 2799 } 2800 2801 if trans == blas.Trans { 2802 // Form x = inv(Aᵀ)*x. 2803 if uplo == blas.Upper { 2804 if incX == 1 { 2805 for j := 0; j < n; j++ { 2806 if diag == blas.NonUnit { 2807 x[j] /= a[j*lda+j] 2808 } 2809 if n-j-1 > 0 { 2810 c128.AxpyUnitary(-x[j], a[j*lda+j+1:j*lda+n], x[j+1:n]) 2811 } 2812 } 2813 } else { 2814 jx := kx 2815 for j := 0; j < n; j++ { 2816 if diag == blas.NonUnit { 2817 x[jx] /= a[j*lda+j] 2818 } 2819 if n-j-1 > 0 { 2820 c128.AxpyInc(-x[jx], a[j*lda+j+1:j*lda+n], x, uintptr(n-j-1), 1, uintptr(incX), 0, uintptr(jx+incX)) 2821 } 2822 jx += incX 2823 } 2824 } 2825 } else { 2826 if incX == 1 { 2827 for j := n - 1; j >= 0; j-- { 2828 if diag == blas.NonUnit { 2829 x[j] /= a[j*lda+j] 2830 } 2831 xj := x[j] 2832 if j > 0 { 2833 c128.AxpyUnitary(-xj, a[j*lda:j*lda+j], x[:j]) 2834 } 2835 } 2836 } else { 2837 jx := kx + (n-1)*incX 2838 for j := n - 1; j >= 0; j-- { 2839 if diag == blas.NonUnit { 2840 x[jx] /= a[j*lda+j] 2841 } 2842 if j > 0 { 2843 c128.AxpyInc(-x[jx], a[j*lda:j*lda+j], x, uintptr(j), 1, uintptr(incX), 0, uintptr(kx)) 2844 } 2845 jx -= incX 2846 } 2847 } 2848 } 2849 return 2850 } 2851 2852 // Form x = inv(Aᴴ)*x. 2853 if uplo == blas.Upper { 2854 if incX == 1 { 2855 for j := 0; j < n; j++ { 2856 if diag == blas.NonUnit { 2857 x[j] /= cmplx.Conj(a[j*lda+j]) 2858 } 2859 xj := x[j] 2860 for i := j + 1; i < n; i++ { 2861 x[i] -= xj * cmplx.Conj(a[j*lda+i]) 2862 } 2863 } 2864 } else { 2865 jx := kx 2866 for j := 0; j < n; j++ { 2867 if diag == blas.NonUnit { 2868 x[jx] /= cmplx.Conj(a[j*lda+j]) 2869 } 2870 xj := x[jx] 2871 ix := jx + incX 2872 for i := j + 1; i < n; i++ { 2873 x[ix] -= xj * cmplx.Conj(a[j*lda+i]) 2874 ix += incX 2875 } 2876 jx += incX 2877 } 2878 } 2879 } else { 2880 if incX == 1 { 2881 for j := n - 1; j >= 0; j-- { 2882 if diag == blas.NonUnit { 2883 x[j] /= cmplx.Conj(a[j*lda+j]) 2884 } 2885 xj := x[j] 2886 for i := 0; i < j; i++ { 2887 x[i] -= xj * cmplx.Conj(a[j*lda+i]) 2888 } 2889 } 2890 } else { 2891 jx := kx + (n-1)*incX 2892 for j := n - 1; j >= 0; j-- { 2893 if diag == blas.NonUnit { 2894 x[jx] /= cmplx.Conj(a[j*lda+j]) 2895 } 2896 xj := x[jx] 2897 ix := kx 2898 for i := 0; i < j; i++ { 2899 x[ix] -= xj * cmplx.Conj(a[j*lda+i]) 2900 ix += incX 2901 } 2902 jx -= incX 2903 } 2904 } 2905 } 2906 }