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.
8 "gonum.org/v1/gonum/blas"
9 "gonum.org/v1/gonum/blas/gonum"
12 var blas32 blas.Float32 = gonum.Implementation{}
14 // Use sets the BLAS float32 implementation to be used by subsequent BLAS calls.
15 // The default implementation is native.Implementation.
16 func Use(b blas.Float32) {
20 // Implementation returns the current BLAS float32 implementation.
22 // Implementation allows direct calls to the current the BLAS float32 implementation
23 // giving finer control of parameters.
24 func Implementation() blas.Float32 {
28 // Vector represents a vector with an associated element increment.
34 // General represents a matrix using the conventional storage scheme.
41 // Band represents a band matrix using the band storage scheme.
49 // Triangular represents a triangular matrix using the conventional storage scheme.
50 type Triangular struct {
58 // TriangularBand represents a triangular matrix using the band storage scheme.
59 type TriangularBand struct {
67 // TriangularPacked represents a triangular matrix using the packed storage scheme.
68 type TriangularPacked struct {
75 // Symmetric represents a symmetric matrix using the conventional storage scheme.
76 type Symmetric struct {
83 // SymmetricBand represents a symmetric matrix using the band storage scheme.
84 type SymmetricBand struct {
91 // SymmetricPacked represents a symmetric matrix using the packed storage scheme.
92 type SymmetricPacked struct {
100 const negInc = "blas32: negative vector increment"
102 // Dot computes the dot product of the two vectors:
104 func Dot(n int, x, y Vector) float32 {
105 return blas32.Sdot(n, x.Data, x.Inc, y.Data, y.Inc)
108 // DDot computes the dot product of the two vectors:
110 func DDot(n int, x, y Vector) float64 {
111 return blas32.Dsdot(n, x.Data, x.Inc, y.Data, y.Inc)
114 // SDDot computes the dot product of the two vectors adding a constant:
115 // alpha + \sum_i x[i]*y[i].
116 func SDDot(n int, alpha float32, x, y Vector) float32 {
117 return blas32.Sdsdot(n, alpha, x.Data, x.Inc, y.Data, y.Inc)
120 // Nrm2 computes the Euclidean norm of the vector x:
121 // sqrt(\sum_i x[i]*x[i]).
123 // Nrm2 will panic if the vector increment is negative.
124 func Nrm2(n int, x Vector) float32 {
128 return blas32.Snrm2(n, x.Data, x.Inc)
131 // Asum computes the sum of the absolute values of the elements of x:
134 // Asum will panic if the vector increment is negative.
135 func Asum(n int, x Vector) float32 {
139 return blas32.Sasum(n, x.Data, x.Inc)
142 // Iamax returns the index of an element of x with the largest absolute value.
143 // If there are multiple such indices the earliest is returned.
144 // Iamax returns -1 if n == 0.
146 // Iamax will panic if the vector increment is negative.
147 func Iamax(n int, x Vector) int {
151 return blas32.Isamax(n, x.Data, x.Inc)
154 // Swap exchanges the elements of the two vectors:
155 // x[i], y[i] = y[i], x[i] for all i.
156 func Swap(n int, x, y Vector) {
157 blas32.Sswap(n, x.Data, x.Inc, y.Data, y.Inc)
160 // Copy copies the elements of x into the elements of y:
161 // y[i] = x[i] for all i.
162 func Copy(n int, x, y Vector) {
163 blas32.Scopy(n, x.Data, x.Inc, y.Data, y.Inc)
166 // Axpy adds x scaled by alpha to y:
167 // y[i] += alpha*x[i] for all i.
168 func Axpy(n int, alpha float32, x, y Vector) {
169 blas32.Saxpy(n, alpha, x.Data, x.Inc, y.Data, y.Inc)
172 // Rotg computes the parameters of a Givens plane rotation so that
174 // ⎣-s c⎦ * ⎣b⎦ = ⎣0⎦
175 // where a and b are the Cartesian coordinates of a given point.
176 // c, s, and r are defined as
177 // r = ±Sqrt(a^2 + b^2),
178 // c = a/r, the cosine of the rotation angle,
179 // s = a/r, the sine of the rotation angle,
180 // and z is defined such that
181 // if |a| > |b|, z = s,
182 // otherwise if c != 0, z = 1/c,
184 func Rotg(a, b float32) (c, s, r, z float32) {
185 return blas32.Srotg(a, b)
188 // Rotmg computes the modified Givens rotation. See
189 // http://www.netlib.org/lapack/explore-html/df/deb/drotmg_8f.html
191 func Rotmg(d1, d2, b1, b2 float32) (p blas.SrotmParams, rd1, rd2, rb1 float32) {
192 return blas32.Srotmg(d1, d2, b1, b2)
195 // Rot applies a plane transformation to n points represented by the vectors x
197 // x[i] = c*x[i] + s*y[i],
198 // y[i] = -s*x[i] + c*y[i], for all i.
199 func Rot(n int, x, y Vector, c, s float32) {
200 blas32.Srot(n, x.Data, x.Inc, y.Data, y.Inc, c, s)
203 // Rotm applies the modified Givens rotation to n points represented by the
205 func Rotm(n int, x, y Vector, p blas.SrotmParams) {
206 blas32.Srotm(n, x.Data, x.Inc, y.Data, y.Inc, p)
209 // Scal scales the vector x by alpha:
210 // x[i] *= alpha for all i.
212 // Scal will panic if the vector increment is negative.
213 func Scal(n int, alpha float32, x Vector) {
217 blas32.Sscal(n, alpha, x.Data, x.Inc)
223 // y = alpha * A * x + beta * y, if t == blas.NoTrans,
224 // y = alpha * A^T * x + beta * y, if t == blas.Trans or blas.ConjTrans,
225 // where A is an m×n dense matrix, x and y are vectors, and alpha and beta are scalars.
226 func Gemv(t blas.Transpose, alpha float32, a General, x Vector, beta float32, y Vector) {
227 blas32.Sgemv(t, a.Rows, a.Cols, alpha, a.Data, a.Stride, x.Data, x.Inc, beta, y.Data, y.Inc)
231 // y = alpha * A * x + beta * y, if t == blas.NoTrans,
232 // y = alpha * A^T * x + beta * y, if t == blas.Trans or blas.ConjTrans,
233 // where A is an m×n band matrix, x and y are vectors, and alpha and beta are scalars.
234 func Gbmv(t blas.Transpose, alpha float32, a Band, x Vector, beta float32, y Vector) {
235 blas32.Sgbmv(t, a.Rows, a.Cols, a.KL, a.KU, alpha, a.Data, a.Stride, x.Data, x.Inc, beta, y.Data, y.Inc)
239 // x = A * x, if t == blas.NoTrans,
240 // x = A^T * x, if t == blas.Trans or blas.ConjTrans,
241 // where A is an n×n triangular matrix, and x is a vector.
242 func Trmv(t blas.Transpose, a Triangular, x Vector) {
243 blas32.Strmv(a.Uplo, t, a.Diag, a.N, a.Data, a.Stride, x.Data, x.Inc)
247 // x = A * x, if t == blas.NoTrans,
248 // x = A^T * x, if t == blas.Trans or blas.ConjTrans,
249 // where A is an n×n triangular band matrix, and x is a vector.
250 func Tbmv(t blas.Transpose, a TriangularBand, x Vector) {
251 blas32.Stbmv(a.Uplo, t, a.Diag, a.N, a.K, a.Data, a.Stride, x.Data, x.Inc)
255 // x = A * x, if t == blas.NoTrans,
256 // x = A^T * x, if t == blas.Trans or blas.ConjTrans,
257 // where A is an n×n triangular matrix in packed format, and x is a vector.
258 func Tpmv(t blas.Transpose, a TriangularPacked, x Vector) {
259 blas32.Stpmv(a.Uplo, t, a.Diag, a.N, a.Data, x.Data, x.Inc)
263 // A * x = b, if t == blas.NoTrans,
264 // A^T * x = b, if t == blas.Trans or blas.ConjTrans,
265 // where A is an n×n triangular matrix, and x and b are vectors.
267 // At entry to the function, x contains the values of b, and the result is
268 // stored in-place into x.
270 // No test for singularity or near-singularity is included in this
271 // routine. Such tests must be performed before calling this routine.
272 func Trsv(t blas.Transpose, a Triangular, x Vector) {
273 blas32.Strsv(a.Uplo, t, a.Diag, a.N, a.Data, a.Stride, x.Data, x.Inc)
277 // A * x = b, if t == blas.NoTrans,
278 // A^T * x = b, if t == blas.Trans or blas.ConjTrans,
279 // where A is an n×n triangular band matrix, and x and b are vectors.
281 // At entry to the function, x contains the values of b, and the result is
282 // stored in place into x.
284 // No test for singularity or near-singularity is included in this
285 // routine. Such tests must be performed before calling this routine.
286 func Tbsv(t blas.Transpose, a TriangularBand, x Vector) {
287 blas32.Stbsv(a.Uplo, t, a.Diag, a.N, a.K, a.Data, a.Stride, x.Data, x.Inc)
291 // A * x = b, if t == blas.NoTrans,
292 // A^T * x = b, if t == blas.Trans or blas.ConjTrans,
293 // where A is an n×n triangular matrix in packed format, and x and b are
296 // At entry to the function, x contains the values of b, and the result is
297 // stored in place into x.
299 // No test for singularity or near-singularity is included in this
300 // routine. Such tests must be performed before calling this routine.
301 func Tpsv(t blas.Transpose, a TriangularPacked, x Vector) {
302 blas32.Stpsv(a.Uplo, t, a.Diag, a.N, a.Data, x.Data, x.Inc)
306 // y = alpha * A * x + beta * y,
307 // where A is an n×n symmetric matrix, x and y are vectors, and alpha and
309 func Symv(alpha float32, a Symmetric, x Vector, beta float32, y Vector) {
310 blas32.Ssymv(a.Uplo, a.N, alpha, a.Data, a.Stride, x.Data, x.Inc, beta, y.Data, y.Inc)
314 // y = alpha * A * x + beta * y,
315 // where A is an n×n symmetric band matrix, x and y are vectors, and alpha
316 // and beta are scalars.
317 func Sbmv(alpha float32, a SymmetricBand, x Vector, beta float32, y Vector) {
318 blas32.Ssbmv(a.Uplo, a.N, a.K, alpha, a.Data, a.Stride, x.Data, x.Inc, beta, y.Data, y.Inc)
322 // y = alpha * A * x + beta * y,
323 // where A is an n×n symmetric matrix in packed format, x and y are vectors,
324 // and alpha and beta are scalars.
325 func Spmv(alpha float32, a SymmetricPacked, x Vector, beta float32, y Vector) {
326 blas32.Sspmv(a.Uplo, a.N, alpha, a.Data, x.Data, x.Inc, beta, y.Data, y.Inc)
329 // Ger performs a rank-1 update
330 // A += alpha * x * y^T,
331 // where A is an m×n dense matrix, x and y are vectors, and alpha is a scalar.
332 func Ger(alpha float32, x, y Vector, a General) {
333 blas32.Sger(a.Rows, a.Cols, alpha, x.Data, x.Inc, y.Data, y.Inc, a.Data, a.Stride)
336 // Syr performs a rank-1 update
337 // A += alpha * x * x^T,
338 // where A is an n×n symmetric matrix, x is a vector, and alpha is a scalar.
339 func Syr(alpha float32, x Vector, a Symmetric) {
340 blas32.Ssyr(a.Uplo, a.N, alpha, x.Data, x.Inc, a.Data, a.Stride)
343 // Spr performs the rank-1 update
344 // A += alpha * x * x^T,
345 // where A is an n×n symmetric matrix in packed format, x is a vector, and
346 // alpha is a scalar.
347 func Spr(alpha float32, x Vector, a SymmetricPacked) {
348 blas32.Sspr(a.Uplo, a.N, alpha, x.Data, x.Inc, a.Data)
351 // Syr2 performs a rank-2 update
352 // A += alpha * x * y^T + alpha * y * x^T,
353 // where A is a symmetric n×n matrix, x and y are vectors, and alpha is a scalar.
354 func Syr2(alpha float32, x, y Vector, a Symmetric) {
355 blas32.Ssyr2(a.Uplo, a.N, alpha, x.Data, x.Inc, y.Data, y.Inc, a.Data, a.Stride)
358 // Spr2 performs a rank-2 update
359 // A += alpha * x * y^T + alpha * y * x^T,
360 // where A is an n×n symmetric matrix in packed format, x and y are vectors,
361 // and alpha is a scalar.
362 func Spr2(alpha float32, x, y Vector, a SymmetricPacked) {
363 blas32.Sspr2(a.Uplo, a.N, alpha, x.Data, x.Inc, y.Data, y.Inc, a.Data)
369 // C = alpha * A * B + beta * C,
370 // where A, B, and C are dense matrices, and alpha and beta are scalars.
371 // tA and tB specify whether A or B are transposed.
372 func Gemm(tA, tB blas.Transpose, alpha float32, a, b General, beta float32, c General) {
374 if tA == blas.NoTrans {
375 m, k = a.Rows, a.Cols
377 m, k = a.Cols, a.Rows
379 if tB == blas.NoTrans {
384 blas32.Sgemm(tA, tB, m, n, k, alpha, a.Data, a.Stride, b.Data, b.Stride, beta, c.Data, c.Stride)
388 // C = alpha * A * B + beta * C, if s == blas.Left,
389 // C = alpha * B * A + beta * C, if s == blas.Right,
390 // where A is an n×n or m×m symmetric matrix, B and C are m×n matrices, and
391 // alpha is a scalar.
392 func Symm(s blas.Side, alpha float32, a Symmetric, b General, beta float32, c General) {
399 blas32.Ssymm(s, a.Uplo, m, n, alpha, a.Data, a.Stride, b.Data, b.Stride, beta, c.Data, c.Stride)
402 // Syrk performs a symmetric rank-k update
403 // C = alpha * A * A^T + beta * C, if t == blas.NoTrans,
404 // C = alpha * A^T * A + beta * C, if t == blas.Trans or blas.ConjTrans,
405 // where C is an n×n symmetric matrix, A is an n×k matrix if t == blas.NoTrans and
406 // a k×n matrix otherwise, and alpha and beta are scalars.
407 func Syrk(t blas.Transpose, alpha float32, a General, beta float32, c Symmetric) {
409 if t == blas.NoTrans {
410 n, k = a.Rows, a.Cols
412 n, k = a.Cols, a.Rows
414 blas32.Ssyrk(c.Uplo, t, n, k, alpha, a.Data, a.Stride, beta, c.Data, c.Stride)
417 // Syr2k performs a symmetric rank-2k update
418 // C = alpha * A * B^T + alpha * B * A^T + beta * C, if t == blas.NoTrans,
419 // C = alpha * A^T * B + alpha * B^T * A + beta * C, if t == blas.Trans or blas.ConjTrans,
420 // where C is an n×n symmetric matrix, A and B are n×k matrices if t == NoTrans
421 // and k×n matrices otherwise, and alpha and beta are scalars.
422 func Syr2k(t blas.Transpose, alpha float32, a, b General, beta float32, c Symmetric) {
424 if t == blas.NoTrans {
425 n, k = a.Rows, a.Cols
427 n, k = a.Cols, a.Rows
429 blas32.Ssyr2k(c.Uplo, t, n, k, alpha, a.Data, a.Stride, b.Data, b.Stride, beta, c.Data, c.Stride)
433 // B = alpha * A * B, if tA == blas.NoTrans and s == blas.Left,
434 // B = alpha * A^T * B, if tA == blas.Trans or blas.ConjTrans, and s == blas.Left,
435 // B = alpha * B * A, if tA == blas.NoTrans and s == blas.Right,
436 // B = alpha * B * A^T, if tA == blas.Trans or blas.ConjTrans, and s == blas.Right,
437 // where A is an n×n or m×m triangular matrix, B is an m×n matrix, and alpha is
439 func Trmm(s blas.Side, tA blas.Transpose, alpha float32, a Triangular, b General) {
440 blas32.Strmm(s, a.Uplo, tA, a.Diag, b.Rows, b.Cols, alpha, a.Data, a.Stride, b.Data, b.Stride)
444 // A * X = alpha * B, if tA == blas.NoTrans and s == blas.Left,
445 // A^T * X = alpha * B, if tA == blas.Trans or blas.ConjTrans, and s == blas.Left,
446 // X * A = alpha * B, if tA == blas.NoTrans and s == blas.Right,
447 // X * A^T = alpha * B, if tA == blas.Trans or blas.ConjTrans, and s == blas.Right,
448 // where A is an n×n or m×m triangular matrix, X and B are m×n matrices, and
449 // alpha is a scalar.
451 // At entry to the function, X contains the values of B, and the result is
452 // stored in-place into X.
454 // No check is made that A is invertible.
455 func Trsm(s blas.Side, tA blas.Transpose, alpha float32, a Triangular, b General) {
456 blas32.Strsm(s, a.Uplo, tA, a.Diag, b.Rows, b.Cols, alpha, a.Data, a.Stride, b.Data, b.Stride)