vendor/gonum.org/v1/gonum/blas/blas64/blas64.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 blas64
   6
   7 import (
   8         "gonum.org/v1/gonum/blas"
   9         "gonum.org/v1/gonum/blas/gonum"
  10 )
  11
  12 var blas64 blas.Float64 = gonum.Implementation{}
  13
  14 // Use sets the BLAS float64 implementation to be used by subsequent BLAS calls.
  15 // The default implementation is native.Implementation.
  16 func Use(b blas.Float64) {
  17         blas64 = b
  18 }
  19
  20 // Implementation returns the current BLAS float64 implementation.
  21 //
  22 // Implementation allows direct calls to the current the BLAS float64 implementation
  23 // giving finer control of parameters.
  24 func Implementation() blas.Float64 {
  25         return blas64
  26 }
  27
  28 // Vector represents a vector with an associated element increment.
  29 type Vector struct {
  30         Inc  int
  31         Data []float64
  32 }
  33
  34 // General represents a matrix using the conventional storage scheme.
  35 type General struct {
  36         Rows, Cols int
  37         Stride     int
  38         Data       []float64
  39 }
  40
  41 // Band represents a band matrix using the band storage scheme.
  42 type Band struct {
  43         Rows, Cols int
  44         KL, KU     int
  45         Stride     int
  46         Data       []float64
  47 }
  48
  49 // Triangular represents a triangular matrix using the conventional storage scheme.
  50 type Triangular struct {
  51         N      int
  52         Stride int
  53         Data   []float64
  54         Uplo   blas.Uplo
  55         Diag   blas.Diag
  56 }
  57
  58 // TriangularBand represents a triangular matrix using the band storage scheme.
  59 type TriangularBand struct {
  60         N, K   int
  61         Stride int
  62         Data   []float64
  63         Uplo   blas.Uplo
  64         Diag   blas.Diag
  65 }
  66
  67 // TriangularPacked represents a triangular matrix using the packed storage scheme.
  68 type TriangularPacked struct {
  69         N    int
  70         Data []float64
  71         Uplo blas.Uplo
  72         Diag blas.Diag
  73 }
  74
  75 // Symmetric represents a symmetric matrix using the conventional storage scheme.
  76 type Symmetric struct {
  77         N      int
  78         Stride int
  79         Data   []float64
  80         Uplo   blas.Uplo
  81 }
  82
  83 // SymmetricBand represents a symmetric matrix using the band storage scheme.
  84 type SymmetricBand struct {
  85         N, K   int
  86         Stride int
  87         Data   []float64
  88         Uplo   blas.Uplo
  89 }
  90
  91 // SymmetricPacked represents a symmetric matrix using the packed storage scheme.
  92 type SymmetricPacked struct {
  93         N    int
  94         Data []float64
  95         Uplo blas.Uplo
  96 }
  97
  98 // Level 1
  99
 100 const negInc = "blas64: negative vector increment"
 101
 102 // Dot computes the dot product of the two vectors:
 103 //  \sum_i x[i]*y[i].
 104 func Dot(n int, x, y Vector) float64 {
 105         return blas64.Ddot(n, x.Data, x.Inc, y.Data, y.Inc)
 106 }
 107
 108 // Nrm2 computes the Euclidean norm of the vector x:
 109 //  sqrt(\sum_i x[i]*x[i]).
 110 //
 111 // Nrm2 will panic if the vector increment is negative.
 112 func Nrm2(n int, x Vector) float64 {
 113         if x.Inc < 0 {
 114                 panic(negInc)
 115         }
 116         return blas64.Dnrm2(n, x.Data, x.Inc)
 117 }
 118
 119 // Asum computes the sum of the absolute values of the elements of x:
 120 //  \sum_i |x[i]|.
 121 //
 122 // Asum will panic if the vector increment is negative.
 123 func Asum(n int, x Vector) float64 {
 124         if x.Inc < 0 {
 125                 panic(negInc)
 126         }
 127         return blas64.Dasum(n, x.Data, x.Inc)
 128 }
 129
 130 // Iamax returns the index of an element of x with the largest absolute value.
 131 // If there are multiple such indices the earliest is returned.
 132 // Iamax returns -1 if n == 0.
 133 //
 134 // Iamax will panic if the vector increment is negative.
 135 func Iamax(n int, x Vector) int {
 136         if x.Inc < 0 {
 137                 panic(negInc)
 138         }
 139         return blas64.Idamax(n, x.Data, x.Inc)
 140 }
 141
 142 // Swap exchanges the elements of the two vectors:
 143 //  x[i], y[i] = y[i], x[i] for all i.
 144 func Swap(n int, x, y Vector) {
 145         blas64.Dswap(n, x.Data, x.Inc, y.Data, y.Inc)
 146 }
 147
 148 // Copy copies the elements of x into the elements of y:
 149 //  y[i] = x[i] for all i.
 150 func Copy(n int, x, y Vector) {
 151         blas64.Dcopy(n, x.Data, x.Inc, y.Data, y.Inc)
 152 }
 153
 154 // Axpy adds x scaled by alpha to y:
 155 //  y[i] += alpha*x[i] for all i.
 156 func Axpy(n int, alpha float64, x, y Vector) {
 157         blas64.Daxpy(n, alpha, x.Data, x.Inc, y.Data, y.Inc)
 158 }
 159
 160 // Rotg computes the parameters of a Givens plane rotation so that
 161 //  ⎡ c s⎤   ⎡a⎤   ⎡r⎤
 162 //  ⎣-s c⎦ * ⎣b⎦ = ⎣0⎦
 163 // where a and b are the Cartesian coordinates of a given point.
 164 // c, s, and r are defined as
 165 //  r = ±Sqrt(a^2 + b^2),
 166 //  c = a/r, the cosine of the rotation angle,
 167 //  s = a/r, the sine of the rotation angle,
 168 // and z is defined such that
 169 //  if |a| > |b|,        z = s,
 170 //  otherwise if c != 0, z = 1/c,
 171 //  otherwise            z = 1.
 172 func Rotg(a, b float64) (c, s, r, z float64) {
 173         return blas64.Drotg(a, b)
 174 }
 175
 176 // Rotmg computes the modified Givens rotation. See
 177 // http://www.netlib.org/lapack/explore-html/df/deb/drotmg_8f.html
 178 // for more details.
 179 func Rotmg(d1, d2, b1, b2 float64) (p blas.DrotmParams, rd1, rd2, rb1 float64) {
 180         return blas64.Drotmg(d1, d2, b1, b2)
 181 }
 182
 183 // Rot applies a plane transformation to n points represented by the vectors x
 184 // and y:
 185 //  x[i] =  c*x[i] + s*y[i],
 186 //  y[i] = -s*x[i] + c*y[i], for all i.
 187 func Rot(n int, x, y Vector, c, s float64) {
 188         blas64.Drot(n, x.Data, x.Inc, y.Data, y.Inc, c, s)
 189 }
 190
 191 // Rotm applies the modified Givens rotation to n points represented by the
 192 // vectors x and y.
 193 func Rotm(n int, x, y Vector, p blas.DrotmParams) {
 194         blas64.Drotm(n, x.Data, x.Inc, y.Data, y.Inc, p)
 195 }
 196
 197 // Scal scales the vector x by alpha:
 198 //  x[i] *= alpha for all i.
 199 //
 200 // Scal will panic if the vector increment is negative.
 201 func Scal(n int, alpha float64, x Vector) {
 202         if x.Inc < 0 {
 203                 panic(negInc)
 204         }
 205         blas64.Dscal(n, alpha, x.Data, x.Inc)
 206 }
 207
 208 // Level 2
 209
 210 // Gemv computes
 211 //  y = alpha * A * x + beta * y,   if t == blas.NoTrans,
 212 //  y = alpha * A^T * x + beta * y, if t == blas.Trans or blas.ConjTrans,
 213 // where A is an m×n dense matrix, x and y are vectors, and alpha and beta are scalars.
 214 func Gemv(t blas.Transpose, alpha float64, a General, x Vector, beta float64, y Vector) {
 215         blas64.Dgemv(t, a.Rows, a.Cols, alpha, a.Data, a.Stride, x.Data, x.Inc, beta, y.Data, y.Inc)
 216 }
 217
 218 // Gbmv computes
 219 //  y = alpha * A * x + beta * y,   if t == blas.NoTrans,
 220 //  y = alpha * A^T * x + beta * y, if t == blas.Trans or blas.ConjTrans,
 221 // where A is an m×n band matrix, x and y are vectors, and alpha and beta are scalars.
 222 func Gbmv(t blas.Transpose, alpha float64, a Band, x Vector, beta float64, y Vector) {
 223         blas64.Dgbmv(t, a.Rows, a.Cols, a.KL, a.KU, alpha, a.Data, a.Stride, x.Data, x.Inc, beta, y.Data, y.Inc)
 224 }
 225
 226 // Trmv computes
 227 //  x = A * x,   if t == blas.NoTrans,
 228 //  x = A^T * x, if t == blas.Trans or blas.ConjTrans,
 229 // where A is an n×n triangular matrix, and x is a vector.
 230 func Trmv(t blas.Transpose, a Triangular, x Vector) {
 231         blas64.Dtrmv(a.Uplo, t, a.Diag, a.N, a.Data, a.Stride, x.Data, x.Inc)
 232 }
 233
 234 // Tbmv computes
 235 //  x = A * x,   if t == blas.NoTrans,
 236 //  x = A^T * x, if t == blas.Trans or blas.ConjTrans,
 237 // where A is an n×n triangular band matrix, and x is a vector.
 238 func Tbmv(t blas.Transpose, a TriangularBand, x Vector) {
 239         blas64.Dtbmv(a.Uplo, t, a.Diag, a.N, a.K, a.Data, a.Stride, x.Data, x.Inc)
 240 }
 241
 242 // Tpmv computes
 243 //  x = A * x,   if t == blas.NoTrans,
 244 //  x = A^T * x, if t == blas.Trans or blas.ConjTrans,
 245 // where A is an n×n triangular matrix in packed format, and x is a vector.
 246 func Tpmv(t blas.Transpose, a TriangularPacked, x Vector) {
 247         blas64.Dtpmv(a.Uplo, t, a.Diag, a.N, a.Data, x.Data, x.Inc)
 248 }
 249
 250 // Trsv solves
 251 //  A * x = b,   if t == blas.NoTrans,
 252 //  A^T * x = b, if t == blas.Trans or blas.ConjTrans,
 253 // where A is an n×n triangular matrix, and x and b are vectors.
 254 //
 255 // At entry to the function, x contains the values of b, and the result is
 256 // stored in-place into x.
 257 //
 258 // No test for singularity or near-singularity is included in this
 259 // routine. Such tests must be performed before calling this routine.
 260 func Trsv(t blas.Transpose, a Triangular, x Vector) {
 261         blas64.Dtrsv(a.Uplo, t, a.Diag, a.N, a.Data, a.Stride, x.Data, x.Inc)
 262 }
 263
 264 // Tbsv solves
 265 //  A * x = b,   if t == blas.NoTrans,
 266 //  A^T * x = b, if t == blas.Trans or blas.ConjTrans,
 267 // where A is an n×n triangular band matrix, and x and b are vectors.
 268 //
 269 // At entry to the function, x contains the values of b, and the result is
 270 // stored in place into x.
 271 //
 272 // No test for singularity or near-singularity is included in this
 273 // routine. Such tests must be performed before calling this routine.
 274 func Tbsv(t blas.Transpose, a TriangularBand, x Vector) {
 275         blas64.Dtbsv(a.Uplo, t, a.Diag, a.N, a.K, a.Data, a.Stride, x.Data, x.Inc)
 276 }
 277
 278 // Tpsv solves
 279 //  A * x = b,   if t == blas.NoTrans,
 280 //  A^T * x = b, if t == blas.Trans or blas.ConjTrans,
 281 // where A is an n×n triangular matrix in packed format, and x and b are
 282 // vectors.
 283 //
 284 // At entry to the function, x contains the values of b, and the result is
 285 // stored in place into x.
 286 //
 287 // No test for singularity or near-singularity is included in this
 288 // routine. Such tests must be performed before calling this routine.
 289 func Tpsv(t blas.Transpose, a TriangularPacked, x Vector) {
 290         blas64.Dtpsv(a.Uplo, t, a.Diag, a.N, a.Data, x.Data, x.Inc)
 291 }
 292
 293 // Symv computes
 294 //    y = alpha * A * x + beta * y,
 295 // where A is an n×n symmetric matrix, x and y are vectors, and alpha and
 296 // beta are scalars.
 297 func Symv(alpha float64, a Symmetric, x Vector, beta float64, y Vector) {
 298         blas64.Dsymv(a.Uplo, a.N, alpha, a.Data, a.Stride, x.Data, x.Inc, beta, y.Data, y.Inc)
 299 }
 300
 301 // Sbmv performs
 302 //  y = alpha * A * x + beta * y,
 303 // where A is an n×n symmetric band matrix, x and y are vectors, and alpha
 304 // and beta are scalars.
 305 func Sbmv(alpha float64, a SymmetricBand, x Vector, beta float64, y Vector) {
 306         blas64.Dsbmv(a.Uplo, a.N, a.K, alpha, a.Data, a.Stride, x.Data, x.Inc, beta, y.Data, y.Inc)
 307 }
 308
 309 // Spmv performs
 310 //    y = alpha * A * x + beta * y,
 311 // where A is an n×n symmetric matrix in packed format, x and y are vectors,
 312 // and alpha and beta are scalars.
 313 func Spmv(alpha float64, a SymmetricPacked, x Vector, beta float64, y Vector) {
 314         blas64.Dspmv(a.Uplo, a.N, alpha, a.Data, x.Data, x.Inc, beta, y.Data, y.Inc)
 315 }
 316
 317 // Ger performs a rank-1 update
 318 //  A += alpha * x * y^T,
 319 // where A is an m×n dense matrix, x and y are vectors, and alpha is a scalar.
 320 func Ger(alpha float64, x, y Vector, a General) {
 321         blas64.Dger(a.Rows, a.Cols, alpha, x.Data, x.Inc, y.Data, y.Inc, a.Data, a.Stride)
 322 }
 323
 324 // Syr performs a rank-1 update
 325 //  A += alpha * x * x^T,
 326 // where A is an n×n symmetric matrix, x is a vector, and alpha is a scalar.
 327 func Syr(alpha float64, x Vector, a Symmetric) {
 328         blas64.Dsyr(a.Uplo, a.N, alpha, x.Data, x.Inc, a.Data, a.Stride)
 329 }
 330
 331 // Spr performs the rank-1 update
 332 //  A += alpha * x * x^T,
 333 // where A is an n×n symmetric matrix in packed format, x is a vector, and
 334 // alpha is a scalar.
 335 func Spr(alpha float64, x Vector, a SymmetricPacked) {
 336         blas64.Dspr(a.Uplo, a.N, alpha, x.Data, x.Inc, a.Data)
 337 }
 338
 339 // Syr2 performs a rank-2 update
 340 //  A += alpha * x * y^T + alpha * y * x^T,
 341 // where A is a symmetric n×n matrix, x and y are vectors, and alpha is a scalar.
 342 func Syr2(alpha float64, x, y Vector, a Symmetric) {
 343         blas64.Dsyr2(a.Uplo, a.N, alpha, x.Data, x.Inc, y.Data, y.Inc, a.Data, a.Stride)
 344 }
 345
 346 // Spr2 performs a rank-2 update
 347 //  A += alpha * x * y^T + alpha * y * x^T,
 348 // where A is an n×n symmetric matrix in packed format, x and y are vectors,
 349 // and alpha is a scalar.
 350 func Spr2(alpha float64, x, y Vector, a SymmetricPacked) {
 351         blas64.Dspr2(a.Uplo, a.N, alpha, x.Data, x.Inc, y.Data, y.Inc, a.Data)
 352 }
 353
 354 // Level 3
 355
 356 // Gemm computes
 357 //  C = alpha * A * B + beta * C,
 358 // where A, B, and C are dense matrices, and alpha and beta are scalars.
 359 // tA and tB specify whether A or B are transposed.
 360 func Gemm(tA, tB blas.Transpose, alpha float64, a, b General, beta float64, c General) {
 361         var m, n, k int
 362         if tA == blas.NoTrans {
 363                 m, k = a.Rows, a.Cols
 364         } else {
 365                 m, k = a.Cols, a.Rows
 366         }
 367         if tB == blas.NoTrans {
 368                 n = b.Cols
 369         } else {
 370                 n = b.Rows
 371         }
 372         blas64.Dgemm(tA, tB, m, n, k, alpha, a.Data, a.Stride, b.Data, b.Stride, beta, c.Data, c.Stride)
 373 }
 374
 375 // Symm performs
 376 //  C = alpha * A * B + beta * C, if s == blas.Left,
 377 //  C = alpha * B * A + beta * C, if s == blas.Right,
 378 // where A is an n×n or m×m symmetric matrix, B and C are m×n matrices, and
 379 // alpha is a scalar.
 380 func Symm(s blas.Side, alpha float64, a Symmetric, b General, beta float64, c General) {
 381         var m, n int
 382         if s == blas.Left {
 383                 m, n = a.N, b.Cols
 384         } else {
 385                 m, n = b.Rows, a.N
 386         }
 387         blas64.Dsymm(s, a.Uplo, m, n, alpha, a.Data, a.Stride, b.Data, b.Stride, beta, c.Data, c.Stride)
 388 }
 389
 390 // Syrk performs a symmetric rank-k update
 391 //  C = alpha * A * A^T + beta * C, if t == blas.NoTrans,
 392 //  C = alpha * A^T * A + beta * C, if t == blas.Trans or blas.ConjTrans,
 393 // where C is an n×n symmetric matrix, A is an n×k matrix if t == blas.NoTrans and
 394 // a k×n matrix otherwise, and alpha and beta are scalars.
 395 func Syrk(t blas.Transpose, alpha float64, a General, beta float64, c Symmetric) {
 396         var n, k int
 397         if t == blas.NoTrans {
 398                 n, k = a.Rows, a.Cols
 399         } else {
 400                 n, k = a.Cols, a.Rows
 401         }
 402         blas64.Dsyrk(c.Uplo, t, n, k, alpha, a.Data, a.Stride, beta, c.Data, c.Stride)
 403 }
 404
 405 // Syr2k performs a symmetric rank-2k update
 406 //  C = alpha * A * B^T + alpha * B * A^T + beta * C, if t == blas.NoTrans,
 407 //  C = alpha * A^T * B + alpha * B^T * A + beta * C, if t == blas.Trans or blas.ConjTrans,
 408 // where C is an n×n symmetric matrix, A and B are n×k matrices if t == NoTrans
 409 // and k×n matrices otherwise, and alpha and beta are scalars.
 410 func Syr2k(t blas.Transpose, alpha float64, a, b General, beta float64, c Symmetric) {
 411         var n, k int
 412         if t == blas.NoTrans {
 413                 n, k = a.Rows, a.Cols
 414         } else {
 415                 n, k = a.Cols, a.Rows
 416         }
 417         blas64.Dsyr2k(c.Uplo, t, n, k, alpha, a.Data, a.Stride, b.Data, b.Stride, beta, c.Data, c.Stride)
 418 }
 419
 420 // Trmm performs
 421 //  B = alpha * A * B,   if tA == blas.NoTrans and s == blas.Left,
 422 //  B = alpha * A^T * B, if tA == blas.Trans or blas.ConjTrans, and s == blas.Left,
 423 //  B = alpha * B * A,   if tA == blas.NoTrans and s == blas.Right,
 424 //  B = alpha * B * A^T, if tA == blas.Trans or blas.ConjTrans, and s == blas.Right,
 425 // where A is an n×n or m×m triangular matrix, B is an m×n matrix, and alpha is
 426 // a scalar.
 427 func Trmm(s blas.Side, tA blas.Transpose, alpha float64, a Triangular, b General) {
 428         blas64.Dtrmm(s, a.Uplo, tA, a.Diag, b.Rows, b.Cols, alpha, a.Data, a.Stride, b.Data, b.Stride)
 429 }
 430
 431 // Trsm solves
 432 //  A * X = alpha * B,   if tA == blas.NoTrans and s == blas.Left,
 433 //  A^T * X = alpha * B, if tA == blas.Trans or blas.ConjTrans, and s == blas.Left,
 434 //  X * A = alpha * B,   if tA == blas.NoTrans and s == blas.Right,
 435 //  X * A^T = alpha * B, if tA == blas.Trans or blas.ConjTrans, and s == blas.Right,
 436 // where A is an n×n or m×m triangular matrix, X and B are m×n matrices, and
 437 // alpha is a scalar.
 438 //
 439 // At entry to the function, X contains the values of B, and the result is
 440 // stored in-place into X.
 441 //
 442 // No check is made that A is invertible.
 443 func Trsm(s blas.Side, tA blas.Transpose, alpha float64, a Triangular, b General) {
 444         blas64.Dtrsm(s, a.Uplo, tA, a.Diag, b.Rows, b.Cols, alpha, a.Data, a.Stride, b.Data, b.Stride)
 445 }