blas64.go 15 KB

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