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