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

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