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

   1 // Code generated by "go generate gonum.org/v1/gonum/blas/gonum”; DO NOT EDIT.
   2
   3 // Copyright ©2014 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         "gonum.org/v1/gonum/blas"
  11         "gonum.org/v1/gonum/internal/asm/f32"
  12 )
  13
  14 var _ blas.Float32Level2 = Implementation{}
  15
  16 // Sgemv computes
  17 //  y = alpha * A * x + beta * y    if tA = blas.NoTrans
  18 //  y = alpha * A^T * x + beta * y  if tA = blas.Trans or blas.ConjTrans
  19 // where A is an m×n dense matrix, x and y are vectors, and alpha and beta are scalars.
  20 //
  21 // Float32 implementations are autogenerated and not directly tested.
  22 func (Implementation) Sgemv(tA blas.Transpose, m, n int, alpha float32, a []float32, lda int, x []float32, incX int, beta float32, y []float32, incY int) {
  23         if tA != blas.NoTrans && tA != blas.Trans && tA != blas.ConjTrans {
  24                 panic(badTranspose)
  25         }
  26         if m < 0 {
  27                 panic(mLT0)
  28         }
  29         if n < 0 {
  30                 panic(nLT0)
  31         }
  32         if lda < max(1, n) {
  33                 panic(badLdA)
  34         }
  35
  36         if incX == 0 {
  37                 panic(zeroIncX)
  38         }
  39         if incY == 0 {
  40                 panic(zeroIncY)
  41         }
  42         // Set up indexes
  43         lenX := m
  44         lenY := n
  45         if tA == blas.NoTrans {
  46                 lenX = n
  47                 lenY = m
  48         }
  49         if (incX > 0 && (lenX-1)*incX >= len(x)) || (incX < 0 && (1-lenX)*incX >= len(x)) {
  50                 panic(badX)
  51         }
  52         if (incY > 0 && (lenY-1)*incY >= len(y)) || (incY < 0 && (1-lenY)*incY >= len(y)) {
  53                 panic(badY)
  54         }
  55         if lda*(m-1)+n > len(a) || lda < max(1, n) {
  56                 panic(badLdA)
  57         }
  58
  59         // Quick return if possible
  60         if m == 0 || n == 0 || (alpha == 0 && beta == 1) {
  61                 return
  62         }
  63
  64         var kx, ky int
  65         if incX > 0 {
  66                 kx = 0
  67         } else {
  68                 kx = -(lenX - 1) * incX
  69         }
  70         if incY > 0 {
  71                 ky = 0
  72         } else {
  73                 ky = -(lenY - 1) * incY
  74         }
  75
  76         // First form y = beta * y
  77         if incY > 0 {
  78                 Implementation{}.Sscal(lenY, beta, y, incY)
  79         } else {
  80                 Implementation{}.Sscal(lenY, beta, y, -incY)
  81         }
  82
  83         if alpha == 0 {
  84                 return
  85         }
  86
  87         // Form y = alpha * A * x + y
  88         if tA == blas.NoTrans {
  89                 if incX == 1 && incY == 1 {
  90                         for i := 0; i < m; i++ {
  91                                 y[i] += alpha * f32.DotUnitary(a[lda*i:lda*i+n], x)
  92                         }
  93                         return
  94                 }
  95                 iy := ky
  96                 for i := 0; i < m; i++ {
  97                         y[iy] += alpha * f32.DotInc(x, a[lda*i:lda*i+n], uintptr(n), uintptr(incX), 1, uintptr(kx), 0)
  98                         iy += incY
  99                 }
 100                 return
 101         }
 102         // Cases where a is transposed.
 103         if incX == 1 && incY == 1 {
 104                 for i := 0; i < m; i++ {
 105                         tmp := alpha * x[i]
 106                         if tmp != 0 {
 107                                 f32.AxpyUnitaryTo(y, tmp, a[lda*i:lda*i+n], y)
 108                         }
 109                 }
 110                 return
 111         }
 112         ix := kx
 113         for i := 0; i < m; i++ {
 114                 tmp := alpha * x[ix]
 115                 if tmp != 0 {
 116                         f32.AxpyInc(tmp, a[lda*i:lda*i+n], y, uintptr(n), 1, uintptr(incY), 0, uintptr(ky))
 117                 }
 118                 ix += incX
 119         }
 120 }
 121
 122 // Sger performs the rank-one operation
 123 //  A += alpha * x * y^T
 124 // where A is an m×n dense matrix, x and y are vectors, and alpha is a scalar.
 125 //
 126 // Float32 implementations are autogenerated and not directly tested.
 127 func (Implementation) Sger(m, n int, alpha float32, x []float32, incX int, y []float32, incY int, a []float32, lda int) {
 128         // Check inputs
 129         if m < 0 {
 130                 panic("m < 0")
 131         }
 132         if n < 0 {
 133                 panic(negativeN)
 134         }
 135         if incX == 0 {
 136                 panic(zeroIncX)
 137         }
 138         if incY == 0 {
 139                 panic(zeroIncY)
 140         }
 141         if (incX > 0 && (m-1)*incX >= len(x)) || (incX < 0 && (1-m)*incX >= len(x)) {
 142                 panic(badX)
 143         }
 144         if (incY > 0 && (n-1)*incY >= len(y)) || (incY < 0 && (1-n)*incY >= len(y)) {
 145                 panic(badY)
 146         }
 147         if lda*(m-1)+n > len(a) || lda < max(1, n) {
 148                 panic(badLdA)
 149         }
 150         if lda < max(1, n) {
 151                 panic(badLdA)
 152         }
 153
 154         // Quick return if possible
 155         if m == 0 || n == 0 || alpha == 0 {
 156                 return
 157         }
 158
 159         var ky, kx int
 160         if incY > 0 {
 161                 ky = 0
 162         } else {
 163                 ky = -(n - 1) * incY
 164         }
 165
 166         if incX > 0 {
 167                 kx = 0
 168         } else {
 169                 kx = -(m - 1) * incX
 170         }
 171
 172         if incX == 1 && incY == 1 {
 173                 x = x[:m]
 174                 y = y[:n]
 175                 for i, xv := range x {
 176                         f32.AxpyUnitary(alpha*xv, y, a[i*lda:i*lda+n])
 177                 }
 178                 return
 179         }
 180
 181         ix := kx
 182         for i := 0; i < m; i++ {
 183                 f32.AxpyInc(alpha*x[ix], y, a[i*lda:i*lda+n], uintptr(n), uintptr(incY), 1, uintptr(ky), 0)
 184                 ix += incX
 185         }
 186 }
 187
 188 // Sgbmv computes
 189 //  y = alpha * A * x + beta * y if tA == blas.NoTrans
 190 //  y = alpha * A^T * x + beta * y if tA == blas.Trans or blas.ConjTrans
 191 // where a is an m×n band matrix kL subdiagonals and kU super-diagonals, and
 192 // m and n refer to the size of the full dense matrix it represents.
 193 // x and y are vectors, and alpha and beta are scalars.
 194 //
 195 // Float32 implementations are autogenerated and not directly tested.
 196 func (Implementation) Sgbmv(tA blas.Transpose, m, n, kL, kU int, alpha float32, a []float32, lda int, x []float32, incX int, beta float32, y []float32, incY int) {
 197         if tA != blas.NoTrans && tA != blas.Trans && tA != blas.ConjTrans {
 198                 panic(badTranspose)
 199         }
 200         if m < 0 {
 201                 panic(mLT0)
 202         }
 203         if n < 0 {
 204                 panic(nLT0)
 205         }
 206         if kL < 0 {
 207                 panic(kLLT0)
 208         }
 209         if kL < 0 {
 210                 panic(kULT0)
 211         }
 212         if lda < kL+kU+1 {
 213                 panic(badLdA)
 214         }
 215         if incX == 0 {
 216                 panic(zeroIncX)
 217         }
 218         if incY == 0 {
 219                 panic(zeroIncY)
 220         }
 221         // Set up indexes
 222         lenX := m
 223         lenY := n
 224         if tA == blas.NoTrans {
 225                 lenX = n
 226                 lenY = m
 227         }
 228         if (incX > 0 && (lenX-1)*incX >= len(x)) || (incX < 0 && (1-lenX)*incX >= len(x)) {
 229                 panic(badX)
 230         }
 231         if (incY > 0 && (lenY-1)*incY >= len(y)) || (incY < 0 && (1-lenY)*incY >= len(y)) {
 232                 panic(badY)
 233         }
 234         if lda*(min(m, n+kL)-1)+kL+kU+1 > len(a) || lda < kL+kU+1 {
 235                 panic(badLdA)
 236         }
 237
 238         // Quick return if possible
 239         if m == 0 || n == 0 || (alpha == 0 && beta == 1) {
 240                 return
 241         }
 242
 243         var kx, ky int
 244         if incX > 0 {
 245                 kx = 0
 246         } else {
 247                 kx = -(lenX - 1) * incX
 248         }
 249         if incY > 0 {
 250                 ky = 0
 251         } else {
 252                 ky = -(lenY - 1) * incY
 253         }
 254
 255         // First form y = beta * y
 256         if incY > 0 {
 257                 Implementation{}.Sscal(lenY, beta, y, incY)
 258         } else {
 259                 Implementation{}.Sscal(lenY, beta, y, -incY)
 260         }
 261
 262         if alpha == 0 {
 263                 return
 264         }
 265
 266         // i and j are indices of the compacted banded matrix.
 267         // off is the offset into the dense matrix (off + j = densej)
 268         ld := min(m, n)
 269         nCol := kU + 1 + kL
 270         if tA == blas.NoTrans {
 271                 iy := ky
 272                 if incX == 1 {
 273                         for i := 0; i < min(m, n+kL); i++ {
 274                                 l := max(0, kL-i)
 275                                 u := min(nCol, ld+kL-i)
 276                                 off := max(0, i-kL)
 277                                 atmp := a[i*lda+l : i*lda+u]
 278                                 xtmp := x[off : off+u-l]
 279                                 var sum float32
 280                                 for j, v := range atmp {
 281                                         sum += xtmp[j] * v
 282                                 }
 283                                 y[iy] += sum * alpha
 284                                 iy += incY
 285                         }
 286                         return
 287                 }
 288                 for i := 0; i < min(m, n+kL); i++ {
 289                         l := max(0, kL-i)
 290                         u := min(nCol, ld+kL-i)
 291                         off := max(0, i-kL)
 292                         atmp := a[i*lda+l : i*lda+u]
 293                         jx := kx
 294                         var sum float32
 295                         for _, v := range atmp {
 296                                 sum += x[off*incX+jx] * v
 297                                 jx += incX
 298                         }
 299                         y[iy] += sum * alpha
 300                         iy += incY
 301                 }
 302                 return
 303         }
 304         if incX == 1 {
 305                 for i := 0; i < min(m, n+kL); i++ {
 306                         l := max(0, kL-i)
 307                         u := min(nCol, ld+kL-i)
 308                         off := max(0, i-kL)
 309                         atmp := a[i*lda+l : i*lda+u]
 310                         tmp := alpha * x[i]
 311                         jy := ky
 312                         for _, v := range atmp {
 313                                 y[jy+off*incY] += tmp * v
 314                                 jy += incY
 315                         }
 316                 }
 317                 return
 318         }
 319         ix := kx
 320         for i := 0; i < min(m, n+kL); i++ {
 321                 l := max(0, kL-i)
 322                 u := min(nCol, ld+kL-i)
 323                 off := max(0, i-kL)
 324                 atmp := a[i*lda+l : i*lda+u]
 325                 tmp := alpha * x[ix]
 326                 jy := ky
 327                 for _, v := range atmp {
 328                         y[jy+off*incY] += tmp * v
 329                         jy += incY
 330                 }
 331                 ix += incX
 332         }
 333 }
 334
 335 // Strmv computes
 336 //  x = A * x if tA == blas.NoTrans
 337 //  x = A^T * x if tA == blas.Trans or blas.ConjTrans
 338 // A is an n×n Triangular matrix and x is a vector.
 339 //
 340 // Float32 implementations are autogenerated and not directly tested.
 341 func (Implementation) Strmv(ul blas.Uplo, tA blas.Transpose, d blas.Diag, n int, a []float32, lda int, x []float32, incX int) {
 342         if ul != blas.Lower && ul != blas.Upper {
 343                 panic(badUplo)
 344         }
 345         if tA != blas.NoTrans && tA != blas.Trans && tA != blas.ConjTrans {
 346                 panic(badTranspose)
 347         }
 348         if d != blas.NonUnit && d != blas.Unit {
 349                 panic(badDiag)
 350         }
 351         if n < 0 {
 352                 panic(nLT0)
 353         }
 354         if lda < n {
 355                 panic(badLdA)
 356         }
 357         if incX == 0 {
 358                 panic(zeroIncX)
 359         }
 360         if (incX > 0 && (n-1)*incX >= len(x)) || (incX < 0 && (1-n)*incX >= len(x)) {
 361                 panic(badX)
 362         }
 363         if lda*(n-1)+n > len(a) || lda < max(1, n) {
 364                 panic(badLdA)
 365         }
 366         if n == 0 {
 367                 return
 368         }
 369         nonUnit := d != blas.Unit
 370         if n == 1 {
 371                 if nonUnit {
 372                         x[0] *= a[0]
 373                 }
 374                 return
 375         }
 376         var kx int
 377         if incX <= 0 {
 378                 kx = -(n - 1) * incX
 379         }
 380         if tA == blas.NoTrans {
 381                 if ul == blas.Upper {
 382                         if incX == 1 {
 383                                 for i := 0; i < n; i++ {
 384                                         ilda := i * lda
 385                                         var tmp float32
 386                                         if nonUnit {
 387                                                 tmp = a[ilda+i] * x[i]
 388                                         } else {
 389                                                 tmp = x[i]
 390                                         }
 391                                         xtmp := x[i+1:]
 392                                         x[i] = tmp + f32.DotUnitary(a[ilda+i+1:ilda+n], xtmp)
 393                                 }
 394                                 return
 395                         }
 396                         ix := kx
 397                         for i := 0; i < n; i++ {
 398                                 ilda := i * lda
 399                                 var tmp float32
 400                                 if nonUnit {
 401                                         tmp = a[ilda+i] * x[ix]
 402                                 } else {
 403                                         tmp = x[ix]
 404                                 }
 405                                 x[ix] = tmp + f32.DotInc(x, a[ilda+i+1:ilda+n], uintptr(n-i-1), uintptr(incX), 1, uintptr(ix+incX), 0)
 406                                 ix += incX
 407                         }
 408                         return
 409                 }
 410                 if incX == 1 {
 411                         for i := n - 1; i >= 0; i-- {
 412                                 ilda := i * lda
 413                                 var tmp float32
 414                                 if nonUnit {
 415                                         tmp += a[ilda+i] * x[i]
 416                                 } else {
 417                                         tmp = x[i]
 418                                 }
 419                                 x[i] = tmp + f32.DotUnitary(a[ilda:ilda+i], x)
 420                         }
 421                         return
 422                 }
 423                 ix := kx + (n-1)*incX
 424                 for i := n - 1; i >= 0; i-- {
 425                         ilda := i * lda
 426                         var tmp float32
 427                         if nonUnit {
 428                                 tmp = a[ilda+i] * x[ix]
 429                         } else {
 430                                 tmp = x[ix]
 431                         }
 432                         x[ix] = tmp + f32.DotInc(x, a[ilda:ilda+i], uintptr(i), uintptr(incX), 1, uintptr(kx), 0)
 433                         ix -= incX
 434                 }
 435                 return
 436         }
 437         // Cases where a is transposed.
 438         if ul == blas.Upper {
 439                 if incX == 1 {
 440                         for i := n - 1; i >= 0; i-- {
 441                                 ilda := i * lda
 442                                 xi := x[i]
 443                                 f32.AxpyUnitary(xi, a[ilda+i+1:ilda+n], x[i+1:n])
 444                                 if nonUnit {
 445                                         x[i] *= a[ilda+i]
 446                                 }
 447                         }
 448                         return
 449                 }
 450                 ix := kx + (n-1)*incX
 451                 for i := n - 1; i >= 0; i-- {
 452                         ilda := i * lda
 453                         xi := x[ix]
 454                         f32.AxpyInc(xi, a[ilda+i+1:ilda+n], x, uintptr(n-i-1), 1, uintptr(incX), 0, uintptr(kx+(i+1)*incX))
 455                         if nonUnit {
 456                                 x[ix] *= a[ilda+i]
 457                         }
 458                         ix -= incX
 459                 }
 460                 return
 461         }
 462         if incX == 1 {
 463                 for i := 0; i < n; i++ {
 464                         ilda := i * lda
 465                         xi := x[i]
 466                         f32.AxpyUnitary(xi, a[ilda:ilda+i], x)
 467                         if nonUnit {
 468                                 x[i] *= a[i*lda+i]
 469                         }
 470                 }
 471                 return
 472         }
 473         ix := kx
 474         for i := 0; i < n; i++ {
 475                 ilda := i * lda
 476                 xi := x[ix]
 477                 f32.AxpyInc(xi, a[ilda:ilda+i], x, uintptr(i), 1, uintptr(incX), 0, uintptr(kx))
 478                 if nonUnit {
 479                         x[ix] *= a[ilda+i]
 480                 }
 481                 ix += incX
 482         }
 483 }
 484
 485 // Strsv solves
 486 //  A * x = b if tA == blas.NoTrans
 487 //  A^T * x = b if tA == blas.Trans or blas.ConjTrans
 488 // A is an n×n triangular matrix and x is a vector.
 489 // At entry to the function, x contains the values of b, and the result is
 490 // stored in place into x.
 491 //
 492 // No test for singularity or near-singularity is included in this
 493 // routine. Such tests must be performed before calling this routine.
 494 //
 495 // Float32 implementations are autogenerated and not directly tested.
 496 func (Implementation) Strsv(ul blas.Uplo, tA blas.Transpose, d blas.Diag, n int, a []float32, lda int, x []float32, incX int) {
 497         // Test the input parameters
 498         // Verify inputs
 499         if ul != blas.Lower && ul != blas.Upper {
 500                 panic(badUplo)
 501         }
 502         if tA != blas.NoTrans && tA != blas.Trans && tA != blas.ConjTrans {
 503                 panic(badTranspose)
 504         }
 505         if d != blas.NonUnit && d != blas.Unit {
 506                 panic(badDiag)
 507         }
 508         if n < 0 {
 509                 panic(nLT0)
 510         }
 511         if lda*(n-1)+n > len(a) || lda < max(1, n) {
 512                 panic(badLdA)
 513         }
 514         if incX == 0 {
 515                 panic(zeroIncX)
 516         }
 517         if (incX > 0 && (n-1)*incX >= len(x)) || (incX < 0 && (1-n)*incX >= len(x)) {
 518                 panic(badX)
 519         }
 520         // Quick return if possible
 521         if n == 0 {
 522                 return
 523         }
 524         if n == 1 {
 525                 if d == blas.NonUnit {
 526                         x[0] /= a[0]
 527                 }
 528                 return
 529         }
 530
 531         var kx int
 532         if incX < 0 {
 533                 kx = -(n - 1) * incX
 534         }
 535         nonUnit := d == blas.NonUnit
 536         if tA == blas.NoTrans {
 537                 if ul == blas.Upper {
 538                         if incX == 1 {
 539                                 for i := n - 1; i >= 0; i-- {
 540                                         var sum float32
 541                                         atmp := a[i*lda+i+1 : i*lda+n]
 542                                         for j, v := range atmp {
 543                                                 jv := i + j + 1
 544                                                 sum += x[jv] * v
 545                                         }
 546                                         x[i] -= sum
 547                                         if nonUnit {
 548                                                 x[i] /= a[i*lda+i]
 549                                         }
 550                                 }
 551                                 return
 552                         }
 553                         ix := kx + (n-1)*incX
 554                         for i := n - 1; i >= 0; i-- {
 555                                 var sum float32
 556                                 jx := ix + incX
 557                                 atmp := a[i*lda+i+1 : i*lda+n]
 558                                 for _, v := range atmp {
 559                                         sum += x[jx] * v
 560                                         jx += incX
 561                                 }
 562                                 x[ix] -= sum
 563                                 if nonUnit {
 564                                         x[ix] /= a[i*lda+i]
 565                                 }
 566                                 ix -= incX
 567                         }
 568                         return
 569                 }
 570                 if incX == 1 {
 571                         for i := 0; i < n; i++ {
 572                                 var sum float32
 573                                 atmp := a[i*lda : i*lda+i]
 574                                 for j, v := range atmp {
 575                                         sum += x[j] * v
 576                                 }
 577                                 x[i] -= sum
 578                                 if nonUnit {
 579                                         x[i] /= a[i*lda+i]
 580                                 }
 581                         }
 582                         return
 583                 }
 584                 ix := kx
 585                 for i := 0; i < n; i++ {
 586                         jx := kx
 587                         var sum float32
 588                         atmp := a[i*lda : i*lda+i]
 589                         for _, v := range atmp {
 590                                 sum += x[jx] * v
 591                                 jx += incX
 592                         }
 593                         x[ix] -= sum
 594                         if nonUnit {
 595                                 x[ix] /= a[i*lda+i]
 596                         }
 597                         ix += incX
 598                 }
 599                 return
 600         }
 601         // Cases where a is transposed.
 602         if ul == blas.Upper {
 603                 if incX == 1 {
 604                         for i := 0; i < n; i++ {
 605                                 if nonUnit {
 606                                         x[i] /= a[i*lda+i]
 607                                 }
 608                                 xi := x[i]
 609                                 atmp := a[i*lda+i+1 : i*lda+n]
 610                                 for j, v := range atmp {
 611                                         jv := j + i + 1
 612                                         x[jv] -= v * xi
 613                                 }
 614                         }
 615                         return
 616                 }
 617                 ix := kx
 618                 for i := 0; i < n; i++ {
 619                         if nonUnit {
 620                                 x[ix] /= a[i*lda+i]
 621                         }
 622                         xi := x[ix]
 623                         jx := kx + (i+1)*incX
 624                         atmp := a[i*lda+i+1 : i*lda+n]
 625                         for _, v := range atmp {
 626                                 x[jx] -= v * xi
 627                                 jx += incX
 628                         }
 629                         ix += incX
 630                 }
 631                 return
 632         }
 633         if incX == 1 {
 634                 for i := n - 1; i >= 0; i-- {
 635                         if nonUnit {
 636                                 x[i] /= a[i*lda+i]
 637                         }
 638                         xi := x[i]
 639                         atmp := a[i*lda : i*lda+i]
 640                         for j, v := range atmp {
 641                                 x[j] -= v * xi
 642                         }
 643                 }
 644                 return
 645         }
 646         ix := kx + (n-1)*incX
 647         for i := n - 1; i >= 0; i-- {
 648                 if nonUnit {
 649                         x[ix] /= a[i*lda+i]
 650                 }
 651                 xi := x[ix]
 652                 jx := kx
 653                 atmp := a[i*lda : i*lda+i]
 654                 for _, v := range atmp {
 655                         x[jx] -= v * xi
 656                         jx += incX
 657                 }
 658                 ix -= incX
 659         }
 660 }
 661
 662 // Ssymv computes
 663 //    y = alpha * A * x + beta * y,
 664 // where a is an n×n symmetric matrix, x and y are vectors, and alpha and
 665 // beta are scalars.
 666 //
 667 // Float32 implementations are autogenerated and not directly tested.
 668 func (Implementation) Ssymv(ul blas.Uplo, n int, alpha float32, a []float32, lda int, x []float32, incX int, beta float32, y []float32, incY int) {
 669         // Check inputs
 670         if ul != blas.Lower && ul != blas.Upper {
 671                 panic(badUplo)
 672         }
 673         if n < 0 {
 674                 panic(negativeN)
 675         }
 676         if lda > 1 && lda < n {
 677                 panic(badLdA)
 678         }
 679         if incX == 0 {
 680                 panic(zeroIncX)
 681         }
 682         if incY == 0 {
 683                 panic(zeroIncY)
 684         }
 685         if (incX > 0 && (n-1)*incX >= len(x)) || (incX < 0 && (1-n)*incX >= len(x)) {
 686                 panic(badX)
 687         }
 688         if (incY > 0 && (n-1)*incY >= len(y)) || (incY < 0 && (1-n)*incY >= len(y)) {
 689                 panic(badY)
 690         }
 691         if lda*(n-1)+n > len(a) || lda < max(1, n) {
 692                 panic(badLdA)
 693         }
 694         // Quick return if possible
 695         if n == 0 || (alpha == 0 && beta == 1) {
 696                 return
 697         }
 698
 699         // Set up start points
 700         var kx, ky int
 701         if incX > 0 {
 702                 kx = 0
 703         } else {
 704                 kx = -(n - 1) * incX
 705         }
 706         if incY > 0 {
 707                 ky = 0
 708         } else {
 709                 ky = -(n - 1) * incY
 710         }
 711
 712         // Form y = beta * y
 713         if beta != 1 {
 714                 if incY > 0 {
 715                         Implementation{}.Sscal(n, beta, y, incY)
 716                 } else {
 717                         Implementation{}.Sscal(n, beta, y, -incY)
 718                 }
 719         }
 720
 721         if alpha == 0 {
 722                 return
 723         }
 724
 725         if n == 1 {
 726                 y[0] += alpha * a[0] * x[0]
 727                 return
 728         }
 729
 730         if ul == blas.Upper {
 731                 if incX == 1 {
 732                         iy := ky
 733                         for i := 0; i < n; i++ {
 734                                 xv := x[i] * alpha
 735                                 sum := x[i] * a[i*lda+i]
 736                                 jy := ky + (i+1)*incY
 737                                 atmp := a[i*lda+i+1 : i*lda+n]
 738                                 for j, v := range atmp {
 739                                         jp := j + i + 1
 740                                         sum += x[jp] * v
 741                                         y[jy] += xv * v
 742                                         jy += incY
 743                                 }
 744                                 y[iy] += alpha * sum
 745                                 iy += incY
 746                         }
 747                         return
 748                 }
 749                 ix := kx
 750                 iy := ky
 751                 for i := 0; i < n; i++ {
 752                         xv := x[ix] * alpha
 753                         sum := x[ix] * a[i*lda+i]
 754                         jx := kx + (i+1)*incX
 755                         jy := ky + (i+1)*incY
 756                         atmp := a[i*lda+i+1 : i*lda+n]
 757                         for _, v := range atmp {
 758                                 sum += x[jx] * v
 759                                 y[jy] += xv * v
 760                                 jx += incX
 761                                 jy += incY
 762                         }
 763                         y[iy] += alpha * sum
 764                         ix += incX
 765                         iy += incY
 766                 }
 767                 return
 768         }
 769         // Cases where a is lower triangular.
 770         if incX == 1 {
 771                 iy := ky
 772                 for i := 0; i < n; i++ {
 773                         jy := ky
 774                         xv := alpha * x[i]
 775                         atmp := a[i*lda : i*lda+i]
 776                         var sum float32
 777                         for j, v := range atmp {
 778                                 sum += x[j] * v
 779                                 y[jy] += xv * v
 780                                 jy += incY
 781                         }
 782                         sum += x[i] * a[i*lda+i]
 783                         sum *= alpha
 784                         y[iy] += sum
 785                         iy += incY
 786                 }
 787                 return
 788         }
 789         ix := kx
 790         iy := ky
 791         for i := 0; i < n; i++ {
 792                 jx := kx
 793                 jy := ky
 794                 xv := alpha * x[ix]
 795                 atmp := a[i*lda : i*lda+i]
 796                 var sum float32
 797                 for _, v := range atmp {
 798                         sum += x[jx] * v
 799                         y[jy] += xv * v
 800                         jx += incX
 801                         jy += incY
 802                 }
 803                 sum += x[ix] * a[i*lda+i]
 804                 sum *= alpha
 805                 y[iy] += sum
 806                 ix += incX
 807                 iy += incY
 808         }
 809 }
 810
 811 // Stbmv computes
 812 //  x = A * x if tA == blas.NoTrans
 813 //  x = A^T * x if tA == blas.Trans or blas.ConjTrans
 814 // where A is an n×n triangular banded matrix with k diagonals, and x is a vector.
 815 //
 816 // Float32 implementations are autogenerated and not directly tested.
 817 func (Implementation) Stbmv(ul blas.Uplo, tA blas.Transpose, d blas.Diag, n, k int, a []float32, lda int, x []float32, incX int) {
 818         if ul != blas.Lower && ul != blas.Upper {
 819                 panic(badUplo)
 820         }
 821         if tA != blas.NoTrans && tA != blas.Trans && tA != blas.ConjTrans {
 822                 panic(badTranspose)
 823         }
 824         if d != blas.NonUnit && d != blas.Unit {
 825                 panic(badDiag)
 826         }
 827         if n < 0 {
 828                 panic(nLT0)
 829         }
 830         if k < 0 {
 831                 panic(kLT0)
 832         }
 833         if lda*(n-1)+k+1 > len(a) || lda < k+1 {
 834                 panic(badLdA)
 835         }
 836         if incX == 0 {
 837                 panic(zeroIncX)
 838         }
 839         if (incX > 0 && (n-1)*incX >= len(x)) || (incX < 0 && (1-n)*incX >= len(x)) {
 840                 panic(badX)
 841         }
 842         if n == 0 {
 843                 return
 844         }
 845         var kx int
 846         if incX <= 0 {
 847                 kx = -(n - 1) * incX
 848         } else if incX != 1 {
 849                 kx = 0
 850         }
 851
 852         nonunit := d != blas.Unit
 853
 854         if tA == blas.NoTrans {
 855                 if ul == blas.Upper {
 856                         if incX == 1 {
 857                                 for i := 0; i < n; i++ {
 858                                         u := min(1+k, n-i)
 859                                         var sum float32
 860                                         atmp := a[i*lda:]
 861                                         xtmp := x[i:]
 862                                         for j := 1; j < u; j++ {
 863                                                 sum += xtmp[j] * atmp[j]
 864                                         }
 865                                         if nonunit {
 866                                                 sum += xtmp[0] * atmp[0]
 867                                         } else {
 868                                                 sum += xtmp[0]
 869                                         }
 870                                         x[i] = sum
 871                                 }
 872                                 return
 873                         }
 874                         ix := kx
 875                         for i := 0; i < n; i++ {
 876                                 u := min(1+k, n-i)
 877                                 var sum float32
 878                                 atmp := a[i*lda:]
 879                                 jx := incX
 880                                 for j := 1; j < u; j++ {
 881                                         sum += x[ix+jx] * atmp[j]
 882                                         jx += incX
 883                                 }
 884                                 if nonunit {
 885                                         sum += x[ix] * atmp[0]
 886                                 } else {
 887                                         sum += x[ix]
 888                                 }
 889                                 x[ix] = sum
 890                                 ix += incX
 891                         }
 892                         return
 893                 }
 894                 if incX == 1 {
 895                         for i := n - 1; i >= 0; i-- {
 896                                 l := max(0, k-i)
 897                                 atmp := a[i*lda:]
 898                                 var sum float32
 899                                 for j := l; j < k; j++ {
 900                                         sum += x[i-k+j] * atmp[j]
 901                                 }
 902                                 if nonunit {
 903                                         sum += x[i] * atmp[k]
 904                                 } else {
 905                                         sum += x[i]
 906                                 }
 907                                 x[i] = sum
 908                         }
 909                         return
 910                 }
 911                 ix := kx + (n-1)*incX
 912                 for i := n - 1; i >= 0; i-- {
 913                         l := max(0, k-i)
 914                         atmp := a[i*lda:]
 915                         var sum float32
 916                         jx := l * incX
 917                         for j := l; j < k; j++ {
 918                                 sum += x[ix-k*incX+jx] * atmp[j]
 919                                 jx += incX
 920                         }
 921                         if nonunit {
 922                                 sum += x[ix] * atmp[k]
 923                         } else {
 924                                 sum += x[ix]
 925                         }
 926                         x[ix] = sum
 927                         ix -= incX
 928                 }
 929                 return
 930         }
 931         if ul == blas.Upper {
 932                 if incX == 1 {
 933                         for i := n - 1; i >= 0; i-- {
 934                                 u := k + 1
 935                                 if i < u {
 936                                         u = i + 1
 937                                 }
 938                                 var sum float32
 939                                 for j := 1; j < u; j++ {
 940                                         sum += x[i-j] * a[(i-j)*lda+j]
 941                                 }
 942                                 if nonunit {
 943                                         sum += x[i] * a[i*lda]
 944                                 } else {
 945                                         sum += x[i]
 946                                 }
 947                                 x[i] = sum
 948                         }
 949                         return
 950                 }
 951                 ix := kx + (n-1)*incX
 952                 for i := n - 1; i >= 0; i-- {
 953                         u := k + 1
 954                         if i < u {
 955                                 u = i + 1
 956                         }
 957                         var sum float32
 958                         jx := incX
 959                         for j := 1; j < u; j++ {
 960                                 sum += x[ix-jx] * a[(i-j)*lda+j]
 961                                 jx += incX
 962                         }
 963                         if nonunit {
 964                                 sum += x[ix] * a[i*lda]
 965                         } else {
 966                                 sum += x[ix]
 967                         }
 968                         x[ix] = sum
 969                         ix -= incX
 970                 }
 971                 return
 972         }
 973         if incX == 1 {
 974                 for i := 0; i < n; i++ {
 975                         u := k
 976                         if i+k >= n {
 977                                 u = n - i - 1
 978                         }
 979                         var sum float32
 980                         for j := 0; j < u; j++ {
 981                                 sum += x[i+j+1] * a[(i+j+1)*lda+k-j-1]
 982                         }
 983                         if nonunit {
 984                                 sum += x[i] * a[i*lda+k]
 985                         } else {
 986                                 sum += x[i]
 987                         }
 988                         x[i] = sum
 989                 }
 990                 return
 991         }
 992         ix := kx
 993         for i := 0; i < n; i++ {
 994                 u := k
 995                 if i+k >= n {
 996                         u = n - i - 1
 997                 }
 998                 var (
 999                         sum float32
1000                         jx  int
1001                 )
1002                 for j := 0; j < u; j++ {
1003                         sum += x[ix+jx+incX] * a[(i+j+1)*lda+k-j-1]
1004                         jx += incX
1005                 }
1006                 if nonunit {
1007                         sum += x[ix] * a[i*lda+k]
1008                 } else {
1009                         sum += x[ix]
1010                 }
1011                 x[ix] = sum
1012                 ix += incX
1013         }
1014 }
1015
1016 // Stpmv computes
1017 //  x = A * x if tA == blas.NoTrans
1018 //  x = A^T * x if tA == blas.Trans or blas.ConjTrans
1019 // where A is an n×n unit triangular matrix in packed format, and x is a vector.
1020 //
1021 // Float32 implementations are autogenerated and not directly tested.
1022 func (Implementation) Stpmv(ul blas.Uplo, tA blas.Transpose, d blas.Diag, n int, ap []float32, x []float32, incX int) {
1023         // Verify inputs
1024         if ul != blas.Lower && ul != blas.Upper {
1025                 panic(badUplo)
1026         }
1027         if tA != blas.NoTrans && tA != blas.Trans && tA != blas.ConjTrans {
1028                 panic(badTranspose)
1029         }
1030         if d != blas.NonUnit && d != blas.Unit {
1031                 panic(badDiag)
1032         }
1033         if n < 0 {
1034                 panic(nLT0)
1035         }
1036         if len(ap) < (n*(n+1))/2 {
1037                 panic(badLdA)
1038         }
1039         if incX == 0 {
1040                 panic(zeroIncX)
1041         }
1042         if (incX > 0 && (n-1)*incX >= len(x)) || (incX < 0 && (1-n)*incX >= len(x)) {
1043                 panic(badX)
1044         }
1045         if n == 0 {
1046                 return
1047         }
1048         var kx int
1049         if incX <= 0 {
1050                 kx = -(n - 1) * incX
1051         }
1052
1053         nonUnit := d == blas.NonUnit
1054         var offset int // Offset is the index of (i,i)
1055         if tA == blas.NoTrans {
1056                 if ul == blas.Upper {
1057                         if incX == 1 {
1058                                 for i := 0; i < n; i++ {
1059                                         xi := x[i]
1060                                         if nonUnit {
1061                                                 xi *= ap[offset]
1062                                         }
1063                                         atmp := ap[offset+1 : offset+n-i]
1064                                         xtmp := x[i+1:]
1065                                         for j, v := range atmp {
1066                                                 xi += v * xtmp[j]
1067                                         }
1068                                         x[i] = xi
1069                                         offset += n - i
1070                                 }
1071                                 return
1072                         }
1073                         ix := kx
1074                         for i := 0; i < n; i++ {
1075                                 xix := x[ix]
1076                                 if nonUnit {
1077                                         xix *= ap[offset]
1078                                 }
1079                                 atmp := ap[offset+1 : offset+n-i]
1080                                 jx := kx + (i+1)*incX
1081                                 for _, v := range atmp {
1082                                         xix += v * x[jx]
1083                                         jx += incX
1084                                 }
1085                                 x[ix] = xix
1086                                 offset += n - i
1087                                 ix += incX
1088                         }
1089                         return
1090                 }
1091                 if incX == 1 {
1092                         offset = n*(n+1)/2 - 1
1093                         for i := n - 1; i >= 0; i-- {
1094                                 xi := x[i]
1095                                 if nonUnit {
1096                                         xi *= ap[offset]
1097                                 }
1098                                 atmp := ap[offset-i : offset]
1099                                 for j, v := range atmp {
1100                                         xi += v * x[j]
1101                                 }
1102                                 x[i] = xi
1103                                 offset -= i + 1
1104                         }
1105                         return
1106                 }
1107                 ix := kx + (n-1)*incX
1108                 offset = n*(n+1)/2 - 1
1109                 for i := n - 1; i >= 0; i-- {
1110                         xix := x[ix]
1111                         if nonUnit {
1112                                 xix *= ap[offset]
1113                         }
1114                         atmp := ap[offset-i : offset]
1115                         jx := kx
1116                         for _, v := range atmp {
1117                                 xix += v * x[jx]
1118                                 jx += incX
1119                         }
1120                         x[ix] = xix
1121                         offset -= i + 1
1122                         ix -= incX
1123                 }
1124                 return
1125         }
1126         // Cases where ap is transposed.
1127         if ul == blas.Upper {
1128                 if incX == 1 {
1129                         offset = n*(n+1)/2 - 1
1130                         for i := n - 1; i >= 0; i-- {
1131                                 xi := x[i]
1132                                 atmp := ap[offset+1 : offset+n-i]
1133                                 xtmp := x[i+1:]
1134                                 for j, v := range atmp {
1135                                         xtmp[j] += v * xi
1136                                 }
1137                                 if nonUnit {
1138                                         x[i] *= ap[offset]
1139                                 }
1140                                 offset -= n - i + 1
1141                         }
1142                         return
1143                 }
1144                 ix := kx + (n-1)*incX
1145                 offset = n*(n+1)/2 - 1
1146                 for i := n - 1; i >= 0; i-- {
1147                         xix := x[ix]
1148                         jx := kx + (i+1)*incX
1149                         atmp := ap[offset+1 : offset+n-i]
1150                         for _, v := range atmp {
1151                                 x[jx] += v * xix
1152                                 jx += incX
1153                         }
1154                         if nonUnit {
1155                                 x[ix] *= ap[offset]
1156                         }
1157                         offset -= n - i + 1
1158                         ix -= incX
1159                 }
1160                 return
1161         }
1162         if incX == 1 {
1163                 for i := 0; i < n; i++ {
1164                         xi := x[i]
1165                         atmp := ap[offset-i : offset]
1166                         for j, v := range atmp {
1167                                 x[j] += v * xi
1168                         }
1169                         if nonUnit {
1170                                 x[i] *= ap[offset]
1171                         }
1172                         offset += i + 2
1173                 }
1174                 return
1175         }
1176         ix := kx
1177         for i := 0; i < n; i++ {
1178                 xix := x[ix]
1179                 jx := kx
1180                 atmp := ap[offset-i : offset]
1181                 for _, v := range atmp {
1182                         x[jx] += v * xix
1183                         jx += incX
1184                 }
1185                 if nonUnit {
1186                         x[ix] *= ap[offset]
1187                 }
1188                 ix += incX
1189                 offset += i + 2
1190         }
1191 }
1192
1193 // Stbsv solves
1194 //  A * x = b
1195 // where A is an n×n triangular banded matrix with k diagonals in packed format,
1196 // and x is a vector.
1197 // At entry to the function, x contains the values of b, and the result is
1198 // stored in place into x.
1199 //
1200 // No test for singularity or near-singularity is included in this
1201 // routine. Such tests must be performed before calling this routine.
1202 //
1203 // Float32 implementations are autogenerated and not directly tested.
1204 func (Implementation) Stbsv(ul blas.Uplo, tA blas.Transpose, d blas.Diag, n, k int, a []float32, lda int, x []float32, incX int) {
1205         if ul != blas.Lower && ul != blas.Upper {
1206                 panic(badUplo)
1207         }
1208         if tA != blas.NoTrans && tA != blas.Trans && tA != blas.ConjTrans {
1209                 panic(badTranspose)
1210         }
1211         if d != blas.NonUnit && d != blas.Unit {
1212                 panic(badDiag)
1213         }
1214         if n < 0 {
1215                 panic(nLT0)
1216         }
1217         if lda*(n-1)+k+1 > len(a) || lda < k+1 {
1218                 panic(badLdA)
1219         }
1220         if incX == 0 {
1221                 panic(zeroIncX)
1222         }
1223         if (incX > 0 && (n-1)*incX >= len(x)) || (incX < 0 && (1-n)*incX >= len(x)) {
1224                 panic(badX)
1225         }
1226         if n == 0 {
1227                 return
1228         }
1229         var kx int
1230         if incX < 0 {
1231                 kx = -(n - 1) * incX
1232         } else {
1233                 kx = 0
1234         }
1235         nonUnit := d == blas.NonUnit
1236         // Form x = A^-1 x.
1237         // Several cases below use subslices for speed improvement.
1238         // The incX != 1 cases usually do not because incX may be negative.
1239         if tA == blas.NoTrans {
1240                 if ul == blas.Upper {
1241                         if incX == 1 {
1242                                 for i := n - 1; i >= 0; i-- {
1243                                         bands := k
1244                                         if i+bands >= n {
1245                                                 bands = n - i - 1
1246                                         }
1247                                         atmp := a[i*lda+1:]
1248                                         xtmp := x[i+1 : i+bands+1]
1249                                         var sum float32
1250                                         for j, v := range xtmp {
1251                                                 sum += v * atmp[j]
1252                                         }
1253                                         x[i] -= sum
1254                                         if nonUnit {
1255                                                 x[i] /= a[i*lda]
1256                                         }
1257                                 }
1258                                 return
1259                         }
1260                         ix := kx + (n-1)*incX
1261                         for i := n - 1; i >= 0; i-- {
1262                                 max := k + 1
1263                                 if i+max > n {
1264                                         max = n - i
1265                                 }
1266                                 atmp := a[i*lda:]
1267                                 var (
1268                                         jx  int
1269                                         sum float32
1270                                 )
1271                                 for j := 1; j < max; j++ {
1272                                         jx += incX
1273                                         sum += x[ix+jx] * atmp[j]
1274                                 }
1275                                 x[ix] -= sum
1276                                 if nonUnit {
1277                                         x[ix] /= atmp[0]
1278                                 }
1279                                 ix -= incX
1280                         }
1281                         return
1282                 }
1283                 if incX == 1 {
1284                         for i := 0; i < n; i++ {
1285                                 bands := k
1286                                 if i-k < 0 {
1287                                         bands = i
1288                                 }
1289                                 atmp := a[i*lda+k-bands:]
1290                                 xtmp := x[i-bands : i]
1291                                 var sum float32
1292                                 for j, v := range xtmp {
1293                                         sum += v * atmp[j]
1294                                 }
1295                                 x[i] -= sum
1296                                 if nonUnit {
1297                                         x[i] /= atmp[bands]
1298                                 }
1299                         }
1300                         return
1301                 }
1302                 ix := kx
1303                 for i := 0; i < n; i++ {
1304                         bands := k
1305                         if i-k < 0 {
1306                                 bands = i
1307                         }
1308                         atmp := a[i*lda+k-bands:]
1309                         var (
1310                                 sum float32
1311                                 jx  int
1312                         )
1313                         for j := 0; j < bands; j++ {
1314                                 sum += x[ix-bands*incX+jx] * atmp[j]
1315                                 jx += incX
1316                         }
1317                         x[ix] -= sum
1318                         if nonUnit {
1319                                 x[ix] /= atmp[bands]
1320                         }
1321                         ix += incX
1322                 }
1323                 return
1324         }
1325         // Cases where a is transposed.
1326         if ul == blas.Upper {
1327                 if incX == 1 {
1328                         for i := 0; i < n; i++ {
1329                                 bands := k
1330                                 if i-k < 0 {
1331                                         bands = i
1332                                 }
1333                                 var sum float32
1334                                 for j := 0; j < bands; j++ {
1335                                         sum += x[i-bands+j] * a[(i-bands+j)*lda+bands-j]
1336                                 }
1337                                 x[i] -= sum
1338                                 if nonUnit {
1339                                         x[i] /= a[i*lda]
1340                                 }
1341                         }
1342                         return
1343                 }
1344                 ix := kx
1345                 for i := 0; i < n; i++ {
1346                         bands := k
1347                         if i-k < 0 {
1348                                 bands = i
1349                         }
1350                         var (
1351                                 sum float32
1352                                 jx  int
1353                         )
1354                         for j := 0; j < bands; j++ {
1355                                 sum += x[ix-bands*incX+jx] * a[(i-bands+j)*lda+bands-j]
1356                                 jx += incX
1357                         }
1358                         x[ix] -= sum
1359                         if nonUnit {
1360                                 x[ix] /= a[i*lda]
1361                         }
1362                         ix += incX
1363                 }
1364                 return
1365         }
1366         if incX == 1 {
1367                 for i := n - 1; i >= 0; i-- {
1368                         bands := k
1369                         if i+bands >= n {
1370                                 bands = n - i - 1
1371                         }
1372                         var sum float32
1373                         xtmp := x[i+1 : i+1+bands]
1374                         for j, v := range xtmp {
1375                                 sum += v * a[(i+j+1)*lda+k-j-1]
1376                         }
1377                         x[i] -= sum
1378                         if nonUnit {
1379                                 x[i] /= a[i*lda+k]
1380                         }
1381                 }
1382                 return
1383         }
1384         ix := kx + (n-1)*incX
1385         for i := n - 1; i >= 0; i-- {
1386                 bands := k
1387                 if i+bands >= n {
1388                         bands = n - i - 1
1389                 }
1390                 var (
1391                         sum float32
1392                         jx  int
1393                 )
1394                 for j := 0; j < bands; j++ {
1395                         sum += x[ix+jx+incX] * a[(i+j+1)*lda+k-j-1]
1396                         jx += incX
1397                 }
1398                 x[ix] -= sum
1399                 if nonUnit {
1400                         x[ix] /= a[i*lda+k]
1401                 }
1402                 ix -= incX
1403         }
1404 }
1405
1406 // Ssbmv performs
1407 //  y = alpha * A * x + beta * y
1408 // where A is an n×n symmetric banded matrix, x and y are vectors, and alpha
1409 // and beta are scalars.
1410 //
1411 // Float32 implementations are autogenerated and not directly tested.
1412 func (Implementation) Ssbmv(ul blas.Uplo, n, k int, alpha float32, a []float32, lda int, x []float32, incX int, beta float32, y []float32, incY int) {
1413         if ul != blas.Lower && ul != blas.Upper {
1414                 panic(badUplo)
1415         }
1416         if n < 0 {
1417                 panic(nLT0)
1418         }
1419
1420         if incX == 0 {
1421                 panic(zeroIncX)
1422         }
1423         if incY == 0 {
1424                 panic(zeroIncY)
1425         }
1426         if (incX > 0 && (n-1)*incX >= len(x)) || (incX < 0 && (1-n)*incX >= len(x)) {
1427                 panic(badX)
1428         }
1429         if (incY > 0 && (n-1)*incY >= len(y)) || (incY < 0 && (1-n)*incY >= len(y)) {
1430                 panic(badY)
1431         }
1432         if lda*(n-1)+k+1 > len(a) || lda < k+1 {
1433                 panic(badLdA)
1434         }
1435
1436         // Quick return if possible
1437         if n == 0 || (alpha == 0 && beta == 1) {
1438                 return
1439         }
1440
1441         // Set up indexes
1442         lenX := n
1443         lenY := n
1444         var kx, ky int
1445         if incX > 0 {
1446                 kx = 0
1447         } else {
1448                 kx = -(lenX - 1) * incX
1449         }
1450         if incY > 0 {
1451                 ky = 0
1452         } else {
1453                 ky = -(lenY - 1) * incY
1454         }
1455
1456         // First form y = beta * y
1457         if incY > 0 {
1458                 Implementation{}.Sscal(lenY, beta, y, incY)
1459         } else {
1460                 Implementation{}.Sscal(lenY, beta, y, -incY)
1461         }
1462
1463         if alpha == 0 {
1464                 return
1465         }
1466
1467         if ul == blas.Upper {
1468                 if incX == 1 {
1469                         iy := ky
1470                         for i := 0; i < n; i++ {
1471                                 atmp := a[i*lda:]
1472                                 tmp := alpha * x[i]
1473                                 sum := tmp * atmp[0]
1474                                 u := min(k, n-i-1)
1475                                 jy := incY
1476                                 for j := 1; j <= u; j++ {
1477                                         v := atmp[j]
1478                                         sum += alpha * x[i+j] * v
1479                                         y[iy+jy] += tmp * v
1480                                         jy += incY
1481                                 }
1482                                 y[iy] += sum
1483                                 iy += incY
1484                         }
1485                         return
1486                 }
1487                 ix := kx
1488                 iy := ky
1489                 for i := 0; i < n; i++ {
1490                         atmp := a[i*lda:]
1491                         tmp := alpha * x[ix]
1492                         sum := tmp * atmp[0]
1493                         u := min(k, n-i-1)
1494                         jx := incX
1495                         jy := incY
1496                         for j := 1; j <= u; j++ {
1497                                 v := atmp[j]
1498                                 sum += alpha * x[ix+jx] * v
1499                                 y[iy+jy] += tmp * v
1500                                 jx += incX
1501                                 jy += incY
1502                         }
1503                         y[iy] += sum
1504                         ix += incX
1505                         iy += incY
1506                 }
1507                 return
1508         }
1509
1510         // Casses where a has bands below the diagonal.
1511         if incX == 1 {
1512                 iy := ky
1513                 for i := 0; i < n; i++ {
1514                         l := max(0, k-i)
1515                         tmp := alpha * x[i]
1516                         jy := l * incY
1517                         atmp := a[i*lda:]
1518                         for j := l; j < k; j++ {
1519                                 v := atmp[j]
1520                                 y[iy] += alpha * v * x[i-k+j]
1521                                 y[iy-k*incY+jy] += tmp * v
1522                                 jy += incY
1523                         }
1524                         y[iy] += tmp * atmp[k]
1525                         iy += incY
1526                 }
1527                 return
1528         }
1529         ix := kx
1530         iy := ky
1531         for i := 0; i < n; i++ {
1532                 l := max(0, k-i)
1533                 tmp := alpha * x[ix]
1534                 jx := l * incX
1535                 jy := l * incY
1536                 atmp := a[i*lda:]
1537                 for j := l; j < k; j++ {
1538                         v := atmp[j]
1539                         y[iy] += alpha * v * x[ix-k*incX+jx]
1540                         y[iy-k*incY+jy] += tmp * v
1541                         jx += incX
1542                         jy += incY
1543                 }
1544                 y[iy] += tmp * atmp[k]
1545                 ix += incX
1546                 iy += incY
1547         }
1548 }
1549
1550 // Ssyr performs the rank-one update
1551 //  a += alpha * x * x^T
1552 // where a is an n×n symmetric matrix, and x is a vector.
1553 //
1554 // Float32 implementations are autogenerated and not directly tested.
1555 func (Implementation) Ssyr(ul blas.Uplo, n int, alpha float32, x []float32, incX int, a []float32, lda int) {
1556         if ul != blas.Lower && ul != blas.Upper {
1557                 panic(badUplo)
1558         }
1559         if n < 0 {
1560                 panic(nLT0)
1561         }
1562         if incX == 0 {
1563                 panic(zeroIncX)
1564         }
1565         if (incX > 0 && (n-1)*incX >= len(x)) || (incX < 0 && (1-n)*incX >= len(x)) {
1566                 panic(badX)
1567         }
1568         if lda*(n-1)+n > len(a) || lda < max(1, n) {
1569                 panic(badLdA)
1570         }
1571         if alpha == 0 || n == 0 {
1572                 return
1573         }
1574
1575         lenX := n
1576         var kx int
1577         if incX > 0 {
1578                 kx = 0
1579         } else {
1580                 kx = -(lenX - 1) * incX
1581         }
1582         if ul == blas.Upper {
1583                 if incX == 1 {
1584                         for i := 0; i < n; i++ {
1585                                 tmp := x[i] * alpha
1586                                 if tmp != 0 {
1587                                         atmp := a[i*lda+i : i*lda+n]
1588                                         xtmp := x[i:n]
1589                                         for j, v := range xtmp {
1590                                                 atmp[j] += v * tmp
1591                                         }
1592                                 }
1593                         }
1594                         return
1595                 }
1596                 ix := kx
1597                 for i := 0; i < n; i++ {
1598                         tmp := x[ix] * alpha
1599                         if tmp != 0 {
1600                                 jx := ix
1601                                 atmp := a[i*lda:]
1602                                 for j := i; j < n; j++ {
1603                                         atmp[j] += x[jx] * tmp
1604                                         jx += incX
1605                                 }
1606                         }
1607                         ix += incX
1608                 }
1609                 return
1610         }
1611         // Cases where a is lower triangular.
1612         if incX == 1 {
1613                 for i := 0; i < n; i++ {
1614                         tmp := x[i] * alpha
1615                         if tmp != 0 {
1616                                 atmp := a[i*lda:]
1617                                 xtmp := x[:i+1]
1618                                 for j, v := range xtmp {
1619                                         atmp[j] += tmp * v
1620                                 }
1621                         }
1622                 }
1623                 return
1624         }
1625         ix := kx
1626         for i := 0; i < n; i++ {
1627                 tmp := x[ix] * alpha
1628                 if tmp != 0 {
1629                         atmp := a[i*lda:]
1630                         jx := kx
1631                         for j := 0; j < i+1; j++ {
1632                                 atmp[j] += tmp * x[jx]
1633                                 jx += incX
1634                         }
1635                 }
1636                 ix += incX
1637         }
1638 }
1639
1640 // Ssyr2 performs the symmetric rank-two update
1641 //  A += alpha * x * y^T + alpha * y * x^T
1642 // where A is a symmetric n×n matrix, x and y are vectors, and alpha is a scalar.
1643 //
1644 // Float32 implementations are autogenerated and not directly tested.
1645 func (Implementation) Ssyr2(ul blas.Uplo, n int, alpha float32, x []float32, incX int, y []float32, incY int, a []float32, lda int) {
1646         if ul != blas.Lower && ul != blas.Upper {
1647                 panic(badUplo)
1648         }
1649         if n < 0 {
1650                 panic(nLT0)
1651         }
1652         if incX == 0 {
1653                 panic(zeroIncX)
1654         }
1655         if incY == 0 {
1656                 panic(zeroIncY)
1657         }
1658         if (incX > 0 && (n-1)*incX >= len(x)) || (incX < 0 && (1-n)*incX >= len(x)) {
1659                 panic(badX)
1660         }
1661         if (incY > 0 && (n-1)*incY >= len(y)) || (incY < 0 && (1-n)*incY >= len(y)) {
1662                 panic(badY)
1663         }
1664         if lda*(n-1)+n > len(a) || lda < max(1, n) {
1665                 panic(badLdA)
1666         }
1667         if alpha == 0 {
1668                 return
1669         }
1670
1671         var ky, kx int
1672         if incY > 0 {
1673                 ky = 0
1674         } else {
1675                 ky = -(n - 1) * incY
1676         }
1677         if incX > 0 {
1678                 kx = 0
1679         } else {
1680                 kx = -(n - 1) * incX
1681         }
1682         if ul == blas.Upper {
1683                 if incX == 1 && incY == 1 {
1684                         for i := 0; i < n; i++ {
1685                                 xi := x[i]
1686                                 yi := y[i]
1687                                 atmp := a[i*lda:]
1688                                 for j := i; j < n; j++ {
1689                                         atmp[j] += alpha * (xi*y[j] + x[j]*yi)
1690                                 }
1691                         }
1692                         return
1693                 }
1694                 ix := kx
1695                 iy := ky
1696                 for i := 0; i < n; i++ {
1697                         jx := kx + i*incX
1698                         jy := ky + i*incY
1699                         xi := x[ix]
1700                         yi := y[iy]
1701                         atmp := a[i*lda:]
1702                         for j := i; j < n; j++ {
1703                                 atmp[j] += alpha * (xi*y[jy] + x[jx]*yi)
1704                                 jx += incX
1705                                 jy += incY
1706                         }
1707                         ix += incX
1708                         iy += incY
1709                 }
1710                 return
1711         }
1712         if incX == 1 && incY == 1 {
1713                 for i := 0; i < n; i++ {
1714                         xi := x[i]
1715                         yi := y[i]
1716                         atmp := a[i*lda:]
1717                         for j := 0; j <= i; j++ {
1718                                 atmp[j] += alpha * (xi*y[j] + x[j]*yi)
1719                         }
1720                 }
1721                 return
1722         }
1723         ix := kx
1724         iy := ky
1725         for i := 0; i < n; i++ {
1726                 jx := kx
1727                 jy := ky
1728                 xi := x[ix]
1729                 yi := y[iy]
1730                 atmp := a[i*lda:]
1731                 for j := 0; j <= i; j++ {
1732                         atmp[j] += alpha * (xi*y[jy] + x[jx]*yi)
1733                         jx += incX
1734                         jy += incY
1735                 }
1736                 ix += incX
1737                 iy += incY
1738         }
1739 }
1740
1741 // Stpsv solves
1742 //  A * x = b if tA == blas.NoTrans
1743 //  A^T * x = b if tA == blas.Trans or blas.ConjTrans
1744 // where A is an n×n triangular matrix in packed format and x is a vector.
1745 // At entry to the function, x contains the values of b, and the result is
1746 // stored in place into x.
1747 //
1748 // No test for singularity or near-singularity is included in this
1749 // routine. Such tests must be performed before calling this routine.
1750 //
1751 // Float32 implementations are autogenerated and not directly tested.
1752 func (Implementation) Stpsv(ul blas.Uplo, tA blas.Transpose, d blas.Diag, n int, ap []float32, x []float32, incX int) {
1753         // Verify inputs
1754         if ul != blas.Lower && ul != blas.Upper {
1755                 panic(badUplo)
1756         }
1757         if tA != blas.NoTrans && tA != blas.Trans && tA != blas.ConjTrans {
1758                 panic(badTranspose)
1759         }
1760         if d != blas.NonUnit && d != blas.Unit {
1761                 panic(badDiag)
1762         }
1763         if n < 0 {
1764                 panic(nLT0)
1765         }
1766         if len(ap) < (n*(n+1))/2 {
1767                 panic(badLdA)
1768         }
1769         if incX == 0 {
1770                 panic(zeroIncX)
1771         }
1772         if (incX > 0 && (n-1)*incX >= len(x)) || (incX < 0 && (1-n)*incX >= len(x)) {
1773                 panic(badX)
1774         }
1775         if n == 0 {
1776                 return
1777         }
1778         var kx int
1779         if incX <= 0 {
1780                 kx = -(n - 1) * incX
1781         }
1782
1783         nonUnit := d == blas.NonUnit
1784         var offset int // Offset is the index of (i,i)
1785         if tA == blas.NoTrans {
1786                 if ul == blas.Upper {
1787                         offset = n*(n+1)/2 - 1
1788                         if incX == 1 {
1789                                 for i := n - 1; i >= 0; i-- {
1790                                         atmp := ap[offset+1 : offset+n-i]
1791                                         xtmp := x[i+1:]
1792                                         var sum float32
1793                                         for j, v := range atmp {
1794                                                 sum += v * xtmp[j]
1795                                         }
1796                                         x[i] -= sum
1797                                         if nonUnit {
1798                                                 x[i] /= ap[offset]
1799                                         }
1800                                         offset -= n - i + 1
1801                                 }
1802                                 return
1803                         }
1804                         ix := kx + (n-1)*incX
1805                         for i := n - 1; i >= 0; i-- {
1806                                 atmp := ap[offset+1 : offset+n-i]
1807                                 jx := kx + (i+1)*incX
1808                                 var sum float32
1809                                 for _, v := range atmp {
1810                                         sum += v * x[jx]
1811                                         jx += incX
1812                                 }
1813                                 x[ix] -= sum
1814                                 if nonUnit {
1815                                         x[ix] /= ap[offset]
1816                                 }
1817                                 ix -= incX
1818                                 offset -= n - i + 1
1819                         }
1820                         return
1821                 }
1822                 if incX == 1 {
1823                         for i := 0; i < n; i++ {
1824                                 atmp := ap[offset-i : offset]
1825                                 var sum float32
1826                                 for j, v := range atmp {
1827                                         sum += v * x[j]
1828                                 }
1829                                 x[i] -= sum
1830                                 if nonUnit {
1831                                         x[i] /= ap[offset]
1832                                 }
1833                                 offset += i + 2
1834                         }
1835                         return
1836                 }
1837                 ix := kx
1838                 for i := 0; i < n; i++ {
1839                         jx := kx
1840                         atmp := ap[offset-i : offset]
1841                         var sum float32
1842                         for _, v := range atmp {
1843                                 sum += v * x[jx]
1844                                 jx += incX
1845                         }
1846                         x[ix] -= sum
1847                         if nonUnit {
1848                                 x[ix] /= ap[offset]
1849                         }
1850                         ix += incX
1851                         offset += i + 2
1852                 }
1853                 return
1854         }
1855         // Cases where ap is transposed.
1856         if ul == blas.Upper {
1857                 if incX == 1 {
1858                         for i := 0; i < n; i++ {
1859                                 if nonUnit {
1860                                         x[i] /= ap[offset]
1861                                 }
1862                                 xi := x[i]
1863                                 atmp := ap[offset+1 : offset+n-i]
1864                                 xtmp := x[i+1:]
1865                                 for j, v := range atmp {
1866                                         xtmp[j] -= v * xi
1867                                 }
1868                                 offset += n - i
1869                         }
1870                         return
1871                 }
1872                 ix := kx
1873                 for i := 0; i < n; i++ {
1874                         if nonUnit {
1875                                 x[ix] /= ap[offset]
1876                         }
1877                         xix := x[ix]
1878                         atmp := ap[offset+1 : offset+n-i]
1879                         jx := kx + (i+1)*incX
1880                         for _, v := range atmp {
1881                                 x[jx] -= v * xix
1882                                 jx += incX
1883                         }
1884                         ix += incX
1885                         offset += n - i
1886                 }
1887                 return
1888         }
1889         if incX == 1 {
1890                 offset = n*(n+1)/2 - 1
1891                 for i := n - 1; i >= 0; i-- {
1892                         if nonUnit {
1893                                 x[i] /= ap[offset]
1894                         }
1895                         xi := x[i]
1896                         atmp := ap[offset-i : offset]
1897                         for j, v := range atmp {
1898                                 x[j] -= v * xi
1899                         }
1900                         offset -= i + 1
1901                 }
1902                 return
1903         }
1904         ix := kx + (n-1)*incX
1905         offset = n*(n+1)/2 - 1
1906         for i := n - 1; i >= 0; i-- {
1907                 if nonUnit {
1908                         x[ix] /= ap[offset]
1909                 }
1910                 xix := x[ix]
1911                 atmp := ap[offset-i : offset]
1912                 jx := kx
1913                 for _, v := range atmp {
1914                         x[jx] -= v * xix
1915                         jx += incX
1916                 }
1917                 ix -= incX
1918                 offset -= i + 1
1919         }
1920 }
1921
1922 // Sspmv performs
1923 //    y = alpha * A * x + beta * y,
1924 // where A is an n×n symmetric matrix in packed format, x and y are vectors
1925 // and alpha and beta are scalars.
1926 //
1927 // Float32 implementations are autogenerated and not directly tested.
1928 func (Implementation) Sspmv(ul blas.Uplo, n int, alpha float32, a []float32, x []float32, incX int, beta float32, y []float32, incY int) {
1929         // Verify inputs
1930         if ul != blas.Lower && ul != blas.Upper {
1931                 panic(badUplo)
1932         }
1933         if n < 0 {
1934                 panic(nLT0)
1935         }
1936         if len(a) < (n*(n+1))/2 {
1937                 panic(badLdA)
1938         }
1939         if incX == 0 {
1940                 panic(zeroIncX)
1941         }
1942         if incY == 0 {
1943                 panic(zeroIncY)
1944         }
1945         if (incX > 0 && (n-1)*incX >= len(x)) || (incX < 0 && (1-n)*incX >= len(x)) {
1946                 panic(badX)
1947         }
1948         if (incY > 0 && (n-1)*incY >= len(y)) || (incY < 0 && (1-n)*incY >= len(y)) {
1949                 panic(badY)
1950         }
1951         // Quick return if possible
1952         if n == 0 || (alpha == 0 && beta == 1) {
1953                 return
1954         }
1955
1956         // Set up start points
1957         var kx, ky int
1958         if incX > 0 {
1959                 kx = 0
1960         } else {
1961                 kx = -(n - 1) * incX
1962         }
1963         if incY > 0 {
1964                 ky = 0
1965         } else {
1966                 ky = -(n - 1) * incY
1967         }
1968
1969         // Form y = beta * y
1970         if beta != 1 {
1971                 if incY > 0 {
1972                         Implementation{}.Sscal(n, beta, y, incY)
1973                 } else {
1974                         Implementation{}.Sscal(n, beta, y, -incY)
1975                 }
1976         }
1977
1978         if alpha == 0 {
1979                 return
1980         }
1981
1982         if n == 1 {
1983                 y[0] += alpha * a[0] * x[0]
1984                 return
1985         }
1986         var offset int // Offset is the index of (i,i).
1987         if ul == blas.Upper {
1988                 if incX == 1 {
1989                         iy := ky
1990                         for i := 0; i < n; i++ {
1991                                 xv := x[i] * alpha
1992                                 sum := a[offset] * x[i]
1993                                 atmp := a[offset+1 : offset+n-i]
1994                                 xtmp := x[i+1:]
1995                                 jy := ky + (i+1)*incY
1996                                 for j, v := range atmp {
1997                                         sum += v * xtmp[j]
1998                                         y[jy] += v * xv
1999                                         jy += incY
2000                                 }
2001                                 y[iy] += alpha * sum
2002                                 iy += incY
2003                                 offset += n - i
2004                         }
2005                         return
2006                 }
2007                 ix := kx
2008                 iy := ky
2009                 for i := 0; i < n; i++ {
2010                         xv := x[ix] * alpha
2011                         sum := a[offset] * x[ix]
2012                         atmp := a[offset+1 : offset+n-i]
2013                         jx := kx + (i+1)*incX
2014                         jy := ky + (i+1)*incY
2015                         for _, v := range atmp {
2016                                 sum += v * x[jx]
2017                                 y[jy] += v * xv
2018                                 jx += incX
2019                                 jy += incY
2020                         }
2021                         y[iy] += alpha * sum
2022                         ix += incX
2023                         iy += incY
2024                         offset += n - i
2025                 }
2026                 return
2027         }
2028         if incX == 1 {
2029                 iy := ky
2030                 for i := 0; i < n; i++ {
2031                         xv := x[i] * alpha
2032                         atmp := a[offset-i : offset]
2033                         jy := ky
2034                         var sum float32
2035                         for j, v := range atmp {
2036                                 sum += v * x[j]
2037                                 y[jy] += v * xv
2038                                 jy += incY
2039                         }
2040                         sum += a[offset] * x[i]
2041                         y[iy] += alpha * sum
2042                         iy += incY
2043                         offset += i + 2
2044                 }
2045                 return
2046         }
2047         ix := kx
2048         iy := ky
2049         for i := 0; i < n; i++ {
2050                 xv := x[ix] * alpha
2051                 atmp := a[offset-i : offset]
2052                 jx := kx
2053                 jy := ky
2054                 var sum float32
2055                 for _, v := range atmp {
2056                         sum += v * x[jx]
2057                         y[jy] += v * xv
2058                         jx += incX
2059                         jy += incY
2060                 }
2061
2062                 sum += a[offset] * x[ix]
2063                 y[iy] += alpha * sum
2064                 ix += incX
2065                 iy += incY
2066                 offset += i + 2
2067         }
2068 }
2069
2070 // Sspr computes the rank-one operation
2071 //  a += alpha * x * x^T
2072 // where a is an n×n symmetric matrix in packed format, x is a vector, and
2073 // alpha is a scalar.
2074 //
2075 // Float32 implementations are autogenerated and not directly tested.
2076 func (Implementation) Sspr(ul blas.Uplo, n int, alpha float32, x []float32, incX int, a []float32) {
2077         if ul != blas.Lower && ul != blas.Upper {
2078                 panic(badUplo)
2079         }
2080         if n < 0 {
2081                 panic(nLT0)
2082         }
2083         if incX == 0 {
2084                 panic(zeroIncX)
2085         }
2086         if (incX > 0 && (n-1)*incX >= len(x)) || (incX < 0 && (1-n)*incX >= len(x)) {
2087                 panic(badX)
2088         }
2089         if len(a) < (n*(n+1))/2 {
2090                 panic(badLdA)
2091         }
2092         if alpha == 0 || n == 0 {
2093                 return
2094         }
2095         lenX := n
2096         var kx int
2097         if incX > 0 {
2098                 kx = 0
2099         } else {
2100                 kx = -(lenX - 1) * incX
2101         }
2102         var offset int // Offset is the index of (i,i).
2103         if ul == blas.Upper {
2104                 if incX == 1 {
2105                         for i := 0; i < n; i++ {
2106                                 atmp := a[offset:]
2107                                 xv := alpha * x[i]
2108                                 xtmp := x[i:n]
2109                                 for j, v := range xtmp {
2110                                         atmp[j] += xv * v
2111                                 }
2112                                 offset += n - i
2113                         }
2114                         return
2115                 }
2116                 ix := kx
2117                 for i := 0; i < n; i++ {
2118                         jx := kx + i*incX
2119                         atmp := a[offset:]
2120                         xv := alpha * x[ix]
2121                         for j := 0; j < n-i; j++ {
2122                                 atmp[j] += xv * x[jx]
2123                                 jx += incX
2124                         }
2125                         ix += incX
2126                         offset += n - i
2127                 }
2128                 return
2129         }
2130         if incX == 1 {
2131                 for i := 0; i < n; i++ {
2132                         atmp := a[offset-i:]
2133                         xv := alpha * x[i]
2134                         xtmp := x[:i+1]
2135                         for j, v := range xtmp {
2136                                 atmp[j] += xv * v
2137                         }
2138                         offset += i + 2
2139                 }
2140                 return
2141         }
2142         ix := kx
2143         for i := 0; i < n; i++ {
2144                 jx := kx
2145                 atmp := a[offset-i:]
2146                 xv := alpha * x[ix]
2147                 for j := 0; j <= i; j++ {
2148                         atmp[j] += xv * x[jx]
2149                         jx += incX
2150                 }
2151                 ix += incX
2152                 offset += i + 2
2153         }
2154 }
2155
2156 // Sspr2 performs the symmetric rank-2 update
2157 //  A += alpha * x * y^T + alpha * y * x^T,
2158 // where A is an n×n symmetric matrix in packed format, x and y are vectors,
2159 // and alpha is a scalar.
2160 //
2161 // Float32 implementations are autogenerated and not directly tested.
2162 func (Implementation) Sspr2(ul blas.Uplo, n int, alpha float32, x []float32, incX int, y []float32, incY int, ap []float32) {
2163         if ul != blas.Lower && ul != blas.Upper {
2164                 panic(badUplo)
2165         }
2166         if n < 0 {
2167                 panic(nLT0)
2168         }
2169         if incX == 0 {
2170                 panic(zeroIncX)
2171         }
2172         if incY == 0 {
2173                 panic(zeroIncY)
2174         }
2175         if (incX > 0 && (n-1)*incX >= len(x)) || (incX < 0 && (1-n)*incX >= len(x)) {
2176                 panic(badX)
2177         }
2178         if (incY > 0 && (n-1)*incY >= len(y)) || (incY < 0 && (1-n)*incY >= len(y)) {
2179                 panic(badY)
2180         }
2181         if len(ap) < (n*(n+1))/2 {
2182                 panic(badLdA)
2183         }
2184         if alpha == 0 {
2185                 return
2186         }
2187         var ky, kx int
2188         if incY > 0 {
2189                 ky = 0
2190         } else {
2191                 ky = -(n - 1) * incY
2192         }
2193         if incX > 0 {
2194                 kx = 0
2195         } else {
2196                 kx = -(n - 1) * incX
2197         }
2198         var offset int // Offset is the index of (i,i).
2199         if ul == blas.Upper {
2200                 if incX == 1 && incY == 1 {
2201                         for i := 0; i < n; i++ {
2202                                 atmp := ap[offset:]
2203                                 xi := x[i]
2204                                 yi := y[i]
2205                                 xtmp := x[i:n]
2206                                 ytmp := y[i:n]
2207                                 for j, v := range xtmp {
2208                                         atmp[j] += alpha * (xi*ytmp[j] + v*yi)
2209                                 }
2210                                 offset += n - i
2211                         }
2212                         return
2213                 }
2214                 ix := kx
2215                 iy := ky
2216                 for i := 0; i < n; i++ {
2217                         jx := kx + i*incX
2218                         jy := ky + i*incY
2219                         atmp := ap[offset:]
2220                         xi := x[ix]
2221                         yi := y[iy]
2222                         for j := 0; j < n-i; j++ {
2223                                 atmp[j] += alpha * (xi*y[jy] + x[jx]*yi)
2224                                 jx += incX
2225                                 jy += incY
2226                         }
2227                         ix += incX
2228                         iy += incY
2229                         offset += n - i
2230                 }
2231                 return
2232         }
2233         if incX == 1 && incY == 1 {
2234                 for i := 0; i < n; i++ {
2235                         atmp := ap[offset-i:]
2236                         xi := x[i]
2237                         yi := y[i]
2238                         xtmp := x[:i+1]
2239                         for j, v := range xtmp {
2240                                 atmp[j] += alpha * (xi*y[j] + v*yi)
2241                         }
2242                         offset += i + 2
2243                 }
2244                 return
2245         }
2246         ix := kx
2247         iy := ky
2248         for i := 0; i < n; i++ {
2249                 jx := kx
2250                 jy := ky
2251                 atmp := ap[offset-i:]
2252                 for j := 0; j <= i; j++ {
2253                         atmp[j] += alpha * (x[ix]*y[jy] + x[jx]*y[iy])
2254                         jx += incX
2255                         jy += incY
2256                 }
2257                 ix += incX
2258                 iy += incY
2259                 offset += i + 2
2260         }
2261 }