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

   1 // Code generated by "go generate gonum.org/v1/gonum/blas/gonum”; DO NOT EDIT.
   2
   3 // Copyright ©2015 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         math "gonum.org/v1/gonum/internal/math32"
  11
  12         "gonum.org/v1/gonum/blas"
  13         "gonum.org/v1/gonum/internal/asm/f32"
  14 )
  15
  16 var _ blas.Float32Level1 = Implementation{}
  17
  18 // Snrm2 computes the Euclidean norm of a vector,
  19 //  sqrt(\sum_i x[i] * x[i]).
  20 // This function returns 0 if incX is negative.
  21 //
  22 // Float32 implementations are autogenerated and not directly tested.
  23 func (Implementation) Snrm2(n int, x []float32, incX int) float32 {
  24         if incX < 1 {
  25                 if incX == 0 {
  26                         panic(zeroIncX)
  27                 }
  28                 return 0
  29         }
  30         if incX > 0 && (n-1)*incX >= len(x) {
  31                 panic(badX)
  32         }
  33         if n < 2 {
  34                 if n == 1 {
  35                         return math.Abs(x[0])
  36                 }
  37                 if n == 0 {
  38                         return 0
  39                 }
  40                 if n < 1 {
  41                         panic(negativeN)
  42                 }
  43         }
  44         var (
  45                 scale      float32 = 0
  46                 sumSquares float32 = 1
  47         )
  48         if incX == 1 {
  49                 x = x[:n]
  50                 for _, v := range x {
  51                         if v == 0 {
  52                                 continue
  53                         }
  54                         absxi := math.Abs(v)
  55                         if math.IsNaN(absxi) {
  56                                 return math.NaN()
  57                         }
  58                         if scale < absxi {
  59                                 sumSquares = 1 + sumSquares*(scale/absxi)*(scale/absxi)
  60                                 scale = absxi
  61                         } else {
  62                                 sumSquares = sumSquares + (absxi/scale)*(absxi/scale)
  63                         }
  64                 }
  65                 if math.IsInf(scale, 1) {
  66                         return math.Inf(1)
  67                 }
  68                 return scale * math.Sqrt(sumSquares)
  69         }
  70         for ix := 0; ix < n*incX; ix += incX {
  71                 val := x[ix]
  72                 if val == 0 {
  73                         continue
  74                 }
  75                 absxi := math.Abs(val)
  76                 if math.IsNaN(absxi) {
  77                         return math.NaN()
  78                 }
  79                 if scale < absxi {
  80                         sumSquares = 1 + sumSquares*(scale/absxi)*(scale/absxi)
  81                         scale = absxi
  82                 } else {
  83                         sumSquares = sumSquares + (absxi/scale)*(absxi/scale)
  84                 }
  85         }
  86         if math.IsInf(scale, 1) {
  87                 return math.Inf(1)
  88         }
  89         return scale * math.Sqrt(sumSquares)
  90 }
  91
  92 // Sasum computes the sum of the absolute values of the elements of x.
  93 //  \sum_i |x[i]|
  94 // Sasum returns 0 if incX is negative.
  95 //
  96 // Float32 implementations are autogenerated and not directly tested.
  97 func (Implementation) Sasum(n int, x []float32, incX int) float32 {
  98         var sum float32
  99         if n < 0 {
 100                 panic(negativeN)
 101         }
 102         if incX < 1 {
 103                 if incX == 0 {
 104                         panic(zeroIncX)
 105                 }
 106                 return 0
 107         }
 108         if incX > 0 && (n-1)*incX >= len(x) {
 109                 panic(badX)
 110         }
 111         if incX == 1 {
 112                 x = x[:n]
 113                 for _, v := range x {
 114                         sum += math.Abs(v)
 115                 }
 116                 return sum
 117         }
 118         for i := 0; i < n; i++ {
 119                 sum += math.Abs(x[i*incX])
 120         }
 121         return sum
 122 }
 123
 124 // Isamax returns the index of an element of x with the largest absolute value.
 125 // If there are multiple such indices the earliest is returned.
 126 // Isamax returns -1 if n == 0.
 127 //
 128 // Float32 implementations are autogenerated and not directly tested.
 129 func (Implementation) Isamax(n int, x []float32, incX int) int {
 130         if incX < 1 {
 131                 if incX == 0 {
 132                         panic(zeroIncX)
 133                 }
 134                 return -1
 135         }
 136         if incX > 0 && (n-1)*incX >= len(x) {
 137                 panic(badX)
 138         }
 139         if n < 2 {
 140                 if n == 1 {
 141                         return 0
 142                 }
 143                 if n == 0 {
 144                         return -1 // Netlib returns invalid index when n == 0
 145                 }
 146                 if n < 1 {
 147                         panic(negativeN)
 148                 }
 149         }
 150         idx := 0
 151         max := math.Abs(x[0])
 152         if incX == 1 {
 153                 for i, v := range x[:n] {
 154                         absV := math.Abs(v)
 155                         if absV > max {
 156                                 max = absV
 157                                 idx = i
 158                         }
 159                 }
 160                 return idx
 161         }
 162         ix := incX
 163         for i := 1; i < n; i++ {
 164                 v := x[ix]
 165                 absV := math.Abs(v)
 166                 if absV > max {
 167                         max = absV
 168                         idx = i
 169                 }
 170                 ix += incX
 171         }
 172         return idx
 173 }
 174
 175 // Sswap exchanges the elements of two vectors.
 176 //  x[i], y[i] = y[i], x[i] for all i
 177 //
 178 // Float32 implementations are autogenerated and not directly tested.
 179 func (Implementation) Sswap(n int, x []float32, incX int, y []float32, incY int) {
 180         if incX == 0 {
 181                 panic(zeroIncX)
 182         }
 183         if incY == 0 {
 184                 panic(zeroIncY)
 185         }
 186         if n < 1 {
 187                 if n == 0 {
 188                         return
 189                 }
 190                 panic(negativeN)
 191         }
 192         if (incX > 0 && (n-1)*incX >= len(x)) || (incX < 0 && (1-n)*incX >= len(x)) {
 193                 panic(badX)
 194         }
 195         if (incY > 0 && (n-1)*incY >= len(y)) || (incY < 0 && (1-n)*incY >= len(y)) {
 196                 panic(badY)
 197         }
 198         if incX == 1 && incY == 1 {
 199                 x = x[:n]
 200                 for i, v := range x {
 201                         x[i], y[i] = y[i], v
 202                 }
 203                 return
 204         }
 205         var ix, iy int
 206         if incX < 0 {
 207                 ix = (-n + 1) * incX
 208         }
 209         if incY < 0 {
 210                 iy = (-n + 1) * incY
 211         }
 212         for i := 0; i < n; i++ {
 213                 x[ix], y[iy] = y[iy], x[ix]
 214                 ix += incX
 215                 iy += incY
 216         }
 217 }
 218
 219 // Scopy copies the elements of x into the elements of y.
 220 //  y[i] = x[i] for all i
 221 //
 222 // Float32 implementations are autogenerated and not directly tested.
 223 func (Implementation) Scopy(n int, x []float32, incX int, y []float32, incY int) {
 224         if incX == 0 {
 225                 panic(zeroIncX)
 226         }
 227         if incY == 0 {
 228                 panic(zeroIncY)
 229         }
 230         if n < 1 {
 231                 if n == 0 {
 232                         return
 233                 }
 234                 panic(negativeN)
 235         }
 236         if (incX > 0 && (n-1)*incX >= len(x)) || (incX < 0 && (1-n)*incX >= len(x)) {
 237                 panic(badX)
 238         }
 239         if (incY > 0 && (n-1)*incY >= len(y)) || (incY < 0 && (1-n)*incY >= len(y)) {
 240                 panic(badY)
 241         }
 242         if incX == 1 && incY == 1 {
 243                 copy(y[:n], x[:n])
 244                 return
 245         }
 246         var ix, iy int
 247         if incX < 0 {
 248                 ix = (-n + 1) * incX
 249         }
 250         if incY < 0 {
 251                 iy = (-n + 1) * incY
 252         }
 253         for i := 0; i < n; i++ {
 254                 y[iy] = x[ix]
 255                 ix += incX
 256                 iy += incY
 257         }
 258 }
 259
 260 // Saxpy adds alpha times x to y
 261 //  y[i] += alpha * x[i] for all i
 262 //
 263 // Float32 implementations are autogenerated and not directly tested.
 264 func (Implementation) Saxpy(n int, alpha float32, x []float32, incX int, y []float32, incY int) {
 265         if incX == 0 {
 266                 panic(zeroIncX)
 267         }
 268         if incY == 0 {
 269                 panic(zeroIncY)
 270         }
 271         if n < 1 {
 272                 if n == 0 {
 273                         return
 274                 }
 275                 panic(negativeN)
 276         }
 277         if (incX > 0 && (n-1)*incX >= len(x)) || (incX < 0 && (1-n)*incX >= len(x)) {
 278                 panic(badX)
 279         }
 280         if (incY > 0 && (n-1)*incY >= len(y)) || (incY < 0 && (1-n)*incY >= len(y)) {
 281                 panic(badY)
 282         }
 283         if alpha == 0 {
 284                 return
 285         }
 286         if incX == 1 && incY == 1 {
 287                 f32.AxpyUnitary(alpha, x[:n], y[:n])
 288                 return
 289         }
 290         var ix, iy int
 291         if incX < 0 {
 292                 ix = (-n + 1) * incX
 293         }
 294         if incY < 0 {
 295                 iy = (-n + 1) * incY
 296         }
 297         f32.AxpyInc(alpha, x, y, uintptr(n), uintptr(incX), uintptr(incY), uintptr(ix), uintptr(iy))
 298 }
 299
 300 // Srotg computes the plane rotation
 301 //   _    _      _ _       _ _
 302 //  |  c s |    | a |     | r |
 303 //  | -s c |  * | b |   = | 0 |
 304 //   ‾    ‾      ‾ ‾       ‾ ‾
 305 // where
 306 //  r = ±√(a^2 + b^2)
 307 //  c = a/r, the cosine of the plane rotation
 308 //  s = b/r, the sine of the plane rotation
 309 //
 310 // NOTE: There is a discrepancy between the refence implementation and the BLAS
 311 // technical manual regarding the sign for r when a or b are zero.
 312 // Srotg agrees with the definition in the manual and other
 313 // common BLAS implementations.
 314 //
 315 // Float32 implementations are autogenerated and not directly tested.
 316 func (Implementation) Srotg(a, b float32) (c, s, r, z float32) {
 317         if b == 0 && a == 0 {
 318                 return 1, 0, a, 0
 319         }
 320         absA := math.Abs(a)
 321         absB := math.Abs(b)
 322         aGTb := absA > absB
 323         r = math.Hypot(a, b)
 324         if aGTb {
 325                 r = math.Copysign(r, a)
 326         } else {
 327                 r = math.Copysign(r, b)
 328         }
 329         c = a / r
 330         s = b / r
 331         if aGTb {
 332                 z = s
 333         } else if c != 0 { // r == 0 case handled above
 334                 z = 1 / c
 335         } else {
 336                 z = 1
 337         }
 338         return
 339 }
 340
 341 // Srotmg computes the modified Givens rotation. See
 342 // http://www.netlib.org/lapack/explore-html/df/deb/drotmg_8f.html
 343 // for more details.
 344 //
 345 // Float32 implementations are autogenerated and not directly tested.
 346 func (Implementation) Srotmg(d1, d2, x1, y1 float32) (p blas.SrotmParams, rd1, rd2, rx1 float32) {
 347         var p1, p2, q1, q2, u float32
 348
 349         const (
 350                 gam    = 4096.0
 351                 gamsq  = 16777216.0
 352                 rgamsq = 5.9604645e-8
 353         )
 354
 355         if d1 < 0 {
 356                 p.Flag = blas.Rescaling
 357                 return
 358         }
 359
 360         p2 = d2 * y1
 361         if p2 == 0 {
 362                 p.Flag = blas.Identity
 363                 rd1 = d1
 364                 rd2 = d2
 365                 rx1 = x1
 366                 return
 367         }
 368         p1 = d1 * x1
 369         q2 = p2 * y1
 370         q1 = p1 * x1
 371
 372         absQ1 := math.Abs(q1)
 373         absQ2 := math.Abs(q2)
 374
 375         if absQ1 < absQ2 && q2 < 0 {
 376                 p.Flag = blas.Rescaling
 377                 return
 378         }
 379
 380         if d1 == 0 {
 381                 p.Flag = blas.Diagonal
 382                 p.H[0] = p1 / p2
 383                 p.H[3] = x1 / y1
 384                 u = 1 + p.H[0]*p.H[3]
 385                 rd1, rd2 = d2/u, d1/u
 386                 rx1 = y1 / u
 387                 return
 388         }
 389
 390         // Now we know that d1 != 0, and d2 != 0. If d2 == 0, it would be caught
 391         // when p2 == 0, and if d1 == 0, then it is caught above
 392
 393         if absQ1 > absQ2 {
 394                 p.H[1] = -y1 / x1
 395                 p.H[2] = p2 / p1
 396                 u = 1 - p.H[2]*p.H[1]
 397                 rd1 = d1
 398                 rd2 = d2
 399                 rx1 = x1
 400                 p.Flag = blas.OffDiagonal
 401                 // u must be greater than zero because |q1| > |q2|, so check from netlib
 402                 // is unnecessary
 403                 // This is left in for ease of comparison with complex routines
 404                 //if u > 0 {
 405                 rd1 /= u
 406                 rd2 /= u
 407                 rx1 *= u
 408                 //}
 409         } else {
 410                 p.Flag = blas.Diagonal
 411                 p.H[0] = p1 / p2
 412                 p.H[3] = x1 / y1
 413                 u = 1 + p.H[0]*p.H[3]
 414                 rd1 = d2 / u
 415                 rd2 = d1 / u
 416                 rx1 = y1 * u
 417         }
 418
 419         for rd1 <= rgamsq || rd1 >= gamsq {
 420                 if p.Flag == blas.OffDiagonal {
 421                         p.H[0] = 1
 422                         p.H[3] = 1
 423                         p.Flag = blas.Rescaling
 424                 } else if p.Flag == blas.Diagonal {
 425                         p.H[1] = -1
 426                         p.H[2] = 1
 427                         p.Flag = blas.Rescaling
 428                 }
 429                 if rd1 <= rgamsq {
 430                         rd1 *= gam * gam
 431                         rx1 /= gam
 432                         p.H[0] /= gam
 433                         p.H[2] /= gam
 434                 } else {
 435                         rd1 /= gam * gam
 436                         rx1 *= gam
 437                         p.H[0] *= gam
 438                         p.H[2] *= gam
 439                 }
 440         }
 441
 442         for math.Abs(rd2) <= rgamsq || math.Abs(rd2) >= gamsq {
 443                 if p.Flag == blas.OffDiagonal {
 444                         p.H[0] = 1
 445                         p.H[3] = 1
 446                         p.Flag = blas.Rescaling
 447                 } else if p.Flag == blas.Diagonal {
 448                         p.H[1] = -1
 449                         p.H[2] = 1
 450                         p.Flag = blas.Rescaling
 451                 }
 452                 if math.Abs(rd2) <= rgamsq {
 453                         rd2 *= gam * gam
 454                         p.H[1] /= gam
 455                         p.H[3] /= gam
 456                 } else {
 457                         rd2 /= gam * gam
 458                         p.H[1] *= gam
 459                         p.H[3] *= gam
 460                 }
 461         }
 462         return
 463 }
 464
 465 // Srot applies a plane transformation.
 466 //  x[i] = c * x[i] + s * y[i]
 467 //  y[i] = c * y[i] - s * x[i]
 468 //
 469 // Float32 implementations are autogenerated and not directly tested.
 470 func (Implementation) Srot(n int, x []float32, incX int, y []float32, incY int, c float32, s float32) {
 471         if incX == 0 {
 472                 panic(zeroIncX)
 473         }
 474         if incY == 0 {
 475                 panic(zeroIncY)
 476         }
 477         if n < 1 {
 478                 if n == 0 {
 479                         return
 480                 }
 481                 panic(negativeN)
 482         }
 483         if (incX > 0 && (n-1)*incX >= len(x)) || (incX < 0 && (1-n)*incX >= len(x)) {
 484                 panic(badX)
 485         }
 486         if (incY > 0 && (n-1)*incY >= len(y)) || (incY < 0 && (1-n)*incY >= len(y)) {
 487                 panic(badY)
 488         }
 489         if incX == 1 && incY == 1 {
 490                 x = x[:n]
 491                 for i, vx := range x {
 492                         vy := y[i]
 493                         x[i], y[i] = c*vx+s*vy, c*vy-s*vx
 494                 }
 495                 return
 496         }
 497         var ix, iy int
 498         if incX < 0 {
 499                 ix = (-n + 1) * incX
 500         }
 501         if incY < 0 {
 502                 iy = (-n + 1) * incY
 503         }
 504         for i := 0; i < n; i++ {
 505                 vx := x[ix]
 506                 vy := y[iy]
 507                 x[ix], y[iy] = c*vx+s*vy, c*vy-s*vx
 508                 ix += incX
 509                 iy += incY
 510         }
 511 }
 512
 513 // Srotm applies the modified Givens rotation to the 2×n matrix.
 514 //
 515 // Float32 implementations are autogenerated and not directly tested.
 516 func (Implementation) Srotm(n int, x []float32, incX int, y []float32, incY int, p blas.SrotmParams) {
 517         if incX == 0 {
 518                 panic(zeroIncX)
 519         }
 520         if incY == 0 {
 521                 panic(zeroIncY)
 522         }
 523         if n <= 0 {
 524                 if n == 0 {
 525                         return
 526                 }
 527                 panic(negativeN)
 528         }
 529         if (incX > 0 && (n-1)*incX >= len(x)) || (incX < 0 && (1-n)*incX >= len(x)) {
 530                 panic(badX)
 531         }
 532         if (incY > 0 && (n-1)*incY >= len(y)) || (incY < 0 && (1-n)*incY >= len(y)) {
 533                 panic(badY)
 534         }
 535
 536         var h11, h12, h21, h22 float32
 537         var ix, iy int
 538         switch p.Flag {
 539         case blas.Identity:
 540                 return
 541         case blas.Rescaling:
 542                 h11 = p.H[0]
 543                 h12 = p.H[2]
 544                 h21 = p.H[1]
 545                 h22 = p.H[3]
 546         case blas.OffDiagonal:
 547                 h11 = 1
 548                 h12 = p.H[2]
 549                 h21 = p.H[1]
 550                 h22 = 1
 551         case blas.Diagonal:
 552                 h11 = p.H[0]
 553                 h12 = 1
 554                 h21 = -1
 555                 h22 = p.H[3]
 556         }
 557         if incX < 0 {
 558                 ix = (-n + 1) * incX
 559         }
 560         if incY < 0 {
 561                 iy = (-n + 1) * incY
 562         }
 563         if incX == 1 && incY == 1 {
 564                 x = x[:n]
 565                 for i, vx := range x {
 566                         vy := y[i]
 567                         x[i], y[i] = vx*h11+vy*h12, vx*h21+vy*h22
 568                 }
 569                 return
 570         }
 571         for i := 0; i < n; i++ {
 572                 vx := x[ix]
 573                 vy := y[iy]
 574                 x[ix], y[iy] = vx*h11+vy*h12, vx*h21+vy*h22
 575                 ix += incX
 576                 iy += incY
 577         }
 578 }
 579
 580 // Sscal scales x by alpha.
 581 //  x[i] *= alpha
 582 // Sscal has no effect if incX < 0.
 583 //
 584 // Float32 implementations are autogenerated and not directly tested.
 585 func (Implementation) Sscal(n int, alpha float32, x []float32, incX int) {
 586         if incX < 1 {
 587                 if incX == 0 {
 588                         panic(zeroIncX)
 589                 }
 590                 return
 591         }
 592         if (n-1)*incX >= len(x) {
 593                 panic(badX)
 594         }
 595         if n < 1 {
 596                 if n == 0 {
 597                         return
 598                 }
 599                 panic(negativeN)
 600         }
 601         if alpha == 0 {
 602                 if incX == 1 {
 603                         x = x[:n]
 604                         for i := range x {
 605                                 x[i] = 0
 606                         }
 607                         return
 608                 }
 609                 for ix := 0; ix < n*incX; ix += incX {
 610                         x[ix] = 0
 611                 }
 612                 return
 613         }
 614         if incX == 1 {
 615                 f32.ScalUnitary(alpha, x[:n])
 616                 return
 617         }
 618         f32.ScalInc(alpha, x, uintptr(n), uintptr(incX))
 619 }