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

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