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 cblas64 blas.Complex64 = gonum.Implementation{}
14 // Use sets the BLAS complex64 implementation to be used by subsequent BLAS calls.
15 // The default implementation is cgo.Implementation.
16 func Use(b blas.Complex64) {
20 // Implementation returns the current BLAS complex64 implementation.
22 // Implementation allows direct calls to the current the BLAS complex64 implementation
23 // giving finer control of parameters.
24 func Implementation() blas.Complex64 {
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 {
98 // Hermitian represents an Hermitian matrix using the conventional storage scheme.
99 type Hermitian Symmetric
101 // HermitianBand represents an Hermitian matrix using the band storage scheme.
102 type HermitianBand SymmetricBand
104 // HermitianPacked represents an Hermitian matrix using the packed storage scheme.
105 type HermitianPacked SymmetricPacked
109 const negInc = "cblas64: negative vector increment"
111 // Dotu computes the dot product of the two vectors without
112 // complex conjugation:
114 func Dotu(n int, x, y Vector) complex64 {
115 return cblas64.Cdotu(n, x.Data, x.Inc, y.Data, y.Inc)
118 // Dotc computes the dot product of the two vectors with
119 // complex conjugation:
121 func Dotc(n int, x, y Vector) complex64 {
122 return cblas64.Cdotc(n, x.Data, x.Inc, y.Data, y.Inc)
125 // Nrm2 computes the Euclidean norm of the vector x:
126 // sqrt(\sum_i x[i] * x[i]).
128 // Nrm2 will panic if the vector increment is negative.
129 func Nrm2(n int, x Vector) float32 {
133 return cblas64.Scnrm2(n, x.Data, x.Inc)
136 // Asum computes the sum of magnitudes of the real and imaginary parts of
137 // elements of the vector x:
138 // \sum_i (|Re x[i]| + |Im x[i]|).
140 // Asum will panic if the vector increment is negative.
141 func Asum(n int, x Vector) float32 {
145 return cblas64.Scasum(n, x.Data, x.Inc)
148 // Iamax returns the index of an element of x with the largest sum of
149 // magnitudes of the real and imaginary parts (|Re x[i]|+|Im x[i]|).
150 // If there are multiple such indices, the earliest is returned.
152 // Iamax returns -1 if n == 0.
154 // Iamax will panic if the vector increment is negative.
155 func Iamax(n int, x Vector) int {
159 return cblas64.Icamax(n, x.Data, x.Inc)
162 // Swap exchanges the elements of two vectors:
163 // x[i], y[i] = y[i], x[i] for all i.
164 func Swap(n int, x, y Vector) {
165 cblas64.Cswap(n, x.Data, x.Inc, y.Data, y.Inc)
168 // Copy copies the elements of x into the elements of y:
169 // y[i] = x[i] for all i.
170 func Copy(n int, x, y Vector) {
171 cblas64.Ccopy(n, x.Data, x.Inc, y.Data, y.Inc)
175 // y = alpha * x + y,
176 // where x and y are vectors, and alpha is a scalar.
177 func Axpy(n int, alpha complex64, x, y Vector) {
178 cblas64.Caxpy(n, alpha, x.Data, x.Inc, y.Data, y.Inc)
183 // where x is a vector, and alpha is a scalar.
185 // Scal will panic if the vector increment is negative.
186 func Scal(n int, alpha complex64, x Vector) {
190 cblas64.Cscal(n, alpha, x.Data, x.Inc)
195 // where x is a vector, and alpha is a real scalar.
197 // Dscal will panic if the vector increment is negative.
198 func Dscal(n int, alpha float32, x Vector) {
202 cblas64.Csscal(n, alpha, x.Data, x.Inc)
208 // y = alpha * A * x + beta * y, if t == blas.NoTrans,
209 // y = alpha * A^T * x + beta * y, if t == blas.Trans,
210 // y = alpha * A^H * x + beta * y, if t == blas.ConjTrans,
211 // where A is an m×n dense matrix, x and y are vectors, and alpha and beta are
213 func Gemv(t blas.Transpose, alpha complex64, a General, x Vector, beta complex64, y Vector) {
214 cblas64.Cgemv(t, a.Rows, a.Cols, alpha, a.Data, a.Stride, x.Data, x.Inc, beta, y.Data, y.Inc)
218 // y = alpha * A * x + beta * y, if t == blas.NoTrans,
219 // y = alpha * A^T * x + beta * y, if t == blas.Trans,
220 // y = alpha * A^H * x + beta * y, if t == blas.ConjTrans,
221 // where A is an m×n band matrix, x and y are vectors, and alpha and beta are
223 func Gbmv(t blas.Transpose, alpha complex64, a Band, x Vector, beta complex64, y Vector) {
224 cblas64.Cgbmv(t, a.Rows, a.Cols, a.KL, a.KU, alpha, a.Data, a.Stride, x.Data, x.Inc, beta, y.Data, y.Inc)
228 // x = A * x, if t == blas.NoTrans,
229 // x = A^T * x, if t == blas.Trans,
230 // x = A^H * x, if t == blas.ConjTrans,
231 // where A is an n×n triangular matrix, and x is a vector.
232 func Trmv(t blas.Transpose, a Triangular, x Vector) {
233 cblas64.Ctrmv(a.Uplo, t, a.Diag, a.N, a.Data, a.Stride, x.Data, x.Inc)
237 // x = A * x, if t == blas.NoTrans,
238 // x = A^T * x, if t == blas.Trans,
239 // x = A^H * x, if t == blas.ConjTrans,
240 // where A is an n×n triangular band matrix, and x is a vector.
241 func Tbmv(t blas.Transpose, a TriangularBand, x Vector) {
242 cblas64.Ctbmv(a.Uplo, t, a.Diag, a.N, a.K, a.Data, a.Stride, x.Data, x.Inc)
246 // x = A * x, if t == blas.NoTrans,
247 // x = A^T * x, if t == blas.Trans,
248 // x = A^H * x, if t == blas.ConjTrans,
249 // where A is an n×n triangular matrix in packed format, and x is a vector.
250 func Tpmv(t blas.Transpose, a TriangularPacked, x Vector) {
251 cblas64.Ctpmv(a.Uplo, t, a.Diag, a.N, a.Data, x.Data, x.Inc)
255 // A * x = b, if t == blas.NoTrans,
256 // A^T * x = b, if t == blas.Trans,
257 // A^H * x = b, if t == blas.ConjTrans,
258 // where A is an n×n triangular matrix and x is a vector.
260 // At entry to the function, x contains the values of b, and the result is
261 // stored in-place into x.
263 // No test for singularity or near-singularity is included in this
264 // routine. Such tests must be performed before calling this routine.
265 func Trsv(t blas.Transpose, a Triangular, x Vector) {
266 cblas64.Ctrsv(a.Uplo, t, a.Diag, a.N, a.Data, a.Stride, x.Data, x.Inc)
270 // A * x = b, if t == blas.NoTrans,
271 // A^T * x = b, if t == blas.Trans,
272 // A^H * x = b, if t == blas.ConjTrans,
273 // where A is an n×n triangular band matrix, and x is a vector.
275 // At entry to the function, x contains the values of b, and the result is
276 // stored in-place into x.
278 // No test for singularity or near-singularity is included in this
279 // routine. Such tests must be performed before calling this routine.
280 func Tbsv(t blas.Transpose, a TriangularBand, x Vector) {
281 cblas64.Ctbsv(a.Uplo, t, a.Diag, a.N, a.K, a.Data, a.Stride, x.Data, x.Inc)
285 // A * x = b, if t == blas.NoTrans,
286 // A^T * x = b, if t == blas.Trans,
287 // A^H * x = b, if t == blas.ConjTrans,
288 // where A is an n×n triangular matrix in packed format and x is a vector.
290 // At entry to the function, x contains the values of b, and the result is
291 // stored in-place into x.
293 // No test for singularity or near-singularity is included in this
294 // routine. Such tests must be performed before calling this routine.
295 func Tpsv(t blas.Transpose, a TriangularPacked, x Vector) {
296 cblas64.Ctpsv(a.Uplo, t, a.Diag, a.N, a.Data, x.Data, x.Inc)
300 // y = alpha * A * x + beta * y,
301 // where A is an n×n Hermitian matrix, x and y are vectors, and alpha and
303 func Hemv(alpha complex64, a Hermitian, x Vector, beta complex64, y Vector) {
304 cblas64.Chemv(a.Uplo, a.N, alpha, a.Data, a.Stride, x.Data, x.Inc, beta, y.Data, y.Inc)
308 // y = alpha * A * x + beta * y,
309 // where A is an n×n Hermitian band matrix, x and y are vectors, and alpha
310 // and beta are scalars.
311 func Hbmv(alpha complex64, a HermitianBand, x Vector, beta complex64, y Vector) {
312 cblas64.Chbmv(a.Uplo, a.N, a.K, alpha, a.Data, a.Stride, x.Data, x.Inc, beta, y.Data, y.Inc)
316 // y = alpha * A * x + beta * y,
317 // where A is an n×n Hermitian matrix in packed format, x and y are vectors,
318 // and alpha and beta are scalars.
319 func Hpmv(alpha complex64, a HermitianPacked, x Vector, beta complex64, y Vector) {
320 cblas64.Chpmv(a.Uplo, a.N, alpha, a.Data, x.Data, x.Inc, beta, y.Data, y.Inc)
323 // Geru performs a rank-1 update
324 // A += alpha * x * y^T,
325 // where A is an m×n dense matrix, x and y are vectors, and alpha is a scalar.
326 func Geru(alpha complex64, x, y Vector, a General) {
327 cblas64.Cgeru(a.Rows, a.Cols, alpha, x.Data, x.Inc, y.Data, y.Inc, a.Data, a.Stride)
330 // Gerc performs a rank-1 update
331 // A += alpha * x * y^H,
332 // where A is an m×n dense matrix, x and y are vectors, and alpha is a scalar.
333 func Gerc(alpha complex64, x, y Vector, a General) {
334 cblas64.Cgerc(a.Rows, a.Cols, alpha, x.Data, x.Inc, y.Data, y.Inc, a.Data, a.Stride)
337 // Her performs a rank-1 update
338 // A += alpha * x * y^T,
339 // where A is an m×n Hermitian matrix, x and y are vectors, and alpha is a scalar.
340 func Her(alpha float32, x Vector, a Hermitian) {
341 cblas64.Cher(a.Uplo, a.N, alpha, x.Data, x.Inc, a.Data, a.Stride)
344 // Hpr performs a rank-1 update
345 // A += alpha * x * x^H,
346 // where A is an n×n Hermitian matrix in packed format, x is a vector, and
347 // alpha is a scalar.
348 func Hpr(alpha float32, x Vector, a HermitianPacked) {
349 cblas64.Chpr(a.Uplo, a.N, alpha, x.Data, x.Inc, a.Data)
352 // Her2 performs a rank-2 update
353 // A += alpha * x * y^H + conj(alpha) * y * x^H,
354 // where A is an n×n Hermitian matrix, x and y are vectors, and alpha is a scalar.
355 func Her2(alpha complex64, x, y Vector, a Hermitian) {
356 cblas64.Cher2(a.Uplo, a.N, alpha, x.Data, x.Inc, y.Data, y.Inc, a.Data, a.Stride)
359 // Hpr2 performs a rank-2 update
360 // A += alpha * x * y^H + conj(alpha) * y * x^H,
361 // where A is an n×n Hermitian matrix in packed format, x and y are vectors,
362 // and alpha is a scalar.
363 func Hpr2(alpha complex64, x, y Vector, a HermitianPacked) {
364 cblas64.Chpr2(a.Uplo, a.N, alpha, x.Data, x.Inc, y.Data, y.Inc, a.Data)
370 // C = alpha * A * B + beta * C,
371 // where A, B, and C are dense matrices, and alpha and beta are scalars.
372 // tA and tB specify whether A or B are transposed or conjugated.
373 func Gemm(tA, tB blas.Transpose, alpha complex64, a, b General, beta complex64, c General) {
375 if tA == blas.NoTrans {
376 m, k = a.Rows, a.Cols
378 m, k = a.Cols, a.Rows
380 if tB == blas.NoTrans {
385 cblas64.Cgemm(tA, tB, m, n, k, alpha, a.Data, a.Stride, b.Data, b.Stride, beta, c.Data, c.Stride)
389 // C = alpha * A * B + beta * C, if s == blas.Left,
390 // C = alpha * B * A + beta * C, if s == blas.Right,
391 // where A is an n×n or m×m symmetric matrix, B and C are m×n matrices, and
392 // alpha and beta are scalars.
393 func Symm(s blas.Side, alpha complex64, a Symmetric, b General, beta complex64, c General) {
400 cblas64.Csymm(s, a.Uplo, m, n, alpha, a.Data, a.Stride, b.Data, b.Stride, beta, c.Data, c.Stride)
403 // Syrk performs a symmetric rank-k update
404 // C = alpha * A * A^T + beta * C, if t == blas.NoTrans,
405 // C = alpha * A^T * A + beta * C, if t == blas.Trans,
406 // where C is an n×n symmetric matrix, A is an n×k matrix if t == blas.NoTrans
407 // and a k×n matrix otherwise, and alpha and beta are scalars.
408 func Syrk(t blas.Transpose, alpha complex64, a General, beta complex64, c Symmetric) {
410 if t == blas.NoTrans {
411 n, k = a.Rows, a.Cols
413 n, k = a.Cols, a.Rows
415 cblas64.Csyrk(c.Uplo, t, n, k, alpha, a.Data, a.Stride, beta, c.Data, c.Stride)
418 // Syr2k performs a symmetric rank-2k update
419 // C = alpha * A * B^T + alpha * B * A^T + beta * C, if t == blas.NoTrans,
420 // C = alpha * A^T * B + alpha * B^T * A + beta * C, if t == blas.Trans,
421 // where C is an n×n symmetric matrix, A and B are n×k matrices if
422 // t == blas.NoTrans and k×n otherwise, and alpha and beta are scalars.
423 func Syr2k(t blas.Transpose, alpha complex64, a, b General, beta complex64, c Symmetric) {
425 if t == blas.NoTrans {
426 n, k = a.Rows, a.Cols
428 n, k = a.Cols, a.Rows
430 cblas64.Csyr2k(c.Uplo, t, n, k, alpha, a.Data, a.Stride, b.Data, b.Stride, beta, c.Data, c.Stride)
434 // B = alpha * A * B, if tA == blas.NoTrans and s == blas.Left,
435 // B = alpha * A^T * B, if tA == blas.Trans and s == blas.Left,
436 // B = alpha * A^H * B, if tA == blas.ConjTrans and s == blas.Left,
437 // B = alpha * B * A, if tA == blas.NoTrans and s == blas.Right,
438 // B = alpha * B * A^T, if tA == blas.Trans and s == blas.Right,
439 // B = alpha * B * A^H, if tA == blas.ConjTrans and s == blas.Right,
440 // where A is an n×n or m×m triangular matrix, B is an m×n matrix, and alpha is
442 func Trmm(s blas.Side, tA blas.Transpose, alpha complex64, a Triangular, b General) {
443 cblas64.Ctrmm(s, a.Uplo, tA, a.Diag, b.Rows, b.Cols, alpha, a.Data, a.Stride, b.Data, b.Stride)
447 // A * X = alpha * B, if tA == blas.NoTrans and s == blas.Left,
448 // A^T * X = alpha * B, if tA == blas.Trans and s == blas.Left,
449 // A^H * X = alpha * B, if tA == blas.ConjTrans and s == blas.Left,
450 // X * A = alpha * B, if tA == blas.NoTrans and s == blas.Right,
451 // X * A^T = alpha * B, if tA == blas.Trans and s == blas.Right,
452 // X * A^H = alpha * B, if tA == blas.ConjTrans and s == blas.Right,
453 // where A is an n×n or m×m triangular matrix, X and B are m×n matrices, and
454 // alpha is a scalar.
456 // At entry to the function, b contains the values of B, and the result is
457 // stored in-place into b.
459 // No check is made that A is invertible.
460 func Trsm(s blas.Side, tA blas.Transpose, alpha complex64, a Triangular, b General) {
461 cblas64.Ctrsm(s, a.Uplo, tA, a.Diag, b.Rows, b.Cols, alpha, a.Data, a.Stride, b.Data, b.Stride)
465 // C = alpha * A * B + beta * C, if s == blas.Left,
466 // C = alpha * B * A + beta * C, if s == blas.Right,
467 // where A is an n×n or m×m Hermitian matrix, B and C are m×n matrices, and
468 // alpha and beta are scalars.
469 func Hemm(s blas.Side, alpha complex64, a Hermitian, b General, beta complex64, c General) {
476 cblas64.Chemm(s, a.Uplo, m, n, alpha, a.Data, a.Stride, b.Data, b.Stride, beta, c.Data, c.Stride)
479 // Herk performs the Hermitian rank-k update
480 // C = alpha * A * A^H + beta*C, if t == blas.NoTrans,
481 // C = alpha * A^H * A + beta*C, if t == blas.ConjTrans,
482 // where C is an n×n Hermitian matrix, A is an n×k matrix if t == blas.NoTrans
483 // and a k×n matrix otherwise, and alpha and beta are scalars.
484 func Herk(t blas.Transpose, alpha float32, a General, beta float32, c Hermitian) {
486 if t == blas.NoTrans {
487 n, k = a.Rows, a.Cols
489 n, k = a.Cols, a.Rows
491 cblas64.Cherk(c.Uplo, t, n, k, alpha, a.Data, a.Stride, beta, c.Data, c.Stride)
494 // Her2k performs the Hermitian rank-2k update
495 // C = alpha * A * B^H + conj(alpha) * B * A^H + beta * C, if t == blas.NoTrans,
496 // C = alpha * A^H * B + conj(alpha) * B^H * A + beta * C, if t == blas.ConjTrans,
497 // where C is an n×n Hermitian matrix, A and B are n×k matrices if t == NoTrans
498 // and k×n matrices otherwise, and alpha and beta are scalars.
499 func Her2k(t blas.Transpose, alpha complex64, a, b General, beta float32, c Hermitian) {
501 if t == blas.NoTrans {
502 n, k = a.Rows, a.Cols
504 n, k = a.Cols, a.Rows
506 cblas64.Cher2k(c.Uplo, t, n, k, alpha, a.Data, a.Stride, b.Data, b.Stride, beta, c.Data, c.Stride)