cblas128.go 17 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509
  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 cblas128
  5. import (
  6. "gonum.org/v1/gonum/blas"
  7. "gonum.org/v1/gonum/blas/gonum"
  8. )
  9. var cblas128 blas.Complex128 = gonum.Implementation{}
  10. // Use sets the BLAS complex128 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.Complex128) {
  14. cblas128 = b
  15. }
  16. // Implementation returns the current BLAS complex128 implementation.
  17. //
  18. // Implementation allows direct calls to the current the BLAS complex128 implementation
  19. // giving finer control of parameters.
  20. func Implementation() blas.Complex128 {
  21. return cblas128
  22. }
  23. // Vector represents a vector with an associated element increment.
  24. type Vector struct {
  25. Inc int
  26. Data []complex128
  27. }
  28. // General represents a matrix using the conventional storage scheme.
  29. type General struct {
  30. Rows, Cols int
  31. Stride int
  32. Data []complex128
  33. }
  34. // Band represents a band matrix using the band storage scheme.
  35. type Band struct {
  36. Rows, Cols int
  37. KL, KU int
  38. Stride int
  39. Data []complex128
  40. }
  41. // Triangular represents a triangular matrix using the conventional storage scheme.
  42. type Triangular struct {
  43. N int
  44. Stride int
  45. Data []complex128
  46. Uplo blas.Uplo
  47. Diag blas.Diag
  48. }
  49. // TriangularBand represents a triangular matrix using the band storage scheme.
  50. type TriangularBand struct {
  51. N, K int
  52. Stride int
  53. Data []complex128
  54. Uplo blas.Uplo
  55. Diag blas.Diag
  56. }
  57. // TriangularPacked represents a triangular matrix using the packed storage scheme.
  58. type TriangularPacked struct {
  59. N int
  60. Data []complex128
  61. Uplo blas.Uplo
  62. Diag blas.Diag
  63. }
  64. // Symmetric represents a symmetric matrix using the conventional storage scheme.
  65. type Symmetric struct {
  66. N int
  67. Stride int
  68. Data []complex128
  69. Uplo blas.Uplo
  70. }
  71. // SymmetricBand represents a symmetric matrix using the band storage scheme.
  72. type SymmetricBand struct {
  73. N, K int
  74. Stride int
  75. Data []complex128
  76. Uplo blas.Uplo
  77. }
  78. // SymmetricPacked represents a symmetric matrix using the packed storage scheme.
  79. type SymmetricPacked struct {
  80. N int
  81. Data []complex128
  82. Uplo blas.Uplo
  83. }
  84. // Hermitian represents an Hermitian matrix using the conventional storage scheme.
  85. type Hermitian Symmetric
  86. // HermitianBand represents an Hermitian matrix using the band storage scheme.
  87. type HermitianBand SymmetricBand
  88. // HermitianPacked represents an Hermitian matrix using the packed storage scheme.
  89. type HermitianPacked SymmetricPacked
  90. // Level 1
  91. const negInc = "cblas128: negative vector increment"
  92. // Dotu computes the dot product of the two vectors without
  93. // complex conjugation:
  94. // x^T * y.
  95. func Dotu(n int, x, y Vector) complex128 {
  96. return cblas128.Zdotu(n, x.Data, x.Inc, y.Data, y.Inc)
  97. }
  98. // Dotc computes the dot product of the two vectors with
  99. // complex conjugation:
  100. // x^H * y.
  101. func Dotc(n int, x, y Vector) complex128 {
  102. return cblas128.Zdotc(n, x.Data, x.Inc, y.Data, y.Inc)
  103. }
  104. // Nrm2 computes the Euclidean norm of the vector x:
  105. // sqrt(\sum_i x[i] * x[i]).
  106. //
  107. // Nrm2 will panic if the vector increment is negative.
  108. func Nrm2(n int, x Vector) float64 {
  109. if x.Inc < 0 {
  110. panic(negInc)
  111. }
  112. return cblas128.Dznrm2(n, x.Data, x.Inc)
  113. }
  114. // Asum computes the sum of magnitudes of the real and imaginary parts of
  115. // elements of the vector x:
  116. // \sum_i (|Re x[i]| + |Im x[i]|).
  117. //
  118. // Asum will panic if the vector increment is negative.
  119. func Asum(n int, x Vector) float64 {
  120. if x.Inc < 0 {
  121. panic(negInc)
  122. }
  123. return cblas128.Dzasum(n, x.Data, x.Inc)
  124. }
  125. // Iamax returns the index of an element of x with the largest sum of
  126. // magnitudes of the real and imaginary parts (|Re x[i]|+|Im x[i]|).
  127. // If there are multiple such indices, the earliest is returned.
  128. //
  129. // Iamax returns -1 if n == 0.
  130. //
  131. // Iamax will panic if the vector increment is negative.
  132. func Iamax(n int, x Vector) int {
  133. if x.Inc < 0 {
  134. panic(negInc)
  135. }
  136. return cblas128.Izamax(n, x.Data, x.Inc)
  137. }
  138. // Swap exchanges the elements of two vectors:
  139. // x[i], y[i] = y[i], x[i] for all i.
  140. func Swap(n int, x, y Vector) {
  141. cblas128.Zswap(n, x.Data, x.Inc, y.Data, y.Inc)
  142. }
  143. // Copy copies the elements of x into the elements of y:
  144. // y[i] = x[i] for all i.
  145. func Copy(n int, x, y Vector) {
  146. cblas128.Zcopy(n, x.Data, x.Inc, y.Data, y.Inc)
  147. }
  148. // Axpy computes
  149. // y = alpha * x + y,
  150. // where x and y are vectors, and alpha is a scalar.
  151. func Axpy(n int, alpha complex128, x, y Vector) {
  152. cblas128.Zaxpy(n, alpha, x.Data, x.Inc, y.Data, y.Inc)
  153. }
  154. // Scal computes
  155. // x = alpha * x,
  156. // where x is a vector, and alpha is a scalar.
  157. //
  158. // Scal will panic if the vector increment is negative.
  159. func Scal(n int, alpha complex128, x Vector) {
  160. if x.Inc < 0 {
  161. panic(negInc)
  162. }
  163. cblas128.Zscal(n, alpha, x.Data, x.Inc)
  164. }
  165. // Dscal computes
  166. // x = alpha * x,
  167. // where x is a vector, and alpha is a real scalar.
  168. //
  169. // Dscal will panic if the vector increment is negative.
  170. func Dscal(n int, alpha float64, x Vector) {
  171. if x.Inc < 0 {
  172. panic(negInc)
  173. }
  174. cblas128.Zdscal(n, alpha, x.Data, x.Inc)
  175. }
  176. // Level 2
  177. // Gemv computes
  178. // y = alpha * A * x + beta * y, if t == blas.NoTrans,
  179. // y = alpha * A^T * x + beta * y, if t == blas.Trans,
  180. // y = alpha * A^H * x + beta * y, if t == blas.ConjTrans,
  181. // where A is an m×n dense matrix, x and y are vectors, and alpha and beta are
  182. // scalars.
  183. func Gemv(t blas.Transpose, alpha complex128, a General, x Vector, beta complex128, y Vector) {
  184. cblas128.Zgemv(t, a.Rows, a.Cols, alpha, a.Data, a.Stride, x.Data, x.Inc, beta, y.Data, y.Inc)
  185. }
  186. // Gbmv computes
  187. // y = alpha * A * x + beta * y, if t == blas.NoTrans,
  188. // y = alpha * A^T * x + beta * y, if t == blas.Trans,
  189. // y = alpha * A^H * x + beta * y, if t == blas.ConjTrans,
  190. // where A is an m×n band matrix, x and y are vectors, and alpha and beta are
  191. // scalars.
  192. func Gbmv(t blas.Transpose, alpha complex128, a Band, x Vector, beta complex128, y Vector) {
  193. cblas128.Zgbmv(t, a.Rows, a.Cols, a.KL, a.KU, alpha, a.Data, a.Stride, x.Data, x.Inc, beta, y.Data, y.Inc)
  194. }
  195. // Trmv computes
  196. // x = A * x, if t == blas.NoTrans,
  197. // x = A^T * x, if t == blas.Trans,
  198. // x = A^H * x, if t == blas.ConjTrans,
  199. // where A is an n×n triangular matrix, and x is a vector.
  200. func Trmv(t blas.Transpose, a Triangular, x Vector) {
  201. cblas128.Ztrmv(a.Uplo, t, a.Diag, a.N, a.Data, a.Stride, x.Data, x.Inc)
  202. }
  203. // Tbmv computes
  204. // x = A * x, if t == blas.NoTrans,
  205. // x = A^T * x, if t == blas.Trans,
  206. // x = A^H * x, if t == blas.ConjTrans,
  207. // where A is an n×n triangular band matrix, and x is a vector.
  208. func Tbmv(t blas.Transpose, a TriangularBand, x Vector) {
  209. cblas128.Ztbmv(a.Uplo, t, a.Diag, a.N, a.K, a.Data, a.Stride, x.Data, x.Inc)
  210. }
  211. // Tpmv computes
  212. // x = A * x, if t == blas.NoTrans,
  213. // x = A^T * x, if t == blas.Trans,
  214. // x = A^H * x, if t == blas.ConjTrans,
  215. // where A is an n×n triangular matrix in packed format, and x is a vector.
  216. func Tpmv(t blas.Transpose, a TriangularPacked, x Vector) {
  217. cblas128.Ztpmv(a.Uplo, t, a.Diag, a.N, a.Data, x.Data, x.Inc)
  218. }
  219. // Trsv solves
  220. // A * x = b, if t == blas.NoTrans,
  221. // A^T * x = b, if t == blas.Trans,
  222. // A^H * x = b, if t == blas.ConjTrans,
  223. // where A is an n×n triangular matrix and x is a vector.
  224. //
  225. // At entry to the function, x contains the values of b, and the result is
  226. // stored in-place into x.
  227. //
  228. // No test for singularity or near-singularity is included in this
  229. // routine. Such tests must be performed before calling this routine.
  230. func Trsv(t blas.Transpose, a Triangular, x Vector) {
  231. cblas128.Ztrsv(a.Uplo, t, a.Diag, a.N, a.Data, a.Stride, x.Data, x.Inc)
  232. }
  233. // Tbsv solves
  234. // A * x = b, if t == blas.NoTrans,
  235. // A^T * x = b, if t == blas.Trans,
  236. // A^H * x = b, if t == blas.ConjTrans,
  237. // where A is an n×n triangular band matrix, and x is a vector.
  238. //
  239. // At entry to the function, x contains the values of b, and the result is
  240. // stored in-place into x.
  241. //
  242. // No test for singularity or near-singularity is included in this
  243. // routine. Such tests must be performed before calling this routine.
  244. func Tbsv(t blas.Transpose, a TriangularBand, x Vector) {
  245. cblas128.Ztbsv(a.Uplo, t, a.Diag, a.N, a.K, a.Data, a.Stride, x.Data, x.Inc)
  246. }
  247. // Tpsv solves
  248. // A * x = b, if t == blas.NoTrans,
  249. // A^T * x = b, if t == blas.Trans,
  250. // A^H * x = b, if t == blas.ConjTrans,
  251. // where A is an n×n triangular matrix in packed format and x is a vector.
  252. //
  253. // At entry to the function, x contains the values of b, and the result is
  254. // stored in-place into x.
  255. //
  256. // No test for singularity or near-singularity is included in this
  257. // routine. Such tests must be performed before calling this routine.
  258. func Tpsv(t blas.Transpose, a TriangularPacked, x Vector) {
  259. cblas128.Ztpsv(a.Uplo, t, a.Diag, a.N, a.Data, x.Data, x.Inc)
  260. }
  261. // Hemv computes
  262. // y = alpha * A * x + beta * y,
  263. // where A is an n×n Hermitian matrix, x and y are vectors, and alpha and
  264. // beta are scalars.
  265. func Hemv(alpha complex128, a Hermitian, x Vector, beta complex128, y Vector) {
  266. cblas128.Zhemv(a.Uplo, a.N, alpha, a.Data, a.Stride, x.Data, x.Inc, beta, y.Data, y.Inc)
  267. }
  268. // Hbmv performs
  269. // y = alpha * A * x + beta * y,
  270. // where A is an n×n Hermitian band matrix, x and y are vectors, and alpha
  271. // and beta are scalars.
  272. func Hbmv(alpha complex128, a HermitianBand, x Vector, beta complex128, y Vector) {
  273. cblas128.Zhbmv(a.Uplo, a.N, a.K, alpha, a.Data, a.Stride, x.Data, x.Inc, beta, y.Data, y.Inc)
  274. }
  275. // Hpmv performs
  276. // y = alpha * A * x + beta * y,
  277. // where A is an n×n Hermitian matrix in packed format, x and y are vectors,
  278. // and alpha and beta are scalars.
  279. func Hpmv(alpha complex128, a HermitianPacked, x Vector, beta complex128, y Vector) {
  280. cblas128.Zhpmv(a.Uplo, a.N, alpha, a.Data, x.Data, x.Inc, beta, y.Data, y.Inc)
  281. }
  282. // Geru performs a rank-1 update
  283. // A += alpha * x * y^T,
  284. // where A is an m×n dense matrix, x and y are vectors, and alpha is a scalar.
  285. func Geru(alpha complex128, x, y Vector, a General) {
  286. cblas128.Zgeru(a.Rows, a.Cols, alpha, x.Data, x.Inc, y.Data, y.Inc, a.Data, a.Stride)
  287. }
  288. // Gerc performs a rank-1 update
  289. // A += alpha * x * y^H,
  290. // where A is an m×n dense matrix, x and y are vectors, and alpha is a scalar.
  291. func Gerc(alpha complex128, x, y Vector, a General) {
  292. cblas128.Zgerc(a.Rows, a.Cols, alpha, x.Data, x.Inc, y.Data, y.Inc, a.Data, a.Stride)
  293. }
  294. // Her performs a rank-1 update
  295. // A += alpha * x * y^T,
  296. // where A is an m×n Hermitian matrix, x and y are vectors, and alpha is a scalar.
  297. func Her(alpha float64, x Vector, a Hermitian) {
  298. cblas128.Zher(a.Uplo, a.N, alpha, x.Data, x.Inc, a.Data, a.Stride)
  299. }
  300. // Hpr performs a rank-1 update
  301. // A += alpha * x * x^H,
  302. // where A is an n×n Hermitian matrix in packed format, x is a vector, and
  303. // alpha is a scalar.
  304. func Hpr(alpha float64, x Vector, a HermitianPacked) {
  305. cblas128.Zhpr(a.Uplo, a.N, alpha, x.Data, x.Inc, a.Data)
  306. }
  307. // Her2 performs a rank-2 update
  308. // A += alpha * x * y^H + conj(alpha) * y * x^H,
  309. // where A is an n×n Hermitian matrix, x and y are vectors, and alpha is a scalar.
  310. func Her2(alpha complex128, x, y Vector, a Hermitian) {
  311. cblas128.Zher2(a.Uplo, a.N, alpha, x.Data, x.Inc, y.Data, y.Inc, a.Data, a.Stride)
  312. }
  313. // Hpr2 performs a rank-2 update
  314. // A += alpha * x * y^H + conj(alpha) * y * x^H,
  315. // where A is an n×n Hermitian matrix in packed format, x and y are vectors,
  316. // and alpha is a scalar.
  317. func Hpr2(alpha complex128, x, y Vector, a HermitianPacked) {
  318. cblas128.Zhpr2(a.Uplo, a.N, alpha, x.Data, x.Inc, y.Data, y.Inc, a.Data)
  319. }
  320. // Level 3
  321. // Gemm computes
  322. // C = alpha * A * B + beta * C,
  323. // where A, B, and C are dense matrices, and alpha and beta are scalars.
  324. // tA and tB specify whether A or B are transposed or conjugated.
  325. func Gemm(tA, tB blas.Transpose, alpha complex128, a, b General, beta complex128, c General) {
  326. var m, n, k int
  327. if tA == blas.NoTrans {
  328. m, k = a.Rows, a.Cols
  329. } else {
  330. m, k = a.Cols, a.Rows
  331. }
  332. if tB == blas.NoTrans {
  333. n = b.Cols
  334. } else {
  335. n = b.Rows
  336. }
  337. cblas128.Zgemm(tA, tB, m, n, k, alpha, a.Data, a.Stride, b.Data, b.Stride, beta, c.Data, c.Stride)
  338. }
  339. // Symm performs
  340. // C = alpha * A * B + beta * C, if s == blas.Left,
  341. // C = alpha * B * A + beta * C, if s == blas.Right,
  342. // where A is an n×n or m×m symmetric matrix, B and C are m×n matrices, and
  343. // alpha and beta are scalars.
  344. func Symm(s blas.Side, alpha complex128, a Symmetric, b General, beta complex128, c General) {
  345. var m, n int
  346. if s == blas.Left {
  347. m, n = a.N, b.Cols
  348. } else {
  349. m, n = b.Rows, a.N
  350. }
  351. cblas128.Zsymm(s, a.Uplo, m, n, alpha, a.Data, a.Stride, b.Data, b.Stride, beta, c.Data, c.Stride)
  352. }
  353. // Syrk performs a symmetric rank-k update
  354. // C = alpha * A * A^T + beta * C, if t == blas.NoTrans,
  355. // C = alpha * A^T * A + beta * C, if t == blas.Trans,
  356. // where C is an n×n symmetric matrix, A is an n×k matrix if t == blas.NoTrans
  357. // and a k×n matrix otherwise, and alpha and beta are scalars.
  358. func Syrk(t blas.Transpose, alpha complex128, a General, beta complex128, c Symmetric) {
  359. var n, k int
  360. if t == blas.NoTrans {
  361. n, k = a.Rows, a.Cols
  362. } else {
  363. n, k = a.Cols, a.Rows
  364. }
  365. cblas128.Zsyrk(c.Uplo, t, n, k, alpha, a.Data, a.Stride, beta, c.Data, c.Stride)
  366. }
  367. // Syr2k performs a symmetric rank-2k update
  368. // C = alpha * A * B^T + alpha * B * A^T + beta * C, if t == blas.NoTrans,
  369. // C = alpha * A^T * B + alpha * B^T * A + beta * C, if t == blas.Trans,
  370. // where C is an n×n symmetric matrix, A and B are n×k matrices if
  371. // t == blas.NoTrans and k×n otherwise, and alpha and beta are scalars.
  372. func Syr2k(t blas.Transpose, alpha complex128, a, b General, beta complex128, c Symmetric) {
  373. var n, k int
  374. if t == blas.NoTrans {
  375. n, k = a.Rows, a.Cols
  376. } else {
  377. n, k = a.Cols, a.Rows
  378. }
  379. cblas128.Zsyr2k(c.Uplo, t, n, k, alpha, a.Data, a.Stride, b.Data, b.Stride, beta, c.Data, c.Stride)
  380. }
  381. // Trmm performs
  382. // B = alpha * A * B, if tA == blas.NoTrans and s == blas.Left,
  383. // B = alpha * A^T * B, if tA == blas.Trans and s == blas.Left,
  384. // B = alpha * A^H * B, if tA == blas.ConjTrans and s == blas.Left,
  385. // B = alpha * B * A, if tA == blas.NoTrans and s == blas.Right,
  386. // B = alpha * B * A^T, if tA == blas.Trans and s == blas.Right,
  387. // B = alpha * B * A^H, if tA == blas.ConjTrans and s == blas.Right,
  388. // where A is an n×n or m×m triangular matrix, B is an m×n matrix, and alpha is
  389. // a scalar.
  390. func Trmm(s blas.Side, tA blas.Transpose, alpha complex128, a Triangular, b General) {
  391. cblas128.Ztrmm(s, a.Uplo, tA, a.Diag, b.Rows, b.Cols, alpha, a.Data, a.Stride, b.Data, b.Stride)
  392. }
  393. // Trsm solves
  394. // A * X = alpha * B, if tA == blas.NoTrans and s == blas.Left,
  395. // A^T * X = alpha * B, if tA == blas.Trans and s == blas.Left,
  396. // A^H * X = alpha * B, if tA == blas.ConjTrans and s == blas.Left,
  397. // X * A = alpha * B, if tA == blas.NoTrans and s == blas.Right,
  398. // X * A^T = alpha * B, if tA == blas.Trans and s == blas.Right,
  399. // X * A^H = alpha * B, if tA == blas.ConjTrans and s == blas.Right,
  400. // where A is an n×n or m×m triangular matrix, X and B are m×n matrices, and
  401. // alpha is a scalar.
  402. //
  403. // At entry to the function, b contains the values of B, and the result is
  404. // stored in-place into b.
  405. //
  406. // No check is made that A is invertible.
  407. func Trsm(s blas.Side, tA blas.Transpose, alpha complex128, a Triangular, b General) {
  408. cblas128.Ztrsm(s, a.Uplo, tA, a.Diag, b.Rows, b.Cols, alpha, a.Data, a.Stride, b.Data, b.Stride)
  409. }
  410. // Hemm performs
  411. // C = alpha * A * B + beta * C, if s == blas.Left,
  412. // C = alpha * B * A + beta * C, if s == blas.Right,
  413. // where A is an n×n or m×m Hermitian matrix, B and C are m×n matrices, and
  414. // alpha and beta are scalars.
  415. func Hemm(s blas.Side, alpha complex128, a Hermitian, b General, beta complex128, c General) {
  416. var m, n int
  417. if s == blas.Left {
  418. m, n = a.N, b.Cols
  419. } else {
  420. m, n = b.Rows, a.N
  421. }
  422. cblas128.Zhemm(s, a.Uplo, m, n, alpha, a.Data, a.Stride, b.Data, b.Stride, beta, c.Data, c.Stride)
  423. }
  424. // Herk performs the Hermitian rank-k update
  425. // C = alpha * A * A^H + beta*C, if t == blas.NoTrans,
  426. // C = alpha * A^H * A + beta*C, if t == blas.ConjTrans,
  427. // where C is an n×n Hermitian matrix, A is an n×k matrix if t == blas.NoTrans
  428. // and a k×n matrix otherwise, and alpha and beta are scalars.
  429. func Herk(t blas.Transpose, alpha float64, a General, beta float64, c Hermitian) {
  430. var n, k int
  431. if t == blas.NoTrans {
  432. n, k = a.Rows, a.Cols
  433. } else {
  434. n, k = a.Cols, a.Rows
  435. }
  436. cblas128.Zherk(c.Uplo, t, n, k, alpha, a.Data, a.Stride, beta, c.Data, c.Stride)
  437. }
  438. // Her2k performs the Hermitian rank-2k update
  439. // C = alpha * A * B^H + conj(alpha) * B * A^H + beta * C, if t == blas.NoTrans,
  440. // C = alpha * A^H * B + conj(alpha) * B^H * A + beta * C, if t == blas.ConjTrans,
  441. // where C is an n×n Hermitian matrix, A and B are n×k matrices if t == NoTrans
  442. // and k×n matrices otherwise, and alpha and beta are scalars.
  443. func Her2k(t blas.Transpose, alpha complex128, a, b General, beta float64, c Hermitian) {
  444. var n, k int
  445. if t == blas.NoTrans {
  446. n, k = a.Rows, a.Cols
  447. } else {
  448. n, k = a.Cols, a.Rows
  449. }
  450. cblas128.Zher2k(c.Uplo, t, n, k, alpha, a.Data, a.Stride, b.Data, b.Stride, beta, c.Data, c.Stride)
  451. }