| // Code generated by "go generate gonum.org/v1/gonum/blas/gonum”; DO NOT EDIT. |
| |
| // Copyright ©2015 The gonum Authors. All rights reserved. |
| // Use of this source code is governed by a BSD-style |
| // license that can be found in the LICENSE file. |
| |
| package gonum |
| |
| import ( |
| math "gonum.org/v1/gonum/blas/gonum/internal/math32" |
| |
| "gonum.org/v1/gonum/blas" |
| "gonum.org/v1/gonum/internal/asm/f32" |
| ) |
| |
| var _ blas.Float32Level1 = Implementation{} |
| |
| // Snrm2 computes the Euclidean norm of a vector, |
| // sqrt(\sum_i x[i] * x[i]). |
| // This function returns 0 if incX is negative. |
| // |
| // Float32 implementations are autogenerated and not directly tested. |
| func (Implementation) Snrm2(n int, x []float32, incX int) float32 { |
| if incX < 1 { |
| if incX == 0 { |
| panic(zeroIncX) |
| } |
| return 0 |
| } |
| if incX > 0 && (n-1)*incX >= len(x) { |
| panic(badX) |
| } |
| if n < 2 { |
| if n == 1 { |
| return math.Abs(x[0]) |
| } |
| if n == 0 { |
| return 0 |
| } |
| if n < 1 { |
| panic(negativeN) |
| } |
| } |
| var ( |
| scale float32 = 0 |
| sumSquares float32 = 1 |
| ) |
| if incX == 1 { |
| x = x[:n] |
| for _, v := range x { |
| if v == 0 { |
| continue |
| } |
| absxi := math.Abs(v) |
| if math.IsNaN(absxi) { |
| return math.NaN() |
| } |
| if scale < absxi { |
| sumSquares = 1 + sumSquares*(scale/absxi)*(scale/absxi) |
| scale = absxi |
| } else { |
| sumSquares = sumSquares + (absxi/scale)*(absxi/scale) |
| } |
| } |
| if math.IsInf(scale, 1) { |
| return math.Inf(1) |
| } |
| return scale * math.Sqrt(sumSquares) |
| } |
| for ix := 0; ix < n*incX; ix += incX { |
| val := x[ix] |
| if val == 0 { |
| continue |
| } |
| absxi := math.Abs(val) |
| if math.IsNaN(absxi) { |
| return math.NaN() |
| } |
| if scale < absxi { |
| sumSquares = 1 + sumSquares*(scale/absxi)*(scale/absxi) |
| scale = absxi |
| } else { |
| sumSquares = sumSquares + (absxi/scale)*(absxi/scale) |
| } |
| } |
| if math.IsInf(scale, 1) { |
| return math.Inf(1) |
| } |
| return scale * math.Sqrt(sumSquares) |
| } |
| |
| // Sasum computes the sum of the absolute values of the elements of x. |
| // \sum_i |x[i]| |
| // Sasum returns 0 if incX is negative. |
| // |
| // Float32 implementations are autogenerated and not directly tested. |
| func (Implementation) Sasum(n int, x []float32, incX int) float32 { |
| var sum float32 |
| if n < 0 { |
| panic(negativeN) |
| } |
| if incX < 1 { |
| if incX == 0 { |
| panic(zeroIncX) |
| } |
| return 0 |
| } |
| if incX > 0 && (n-1)*incX >= len(x) { |
| panic(badX) |
| } |
| if incX == 1 { |
| x = x[:n] |
| for _, v := range x { |
| sum += math.Abs(v) |
| } |
| return sum |
| } |
| for i := 0; i < n; i++ { |
| sum += math.Abs(x[i*incX]) |
| } |
| return sum |
| } |
| |
| // Isamax returns the index of an element of x with the largest absolute value. |
| // If there are multiple such indices the earliest is returned. |
| // Isamax returns -1 if n == 0. |
| // |
| // Float32 implementations are autogenerated and not directly tested. |
| func (Implementation) Isamax(n int, x []float32, incX int) int { |
| if incX < 1 { |
| if incX == 0 { |
| panic(zeroIncX) |
| } |
| return -1 |
| } |
| if incX > 0 && (n-1)*incX >= len(x) { |
| panic(badX) |
| } |
| if n < 2 { |
| if n == 1 { |
| return 0 |
| } |
| if n == 0 { |
| return -1 // Netlib returns invalid index when n == 0 |
| } |
| if n < 1 { |
| panic(negativeN) |
| } |
| } |
| idx := 0 |
| max := math.Abs(x[0]) |
| if incX == 1 { |
| for i, v := range x[:n] { |
| absV := math.Abs(v) |
| if absV > max { |
| max = absV |
| idx = i |
| } |
| } |
| return idx |
| } |
| ix := incX |
| for i := 1; i < n; i++ { |
| v := x[ix] |
| absV := math.Abs(v) |
| if absV > max { |
| max = absV |
| idx = i |
| } |
| ix += incX |
| } |
| return idx |
| } |
| |
| // Sswap exchanges the elements of two vectors. |
| // x[i], y[i] = y[i], x[i] for all i |
| // |
| // Float32 implementations are autogenerated and not directly tested. |
| func (Implementation) Sswap(n int, x []float32, incX int, y []float32, incY int) { |
| if incX == 0 { |
| panic(zeroIncX) |
| } |
| if incY == 0 { |
| panic(zeroIncY) |
| } |
| if n < 1 { |
| if n == 0 { |
| return |
| } |
| panic(negativeN) |
| } |
| if (incX > 0 && (n-1)*incX >= len(x)) || (incX < 0 && (1-n)*incX >= len(x)) { |
| panic(badX) |
| } |
| if (incY > 0 && (n-1)*incY >= len(y)) || (incY < 0 && (1-n)*incY >= len(y)) { |
| panic(badY) |
| } |
| if incX == 1 && incY == 1 { |
| x = x[:n] |
| for i, v := range x { |
| x[i], y[i] = y[i], v |
| } |
| return |
| } |
| var ix, iy int |
| if incX < 0 { |
| ix = (-n + 1) * incX |
| } |
| if incY < 0 { |
| iy = (-n + 1) * incY |
| } |
| for i := 0; i < n; i++ { |
| x[ix], y[iy] = y[iy], x[ix] |
| ix += incX |
| iy += incY |
| } |
| } |
| |
| // Scopy copies the elements of x into the elements of y. |
| // y[i] = x[i] for all i |
| // |
| // Float32 implementations are autogenerated and not directly tested. |
| func (Implementation) Scopy(n int, x []float32, incX int, y []float32, incY int) { |
| if incX == 0 { |
| panic(zeroIncX) |
| } |
| if incY == 0 { |
| panic(zeroIncY) |
| } |
| if n < 1 { |
| if n == 0 { |
| return |
| } |
| panic(negativeN) |
| } |
| if (incX > 0 && (n-1)*incX >= len(x)) || (incX < 0 && (1-n)*incX >= len(x)) { |
| panic(badX) |
| } |
| if (incY > 0 && (n-1)*incY >= len(y)) || (incY < 0 && (1-n)*incY >= len(y)) { |
| panic(badY) |
| } |
| if incX == 1 && incY == 1 { |
| copy(y[:n], x[:n]) |
| return |
| } |
| var ix, iy int |
| if incX < 0 { |
| ix = (-n + 1) * incX |
| } |
| if incY < 0 { |
| iy = (-n + 1) * incY |
| } |
| for i := 0; i < n; i++ { |
| y[iy] = x[ix] |
| ix += incX |
| iy += incY |
| } |
| } |
| |
| // Saxpy adds alpha times x to y |
| // y[i] += alpha * x[i] for all i |
| // |
| // Float32 implementations are autogenerated and not directly tested. |
| func (Implementation) Saxpy(n int, alpha float32, x []float32, incX int, y []float32, incY int) { |
| if incX == 0 { |
| panic(zeroIncX) |
| } |
| if incY == 0 { |
| panic(zeroIncY) |
| } |
| if n < 1 { |
| if n == 0 { |
| return |
| } |
| panic(negativeN) |
| } |
| if (incX > 0 && (n-1)*incX >= len(x)) || (incX < 0 && (1-n)*incX >= len(x)) { |
| panic(badX) |
| } |
| if (incY > 0 && (n-1)*incY >= len(y)) || (incY < 0 && (1-n)*incY >= len(y)) { |
| panic(badY) |
| } |
| if alpha == 0 { |
| return |
| } |
| if incX == 1 && incY == 1 { |
| f32.AxpyUnitary(alpha, x[:n], y[:n]) |
| return |
| } |
| var ix, iy int |
| if incX < 0 { |
| ix = (-n + 1) * incX |
| } |
| if incY < 0 { |
| iy = (-n + 1) * incY |
| } |
| f32.AxpyInc(alpha, x, y, uintptr(n), uintptr(incX), uintptr(incY), uintptr(ix), uintptr(iy)) |
| } |
| |
| // Srotg computes the plane rotation |
| // _ _ _ _ _ _ |
| // | c s | | a | | r | |
| // | -s c | * | b | = | 0 | |
| // ‾ ‾ ‾ ‾ ‾ ‾ |
| // where |
| // r = ±√(a^2 + b^2) |
| // c = a/r, the cosine of the plane rotation |
| // s = b/r, the sine of the plane rotation |
| // |
| // NOTE: There is a discrepancy between the refence implementation and the BLAS |
| // technical manual regarding the sign for r when a or b are zero. |
| // Srotg agrees with the definition in the manual and other |
| // common BLAS implementations. |
| // |
| // Float32 implementations are autogenerated and not directly tested. |
| func (Implementation) Srotg(a, b float32) (c, s, r, z float32) { |
| if b == 0 && a == 0 { |
| return 1, 0, a, 0 |
| } |
| absA := math.Abs(a) |
| absB := math.Abs(b) |
| aGTb := absA > absB |
| r = math.Hypot(a, b) |
| if aGTb { |
| r = math.Copysign(r, a) |
| } else { |
| r = math.Copysign(r, b) |
| } |
| c = a / r |
| s = b / r |
| if aGTb { |
| z = s |
| } else if c != 0 { // r == 0 case handled above |
| z = 1 / c |
| } else { |
| z = 1 |
| } |
| return |
| } |
| |
| // Srotmg computes the modified Givens rotation. See |
| // http://www.netlib.org/lapack/explore-html/df/deb/drotmg_8f.html |
| // for more details. |
| // |
| // Float32 implementations are autogenerated and not directly tested. |
| func (Implementation) Srotmg(d1, d2, x1, y1 float32) (p blas.SrotmParams, rd1, rd2, rx1 float32) { |
| var p1, p2, q1, q2, u float32 |
| |
| const ( |
| gam = 4096.0 |
| gamsq = 16777216.0 |
| rgamsq = 5.9604645e-8 |
| ) |
| |
| if d1 < 0 { |
| p.Flag = blas.Rescaling |
| return |
| } |
| |
| p2 = d2 * y1 |
| if p2 == 0 { |
| p.Flag = blas.Identity |
| rd1 = d1 |
| rd2 = d2 |
| rx1 = x1 |
| return |
| } |
| p1 = d1 * x1 |
| q2 = p2 * y1 |
| q1 = p1 * x1 |
| |
| absQ1 := math.Abs(q1) |
| absQ2 := math.Abs(q2) |
| |
| if absQ1 < absQ2 && q2 < 0 { |
| p.Flag = blas.Rescaling |
| return |
| } |
| |
| if d1 == 0 { |
| p.Flag = blas.Diagonal |
| p.H[0] = p1 / p2 |
| p.H[3] = x1 / y1 |
| u = 1 + p.H[0]*p.H[3] |
| rd1, rd2 = d2/u, d1/u |
| rx1 = y1 / u |
| return |
| } |
| |
| // Now we know that d1 != 0, and d2 != 0. If d2 == 0, it would be caught |
| // when p2 == 0, and if d1 == 0, then it is caught above |
| |
| if absQ1 > absQ2 { |
| p.H[1] = -y1 / x1 |
| p.H[2] = p2 / p1 |
| u = 1 - p.H[2]*p.H[1] |
| rd1 = d1 |
| rd2 = d2 |
| rx1 = x1 |
| p.Flag = blas.OffDiagonal |
| // u must be greater than zero because |q1| > |q2|, so check from netlib |
| // is unnecessary |
| // This is left in for ease of comparison with complex routines |
| //if u > 0 { |
| rd1 /= u |
| rd2 /= u |
| rx1 *= u |
| //} |
| } else { |
| p.Flag = blas.Diagonal |
| p.H[0] = p1 / p2 |
| p.H[3] = x1 / y1 |
| u = 1 + p.H[0]*p.H[3] |
| rd1 = d2 / u |
| rd2 = d1 / u |
| rx1 = y1 * u |
| } |
| |
| for rd1 <= rgamsq || rd1 >= gamsq { |
| if p.Flag == blas.OffDiagonal { |
| p.H[0] = 1 |
| p.H[3] = 1 |
| p.Flag = blas.Rescaling |
| } else if p.Flag == blas.Diagonal { |
| p.H[1] = -1 |
| p.H[2] = 1 |
| p.Flag = blas.Rescaling |
| } |
| if rd1 <= rgamsq { |
| rd1 *= gam * gam |
| rx1 /= gam |
| p.H[0] /= gam |
| p.H[2] /= gam |
| } else { |
| rd1 /= gam * gam |
| rx1 *= gam |
| p.H[0] *= gam |
| p.H[2] *= gam |
| } |
| } |
| |
| for math.Abs(rd2) <= rgamsq || math.Abs(rd2) >= gamsq { |
| if p.Flag == blas.OffDiagonal { |
| p.H[0] = 1 |
| p.H[3] = 1 |
| p.Flag = blas.Rescaling |
| } else if p.Flag == blas.Diagonal { |
| p.H[1] = -1 |
| p.H[2] = 1 |
| p.Flag = blas.Rescaling |
| } |
| if math.Abs(rd2) <= rgamsq { |
| rd2 *= gam * gam |
| p.H[1] /= gam |
| p.H[3] /= gam |
| } else { |
| rd2 /= gam * gam |
| p.H[1] *= gam |
| p.H[3] *= gam |
| } |
| } |
| return |
| } |
| |
| // Srot applies a plane transformation. |
| // x[i] = c * x[i] + s * y[i] |
| // y[i] = c * y[i] - s * x[i] |
| // |
| // Float32 implementations are autogenerated and not directly tested. |
| func (Implementation) Srot(n int, x []float32, incX int, y []float32, incY int, c float32, s float32) { |
| if incX == 0 { |
| panic(zeroIncX) |
| } |
| if incY == 0 { |
| panic(zeroIncY) |
| } |
| if n < 1 { |
| if n == 0 { |
| return |
| } |
| panic(negativeN) |
| } |
| if (incX > 0 && (n-1)*incX >= len(x)) || (incX < 0 && (1-n)*incX >= len(x)) { |
| panic(badX) |
| } |
| if (incY > 0 && (n-1)*incY >= len(y)) || (incY < 0 && (1-n)*incY >= len(y)) { |
| panic(badY) |
| } |
| if incX == 1 && incY == 1 { |
| x = x[:n] |
| for i, vx := range x { |
| vy := y[i] |
| x[i], y[i] = c*vx+s*vy, c*vy-s*vx |
| } |
| return |
| } |
| var ix, iy int |
| if incX < 0 { |
| ix = (-n + 1) * incX |
| } |
| if incY < 0 { |
| iy = (-n + 1) * incY |
| } |
| for i := 0; i < n; i++ { |
| vx := x[ix] |
| vy := y[iy] |
| x[ix], y[iy] = c*vx+s*vy, c*vy-s*vx |
| ix += incX |
| iy += incY |
| } |
| } |
| |
| // Srotm applies the modified Givens rotation to the 2×n matrix. |
| // |
| // Float32 implementations are autogenerated and not directly tested. |
| func (Implementation) Srotm(n int, x []float32, incX int, y []float32, incY int, p blas.SrotmParams) { |
| if incX == 0 { |
| panic(zeroIncX) |
| } |
| if incY == 0 { |
| panic(zeroIncY) |
| } |
| if n <= 0 { |
| if n == 0 { |
| return |
| } |
| panic(negativeN) |
| } |
| if (incX > 0 && (n-1)*incX >= len(x)) || (incX < 0 && (1-n)*incX >= len(x)) { |
| panic(badX) |
| } |
| if (incY > 0 && (n-1)*incY >= len(y)) || (incY < 0 && (1-n)*incY >= len(y)) { |
| panic(badY) |
| } |
| |
| var h11, h12, h21, h22 float32 |
| var ix, iy int |
| switch p.Flag { |
| case blas.Identity: |
| return |
| case blas.Rescaling: |
| h11 = p.H[0] |
| h12 = p.H[2] |
| h21 = p.H[1] |
| h22 = p.H[3] |
| case blas.OffDiagonal: |
| h11 = 1 |
| h12 = p.H[2] |
| h21 = p.H[1] |
| h22 = 1 |
| case blas.Diagonal: |
| h11 = p.H[0] |
| h12 = 1 |
| h21 = -1 |
| h22 = p.H[3] |
| } |
| if incX < 0 { |
| ix = (-n + 1) * incX |
| } |
| if incY < 0 { |
| iy = (-n + 1) * incY |
| } |
| if incX == 1 && incY == 1 { |
| x = x[:n] |
| for i, vx := range x { |
| vy := y[i] |
| x[i], y[i] = vx*h11+vy*h12, vx*h21+vy*h22 |
| } |
| return |
| } |
| for i := 0; i < n; i++ { |
| vx := x[ix] |
| vy := y[iy] |
| x[ix], y[iy] = vx*h11+vy*h12, vx*h21+vy*h22 |
| ix += incX |
| iy += incY |
| } |
| } |
| |
| // Sscal scales x by alpha. |
| // x[i] *= alpha |
| // Sscal has no effect if incX < 0. |
| // |
| // Float32 implementations are autogenerated and not directly tested. |
| func (Implementation) Sscal(n int, alpha float32, x []float32, incX int) { |
| if incX < 1 { |
| if incX == 0 { |
| panic(zeroIncX) |
| } |
| return |
| } |
| if (n-1)*incX >= len(x) { |
| panic(badX) |
| } |
| if n < 1 { |
| if n == 0 { |
| return |
| } |
| panic(negativeN) |
| } |
| if alpha == 0 { |
| if incX == 1 { |
| x = x[:n] |
| for i := range x { |
| x[i] = 0 |
| } |
| return |
| } |
| for ix := 0; ix < n*incX; ix += incX { |
| x[ix] = 0 |
| } |
| return |
| } |
| if incX == 1 { |
| f32.ScalUnitary(alpha, x[:n]) |
| return |
| } |
| for ix := 0; ix < n*incX; ix += incX { |
| x[ix] *= alpha |
| } |
| } |