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