From 11b721c90c9eb185d83e5a4cde94f33e63f3c292 Mon Sep 17 00:00:00 2001 From: Joe Groff Date: Tue, 1 Jul 2008 22:00:22 -0700 Subject: [PATCH] CBLAS library bindings. Factor-ish bindings to better part of level 1 BLAS in math.blas.vectors. Syntax for building BLAS vectors in math.blas.syntax --- extra/math/blas/cblas/authors.txt | 1 + extra/math/blas/cblas/cblas.factor | 557 +++++++++++++++++++ extra/math/blas/cblas/summary.txt | 1 + extra/math/blas/cblas/tags.txt | 2 + extra/math/blas/syntax/authors.txt | 1 + extra/math/blas/syntax/summary.txt | 1 + extra/math/blas/syntax/syntax.factor | 12 + extra/math/blas/syntax/tags.txt | 1 + extra/math/blas/vectors/authors.txt | 1 + extra/math/blas/vectors/summary.txt | 1 + extra/math/blas/vectors/tags.txt | 1 + extra/math/blas/vectors/vectors-tests.factor | 173 ++++++ extra/math/blas/vectors/vectors.factor | 273 +++++++++ 13 files changed, 1025 insertions(+) create mode 100644 extra/math/blas/cblas/authors.txt create mode 100644 extra/math/blas/cblas/cblas.factor create mode 100644 extra/math/blas/cblas/summary.txt create mode 100644 extra/math/blas/cblas/tags.txt create mode 100644 extra/math/blas/syntax/authors.txt create mode 100644 extra/math/blas/syntax/summary.txt create mode 100644 extra/math/blas/syntax/syntax.factor create mode 100644 extra/math/blas/syntax/tags.txt create mode 100644 extra/math/blas/vectors/authors.txt create mode 100644 extra/math/blas/vectors/summary.txt create mode 100644 extra/math/blas/vectors/tags.txt create mode 100644 extra/math/blas/vectors/vectors-tests.factor create mode 100644 extra/math/blas/vectors/vectors.factor diff --git a/extra/math/blas/cblas/authors.txt b/extra/math/blas/cblas/authors.txt new file mode 100644 index 0000000000..f13c9c1e77 --- /dev/null +++ b/extra/math/blas/cblas/authors.txt @@ -0,0 +1 @@ +Joe Groff diff --git a/extra/math/blas/cblas/cblas.factor b/extra/math/blas/cblas/cblas.factor new file mode 100644 index 0000000000..266972fc99 --- /dev/null +++ b/extra/math/blas/cblas/cblas.factor @@ -0,0 +1,557 @@ +USING: alien alien.c-types alien.syntax kernel system combinators ; +IN: math.blas.cblas + +<< "cblas" { + { [ os macosx? ] [ "libcblas.dylib" "cdecl" add-library ] } + { [ os windows? ] [ "cblas.dll" "cdecl" add-library ] } + [ drop "libcblas.so" "cdecl" add-library ] +} cond >> + +LIBRARY: cblas + +TYPEDEF: int CBLAS_ORDER +: CblasRowMajor 101 ; inline +: CblasColMajor 102 ; inline + +TYPEDEF: int CBLAS_TRANSPOSE +: CblasNoTrans 111 ; inline +: CblasTrans 112 ; inline +: CblasConjTrans 113 ; inline + +TYPEDEF: int CBLAS_UPLO +: CblasUpper 121 ; inline +: CblasLower 122 ; inline + +TYPEDEF: int CBLAS_DIAG +: CblasNonUnit 131 ; inline +: CblasUnit 132 ; inline + +TYPEDEF: int CBLAS_SIDE +: CblasLeft 141 ; inline +: CblasRight 142 ; inline + +TYPEDEF: int CBLAS_INDEX + +C-STRUCT: CBLAS_C + { "float" "real" } + { "float" "imag" } ; +C-STRUCT: CBLAS_Z + { "double" "real" } + { "double" "imag" } ; + +! Level 1 BLAS (scalar-vector and vector-vector) + +FUNCTION: float cblas_sdsdot + ( int N, float alpha, float* X, int incX, float* Y, int incY ) ; +FUNCTION: double cblas_dsdot + ( int N, float* X, int incX, float* Y, int incY ) ; +FUNCTION: float cblas_sdot + ( int N, float* X, int incX, float* Y, int incY ) ; +FUNCTION: double cblas_ddot + ( int N, double* X, int incX, double* Y, int incY ) ; + +FUNCTION: void cblas_cdotu_sub + ( int N, CBLAS_C* X, int incX, CBLAS_C* Y, int incY, CBLAS_C* dotu ) ; +FUNCTION: void cblas_cdotc_sub + ( int N, CBLAS_C* X, int incX, CBLAS_C* Y, int incY, CBLAS_C* dotc ) ; + +FUNCTION: void cblas_zdotu_sub + ( int N, CBLAS_Z* X, int incX, CBLAS_Z* Y, int incY, CBLAS_Z* dotu ) ; +FUNCTION: void cblas_zdotc_sub + ( int N, CBLAS_Z* X, int incX, CBLAS_Z* Y, int incY, CBLAS_Z* dotc ) ; + +FUNCTION: float cblas_snrm2 + ( int N, float* X, int incX ) ; +FUNCTION: float cblas_sasum + ( int N, float* X, int incX ) ; + +FUNCTION: double cblas_dnrm2 + ( int N, double* X, int incX ) ; +FUNCTION: double cblas_dasum + ( int N, double* X, int incX ) ; + +FUNCTION: float cblas_scnrm2 + ( int N, CBLAS_C* X, int incX ) ; +FUNCTION: float cblas_scasum + ( int N, CBLAS_C* X, int incX ) ; + +FUNCTION: double cblas_dznrm2 + ( int N, CBLAS_Z* X, int incX ) ; +FUNCTION: double cblas_dzasum + ( int N, CBLAS_Z* X, int incX ) ; + +FUNCTION: CBLAS_INDEX cblas_isamax + ( int N, float* X, int incX ) ; +FUNCTION: CBLAS_INDEX cblas_idamax + ( int N, double* X, int incX ) ; +FUNCTION: CBLAS_INDEX cblas_icamax + ( int N, CBLAS_C* X, int incX ) ; +FUNCTION: CBLAS_INDEX cblas_izamax + ( int N, CBLAS_Z* X, int incX ) ; + +FUNCTION: void cblas_sswap + ( int N, float* X, int incX, float* Y, int incY ) ; +FUNCTION: void cblas_scopy + ( int N, float* X, int incX, float* Y, int incY ) ; +FUNCTION: void cblas_saxpy + ( int N, float alpha, float* X, int incX, float* Y, int incY ) ; + +FUNCTION: void cblas_dswap + ( int N, double* X, int incX, double* Y, int incY ) ; +FUNCTION: void cblas_dcopy + ( int N, double* X, int incX, double* Y, int incY ) ; +FUNCTION: void cblas_daxpy + ( int N, double alpha, double* X, int incX, double* Y, int incY ) ; + +FUNCTION: void cblas_cswap + ( int N, CBLAS_C* X, int incX, CBLAS_C* Y, int incY ) ; +FUNCTION: void cblas_ccopy + ( int N, CBLAS_C* X, int incX, CBLAS_C* Y, int incY ) ; +FUNCTION: void cblas_caxpy + ( int N, CBLAS_C* alpha, CBLAS_C* X, int incX, CBLAS_C* Y, int incY ) ; + +FUNCTION: void cblas_zswap + ( int N, CBLAS_Z* X, int incX, CBLAS_Z* Y, int incY ) ; +FUNCTION: void cblas_zcopy + ( int N, CBLAS_Z* X, int incX, CBLAS_Z* Y, int incY ) ; +FUNCTION: void cblas_zaxpy + ( int N, CBLAS_Z* alpha, CBLAS_Z* X, int incX, CBLAS_Z* Y, int incY ) ; + +FUNCTION: void cblas_sscal + ( int N, float alpha, float* X, int incX ) ; +FUNCTION: void cblas_dscal + ( int N, double alpha, double* X, int incX ) ; +FUNCTION: void cblas_cscal + ( int N, CBLAS_C* alpha, CBLAS_C* X, int incX ) ; +FUNCTION: void cblas_zscal + ( int N, CBLAS_Z* alpha, CBLAS_Z* X, int incX ) ; +FUNCTION: void cblas_csscal + ( int N, float alpha, CBLAS_C* X, int incX ) ; +FUNCTION: void cblas_zdscal + ( int N, double alpha, CBLAS_Z* X, int incX ) ; + +FUNCTION: void cblas_srotg + ( float* a, float* b, float* c, float* s ) ; +FUNCTION: void cblas_srotmg + ( float* d1, float* d2, float* b1, float b2, float* P ) ; +FUNCTION: void cblas_srot + ( int N, float* X, int incX, float* Y, int incY, float c, float s ) ; +FUNCTION: void cblas_srotm + ( int N, float* X, int incX, float* Y, int incY, float* P ) ; + +FUNCTION: void cblas_drotg + ( double* a, double* b, double* c, double* s ) ; +FUNCTION: void cblas_drotmg + ( double* d1, double* d2, double* b1, double b2, double* P ) ; +FUNCTION: void cblas_drot + ( int N, double* X, int incX, double* Y, int incY, double c, double s ) ; +FUNCTION: void cblas_drotm + ( int N, double* X, int incX, double* Y, int incY, double* P ) ; + +! Level 2 BLAS (matrix-vector) + +FUNCTION: void cblas_sgemv ( CBLAS_ORDER Order, + CBLAS_TRANSPOSE TransA, int M, int N, + float alpha, float* A, int lda, + float* X, int incX, float beta, + float* Y, int incY ) ; +FUNCTION: void cblas_sgbmv ( CBLAS_ORDER Order, + CBLAS_TRANSPOSE TransA, int M, int N, + int KL, int KU, float alpha, + float* A, int lda, float* X, + int incX, float beta, float* Y, int incY ) ; +FUNCTION: void cblas_strmv ( CBLAS_ORDER Order, CBLAS_UPLO Uplo, + CBLAS_TRANSPOSE TransA, CBLAS_DIAG Diag, + int N, float* A, int lda, + float* X, int incX ) ; +FUNCTION: void cblas_stbmv ( CBLAS_ORDER Order, CBLAS_UPLO Uplo, + CBLAS_TRANSPOSE TransA, CBLAS_DIAG Diag, + int N, int K, float* A, int lda, + float* X, int incX ) ; +FUNCTION: void cblas_stpmv ( CBLAS_ORDER Order, CBLAS_UPLO Uplo, + CBLAS_TRANSPOSE TransA, CBLAS_DIAG Diag, + int N, float* Ap, float* X, int incX ) ; +FUNCTION: void cblas_strsv ( CBLAS_ORDER Order, CBLAS_UPLO Uplo, + CBLAS_TRANSPOSE TransA, CBLAS_DIAG Diag, + int N, float* A, int lda, float* X, + int incX ) ; +FUNCTION: void cblas_stbsv ( CBLAS_ORDER Order, CBLAS_UPLO Uplo, + CBLAS_TRANSPOSE TransA, CBLAS_DIAG Diag, + int N, int K, float* A, int lda, + float* X, int incX ) ; +FUNCTION: void cblas_stpsv ( CBLAS_ORDER Order, CBLAS_UPLO Uplo, + CBLAS_TRANSPOSE TransA, CBLAS_DIAG Diag, + int N, float* Ap, float* X, int incX ) ; + +FUNCTION: void cblas_dgemv ( CBLAS_ORDER Order, + CBLAS_TRANSPOSE TransA, int M, int N, + double alpha, double* A, int lda, + double* X, int incX, double beta, + double* Y, int incY ) ; +FUNCTION: void cblas_dgbmv ( CBLAS_ORDER Order, + CBLAS_TRANSPOSE TransA, int M, int N, + int KL, int KU, double alpha, + double* A, int lda, double* X, + int incX, double beta, double* Y, int incY ) ; +FUNCTION: void cblas_dtrmv ( CBLAS_ORDER Order, CBLAS_UPLO Uplo, + CBLAS_TRANSPOSE TransA, CBLAS_DIAG Diag, + int N, double* A, int lda, + double* X, int incX ) ; +FUNCTION: void cblas_dtbmv ( CBLAS_ORDER Order, CBLAS_UPLO Uplo, + CBLAS_TRANSPOSE TransA, CBLAS_DIAG Diag, + int N, int K, double* A, int lda, + double* X, int incX ) ; +FUNCTION: void cblas_dtpmv ( CBLAS_ORDER Order, CBLAS_UPLO Uplo, + CBLAS_TRANSPOSE TransA, CBLAS_DIAG Diag, + int N, double* Ap, double* X, int incX ) ; +FUNCTION: void cblas_dtrsv ( CBLAS_ORDER Order, CBLAS_UPLO Uplo, + CBLAS_TRANSPOSE TransA, CBLAS_DIAG Diag, + int N, double* A, int lda, double* X, + int incX ) ; +FUNCTION: void cblas_dtbsv ( CBLAS_ORDER Order, CBLAS_UPLO Uplo, + CBLAS_TRANSPOSE TransA, CBLAS_DIAG Diag, + int N, int K, double* A, int lda, + double* X, int incX ) ; +FUNCTION: void cblas_dtpsv ( CBLAS_ORDER Order, CBLAS_UPLO Uplo, + CBLAS_TRANSPOSE TransA, CBLAS_DIAG Diag, + int N, double* Ap, double* X, int incX ) ; + +FUNCTION: void cblas_cgemv ( CBLAS_ORDER Order, + CBLAS_TRANSPOSE TransA, int M, int N, + void* alpha, void* A, int lda, + void* X, int incX, void* beta, + void* Y, int incY ) ; +FUNCTION: void cblas_cgbmv ( CBLAS_ORDER Order, + CBLAS_TRANSPOSE TransA, int M, int N, + int KL, int KU, void* alpha, + void* A, int lda, void* X, + int incX, void* beta, void* Y, int incY ) ; +FUNCTION: void cblas_ctrmv ( CBLAS_ORDER Order, CBLAS_UPLO Uplo, + CBLAS_TRANSPOSE TransA, CBLAS_DIAG Diag, + int N, void* A, int lda, + void* X, int incX ) ; +FUNCTION: void cblas_ctbmv ( CBLAS_ORDER Order, CBLAS_UPLO Uplo, + CBLAS_TRANSPOSE TransA, CBLAS_DIAG Diag, + int N, int K, void* A, int lda, + void* X, int incX ) ; +FUNCTION: void cblas_ctpmv ( CBLAS_ORDER Order, CBLAS_UPLO Uplo, + CBLAS_TRANSPOSE TransA, CBLAS_DIAG Diag, + int N, void* Ap, void* X, int incX ) ; +FUNCTION: void cblas_ctrsv ( CBLAS_ORDER Order, CBLAS_UPLO Uplo, + CBLAS_TRANSPOSE TransA, CBLAS_DIAG Diag, + int N, void* A, int lda, void* X, + int incX ) ; +FUNCTION: void cblas_ctbsv ( CBLAS_ORDER Order, CBLAS_UPLO Uplo, + CBLAS_TRANSPOSE TransA, CBLAS_DIAG Diag, + int N, int K, void* A, int lda, + void* X, int incX ) ; +FUNCTION: void cblas_ctpsv ( CBLAS_ORDER Order, CBLAS_UPLO Uplo, + CBLAS_TRANSPOSE TransA, CBLAS_DIAG Diag, + int N, void* Ap, void* X, int incX ) ; + +FUNCTION: void cblas_zgemv ( CBLAS_ORDER Order, + CBLAS_TRANSPOSE TransA, int M, int N, + void* alpha, void* A, int lda, + void* X, int incX, void* beta, + void* Y, int incY ) ; +FUNCTION: void cblas_zgbmv ( CBLAS_ORDER Order, + CBLAS_TRANSPOSE TransA, int M, int N, + int KL, int KU, void* alpha, + void* A, int lda, void* X, + int incX, void* beta, void* Y, int incY ) ; +FUNCTION: void cblas_ztrmv ( CBLAS_ORDER Order, CBLAS_UPLO Uplo, + CBLAS_TRANSPOSE TransA, CBLAS_DIAG Diag, + int N, void* A, int lda, + void* X, int incX ) ; +FUNCTION: void cblas_ztbmv ( CBLAS_ORDER Order, CBLAS_UPLO Uplo, + CBLAS_TRANSPOSE TransA, CBLAS_DIAG Diag, + int N, int K, void* A, int lda, + void* X, int incX ) ; +FUNCTION: void cblas_ztpmv ( CBLAS_ORDER Order, CBLAS_UPLO Uplo, + CBLAS_TRANSPOSE TransA, CBLAS_DIAG Diag, + int N, void* Ap, void* X, int incX ) ; +FUNCTION: void cblas_ztrsv ( CBLAS_ORDER Order, CBLAS_UPLO Uplo, + CBLAS_TRANSPOSE TransA, CBLAS_DIAG Diag, + int N, void* A, int lda, void* X, + int incX ) ; +FUNCTION: void cblas_ztbsv ( CBLAS_ORDER Order, CBLAS_UPLO Uplo, + CBLAS_TRANSPOSE TransA, CBLAS_DIAG Diag, + int N, int K, void* A, int lda, + void* X, int incX ) ; +FUNCTION: void cblas_ztpsv ( CBLAS_ORDER Order, CBLAS_UPLO Uplo, + CBLAS_TRANSPOSE TransA, CBLAS_DIAG Diag, + int N, void* Ap, void* X, int incX ) ; + + +FUNCTION: void cblas_ssymv ( CBLAS_ORDER Order, CBLAS_UPLO Uplo, + int N, float alpha, float* A, + int lda, float* X, int incX, + float beta, float* Y, int incY ) ; +FUNCTION: void cblas_ssbmv ( CBLAS_ORDER Order, CBLAS_UPLO Uplo, + int N, int K, float alpha, float* A, + int lda, float* X, int incX, + float beta, float* Y, int incY ) ; +FUNCTION: void cblas_sspmv ( CBLAS_ORDER Order, CBLAS_UPLO Uplo, + int N, float alpha, float* Ap, + float* X, int incX, + float beta, float* Y, int incY ) ; +FUNCTION: void cblas_sger ( CBLAS_ORDER Order, int M, int N, + float alpha, float* X, int incX, + float* Y, int incY, float* A, int lda ) ; +FUNCTION: void cblas_ssyr ( CBLAS_ORDER Order, CBLAS_UPLO Uplo, + int N, float alpha, float* X, + int incX, float* A, int lda ) ; +FUNCTION: void cblas_sspr ( CBLAS_ORDER Order, CBLAS_UPLO Uplo, + int N, float alpha, float* X, + int incX, float* Ap ) ; +FUNCTION: void cblas_ssyr2 ( CBLAS_ORDER Order, CBLAS_UPLO Uplo, + int N, float alpha, float* X, + int incX, float* Y, int incY, float* A, + int lda ) ; +FUNCTION: void cblas_sspr2 ( CBLAS_ORDER Order, CBLAS_UPLO Uplo, + int N, float alpha, float* X, + int incX, float* Y, int incY, float* A ) ; + +FUNCTION: void cblas_dsymv ( CBLAS_ORDER Order, CBLAS_UPLO Uplo, + int N, double alpha, double* A, + int lda, double* X, int incX, + double beta, double* Y, int incY ) ; +FUNCTION: void cblas_dsbmv ( CBLAS_ORDER Order, CBLAS_UPLO Uplo, + int N, int K, double alpha, double* A, + int lda, double* X, int incX, + double beta, double* Y, int incY ) ; +FUNCTION: void cblas_dspmv ( CBLAS_ORDER Order, CBLAS_UPLO Uplo, + int N, double alpha, double* Ap, + double* X, int incX, + double beta, double* Y, int incY ) ; +FUNCTION: void cblas_dger ( CBLAS_ORDER Order, int M, int N, + double alpha, double* X, int incX, + double* Y, int incY, double* A, int lda ) ; +FUNCTION: void cblas_dsyr ( CBLAS_ORDER Order, CBLAS_UPLO Uplo, + int N, double alpha, double* X, + int incX, double* A, int lda ) ; +FUNCTION: void cblas_dspr ( CBLAS_ORDER Order, CBLAS_UPLO Uplo, + int N, double alpha, double* X, + int incX, double* Ap ) ; +FUNCTION: void cblas_dsyr2 ( CBLAS_ORDER Order, CBLAS_UPLO Uplo, + int N, double alpha, double* X, + int incX, double* Y, int incY, double* A, + int lda ) ; +FUNCTION: void cblas_dspr2 ( CBLAS_ORDER Order, CBLAS_UPLO Uplo, + int N, double alpha, double* X, + int incX, double* Y, int incY, double* A ) ; + + +FUNCTION: void cblas_chemv ( CBLAS_ORDER Order, CBLAS_UPLO Uplo, + int N, void* alpha, void* A, + int lda, void* X, int incX, + void* beta, void* Y, int incY ) ; +FUNCTION: void cblas_chbmv ( CBLAS_ORDER Order, CBLAS_UPLO Uplo, + int N, int K, void* alpha, void* A, + int lda, void* X, int incX, + void* beta, void* Y, int incY ) ; +FUNCTION: void cblas_chpmv ( CBLAS_ORDER Order, CBLAS_UPLO Uplo, + int N, void* alpha, void* Ap, + void* X, int incX, + void* beta, void* Y, int incY ) ; +FUNCTION: void cblas_cgeru ( CBLAS_ORDER Order, int M, int N, + void* alpha, void* X, int incX, + void* Y, int incY, void* A, int lda ) ; +FUNCTION: void cblas_cgerc ( CBLAS_ORDER Order, int M, int N, + void* alpha, void* X, int incX, + void* Y, int incY, void* A, int lda ) ; +FUNCTION: void cblas_cher ( CBLAS_ORDER Order, CBLAS_UPLO Uplo, + int N, float alpha, void* X, int incX, + void* A, int lda ) ; +FUNCTION: void cblas_chpr ( CBLAS_ORDER Order, CBLAS_UPLO Uplo, + int N, float alpha, void* X, + int incX, void* A ) ; +FUNCTION: void cblas_cher2 ( CBLAS_ORDER Order, CBLAS_UPLO Uplo, int N, + void* alpha, void* X, int incX, + void* Y, int incY, void* A, int lda ) ; +FUNCTION: void cblas_chpr2 ( CBLAS_ORDER Order, CBLAS_UPLO Uplo, int N, + void* alpha, void* X, int incX, + void* Y, int incY, void* Ap ) ; + +FUNCTION: void cblas_zhemv ( CBLAS_ORDER Order, CBLAS_UPLO Uplo, + int N, void* alpha, void* A, + int lda, void* X, int incX, + void* beta, void* Y, int incY ) ; +FUNCTION: void cblas_zhbmv ( CBLAS_ORDER Order, CBLAS_UPLO Uplo, + int N, int K, void* alpha, void* A, + int lda, void* X, int incX, + void* beta, void* Y, int incY ) ; +FUNCTION: void cblas_zhpmv ( CBLAS_ORDER Order, CBLAS_UPLO Uplo, + int N, void* alpha, void* Ap, + void* X, int incX, + void* beta, void* Y, int incY ) ; +FUNCTION: void cblas_zgeru ( CBLAS_ORDER Order, int M, int N, + void* alpha, void* X, int incX, + void* Y, int incY, void* A, int lda ) ; +FUNCTION: void cblas_zgerc ( CBLAS_ORDER Order, int M, int N, + void* alpha, void* X, int incX, + void* Y, int incY, void* A, int lda ) ; +FUNCTION: void cblas_zher ( CBLAS_ORDER Order, CBLAS_UPLO Uplo, + int N, double alpha, void* X, int incX, + void* A, int lda ) ; +FUNCTION: void cblas_zhpr ( CBLAS_ORDER Order, CBLAS_UPLO Uplo, + int N, double alpha, void* X, + int incX, void* A ) ; +FUNCTION: void cblas_zher2 ( CBLAS_ORDER Order, CBLAS_UPLO Uplo, int N, + void* alpha, void* X, int incX, + void* Y, int incY, void* A, int lda ) ; +FUNCTION: void cblas_zhpr2 ( CBLAS_ORDER Order, CBLAS_UPLO Uplo, int N, + void* alpha, void* X, int incX, + void* Y, int incY, void* Ap ) ; + +! Level 3 BLAS (matrix-matrix) + +FUNCTION: void cblas_sgemm ( CBLAS_ORDER Order, CBLAS_TRANSPOSE TransA, + CBLAS_TRANSPOSE TransB, int M, int N, + int K, float alpha, float* A, + int lda, float* B, int ldb, + float beta, float* C, int ldc ) ; +FUNCTION: void cblas_ssymm ( CBLAS_ORDER Order, CBLAS_SIDE Side, + CBLAS_UPLO Uplo, int M, int N, + float alpha, float* A, int lda, + float* B, int ldb, float beta, + float* C, int ldc ) ; +FUNCTION: void cblas_ssyrk ( CBLAS_ORDER Order, CBLAS_UPLO Uplo, + CBLAS_TRANSPOSE Trans, int N, int K, + float alpha, float* A, int lda, + float beta, float* C, int ldc ) ; +FUNCTION: void cblas_ssyr2k ( CBLAS_ORDER Order, CBLAS_UPLO Uplo, + CBLAS_TRANSPOSE Trans, int N, int K, + float alpha, float* A, int lda, + float* B, int ldb, float beta, + float* C, int ldc ) ; +FUNCTION: void cblas_strmm ( CBLAS_ORDER Order, CBLAS_SIDE Side, + CBLAS_UPLO Uplo, CBLAS_TRANSPOSE TransA, + CBLAS_DIAG Diag, int M, int N, + float alpha, float* A, int lda, + float* B, int ldb ) ; +FUNCTION: void cblas_strsm ( CBLAS_ORDER Order, CBLAS_SIDE Side, + CBLAS_UPLO Uplo, CBLAS_TRANSPOSE TransA, + CBLAS_DIAG Diag, int M, int N, + float alpha, float* A, int lda, + float* B, int ldb ) ; + +FUNCTION: void cblas_dgemm ( CBLAS_ORDER Order, CBLAS_TRANSPOSE TransA, + CBLAS_TRANSPOSE TransB, int M, int N, + int K, double alpha, double* A, + int lda, double* B, int ldb, + double beta, double* C, int ldc ) ; +FUNCTION: void cblas_dsymm ( CBLAS_ORDER Order, CBLAS_SIDE Side, + CBLAS_UPLO Uplo, int M, int N, + double alpha, double* A, int lda, + double* B, int ldb, double beta, + double* C, int ldc ) ; +FUNCTION: void cblas_dsyrk ( CBLAS_ORDER Order, CBLAS_UPLO Uplo, + CBLAS_TRANSPOSE Trans, int N, int K, + double alpha, double* A, int lda, + double beta, double* C, int ldc ) ; +FUNCTION: void cblas_dsyr2k ( CBLAS_ORDER Order, CBLAS_UPLO Uplo, + CBLAS_TRANSPOSE Trans, int N, int K, + double alpha, double* A, int lda, + double* B, int ldb, double beta, + double* C, int ldc ) ; +FUNCTION: void cblas_dtrmm ( CBLAS_ORDER Order, CBLAS_SIDE Side, + CBLAS_UPLO Uplo, CBLAS_TRANSPOSE TransA, + CBLAS_DIAG Diag, int M, int N, + double alpha, double* A, int lda, + double* B, int ldb ) ; +FUNCTION: void cblas_dtrsm ( CBLAS_ORDER Order, CBLAS_SIDE Side, + CBLAS_UPLO Uplo, CBLAS_TRANSPOSE TransA, + CBLAS_DIAG Diag, int M, int N, + double alpha, double* A, int lda, + double* B, int ldb ) ; + +FUNCTION: void cblas_cgemm ( CBLAS_ORDER Order, CBLAS_TRANSPOSE TransA, + CBLAS_TRANSPOSE TransB, int M, int N, + int K, void* alpha, void* A, + int lda, void* B, int ldb, + void* beta, void* C, int ldc ) ; +FUNCTION: void cblas_csymm ( CBLAS_ORDER Order, CBLAS_SIDE Side, + CBLAS_UPLO Uplo, int M, int N, + void* alpha, void* A, int lda, + void* B, int ldb, void* beta, + void* C, int ldc ) ; +FUNCTION: void cblas_csyrk ( CBLAS_ORDER Order, CBLAS_UPLO Uplo, + CBLAS_TRANSPOSE Trans, int N, int K, + void* alpha, void* A, int lda, + void* beta, void* C, int ldc ) ; +FUNCTION: void cblas_csyr2k ( CBLAS_ORDER Order, CBLAS_UPLO Uplo, + CBLAS_TRANSPOSE Trans, int N, int K, + void* alpha, void* A, int lda, + void* B, int ldb, void* beta, + void* C, int ldc ) ; +FUNCTION: void cblas_ctrmm ( CBLAS_ORDER Order, CBLAS_SIDE Side, + CBLAS_UPLO Uplo, CBLAS_TRANSPOSE TransA, + CBLAS_DIAG Diag, int M, int N, + void* alpha, void* A, int lda, + void* B, int ldb ) ; +FUNCTION: void cblas_ctrsm ( CBLAS_ORDER Order, CBLAS_SIDE Side, + CBLAS_UPLO Uplo, CBLAS_TRANSPOSE TransA, + CBLAS_DIAG Diag, int M, int N, + void* alpha, void* A, int lda, + void* B, int ldb ) ; + +FUNCTION: void cblas_zgemm ( CBLAS_ORDER Order, CBLAS_TRANSPOSE TransA, + CBLAS_TRANSPOSE TransB, int M, int N, + int K, void* alpha, void* A, + int lda, void* B, int ldb, + void* beta, void* C, int ldc ) ; +FUNCTION: void cblas_zsymm ( CBLAS_ORDER Order, CBLAS_SIDE Side, + CBLAS_UPLO Uplo, int M, int N, + void* alpha, void* A, int lda, + void* B, int ldb, void* beta, + void* C, int ldc ) ; +FUNCTION: void cblas_zsyrk ( CBLAS_ORDER Order, CBLAS_UPLO Uplo, + CBLAS_TRANSPOSE Trans, int N, int K, + void* alpha, void* A, int lda, + void* beta, void* C, int ldc ) ; +FUNCTION: void cblas_zsyr2k ( CBLAS_ORDER Order, CBLAS_UPLO Uplo, + CBLAS_TRANSPOSE Trans, int N, int K, + void* alpha, void* A, int lda, + void* B, int ldb, void* beta, + void* C, int ldc ) ; +FUNCTION: void cblas_ztrmm ( CBLAS_ORDER Order, CBLAS_SIDE Side, + CBLAS_UPLO Uplo, CBLAS_TRANSPOSE TransA, + CBLAS_DIAG Diag, int M, int N, + void* alpha, void* A, int lda, + void* B, int ldb ) ; +FUNCTION: void cblas_ztrsm ( CBLAS_ORDER Order, CBLAS_SIDE Side, + CBLAS_UPLO Uplo, CBLAS_TRANSPOSE TransA, + CBLAS_DIAG Diag, int M, int N, + void* alpha, void* A, int lda, + void* B, int ldb ) ; + +FUNCTION: void cblas_chemm ( CBLAS_ORDER Order, CBLAS_SIDE Side, + CBLAS_UPLO Uplo, int M, int N, + void* alpha, void* A, int lda, + void* B, int ldb, void* beta, + void* C, int ldc ) ; +FUNCTION: void cblas_cherk ( CBLAS_ORDER Order, CBLAS_UPLO Uplo, + CBLAS_TRANSPOSE Trans, int N, int K, + float alpha, void* A, int lda, + float beta, void* C, int ldc ) ; +FUNCTION: void cblas_cher2k ( CBLAS_ORDER Order, CBLAS_UPLO Uplo, + CBLAS_TRANSPOSE Trans, int N, int K, + void* alpha, void* A, int lda, + void* B, int ldb, float beta, + void* C, int ldc ) ; +FUNCTION: void cblas_zhemm ( CBLAS_ORDER Order, CBLAS_SIDE Side, + CBLAS_UPLO Uplo, int M, int N, + void* alpha, void* A, int lda, + void* B, int ldb, void* beta, + void* C, int ldc ) ; +FUNCTION: void cblas_zherk ( CBLAS_ORDER Order, CBLAS_UPLO Uplo, + CBLAS_TRANSPOSE Trans, int N, int K, + double alpha, void* A, int lda, + double beta, void* C, int ldc ) ; +FUNCTION: void cblas_zher2k ( CBLAS_ORDER Order, CBLAS_UPLO Uplo, + CBLAS_TRANSPOSE Trans, int N, int K, + void* alpha, void* A, int lda, + void* B, int ldb, double beta, + void* C, int ldc ) ; + diff --git a/extra/math/blas/cblas/summary.txt b/extra/math/blas/cblas/summary.txt new file mode 100644 index 0000000000..c72e78eb0d --- /dev/null +++ b/extra/math/blas/cblas/summary.txt @@ -0,0 +1 @@ +Low-level bindings to the C Basic Linear Algebra Subroutines (BLAS) library diff --git a/extra/math/blas/cblas/tags.txt b/extra/math/blas/cblas/tags.txt new file mode 100644 index 0000000000..241ec1ecda --- /dev/null +++ b/extra/math/blas/cblas/tags.txt @@ -0,0 +1,2 @@ +math +bindings diff --git a/extra/math/blas/syntax/authors.txt b/extra/math/blas/syntax/authors.txt new file mode 100644 index 0000000000..f13c9c1e77 --- /dev/null +++ b/extra/math/blas/syntax/authors.txt @@ -0,0 +1 @@ +Joe Groff diff --git a/extra/math/blas/syntax/summary.txt b/extra/math/blas/syntax/summary.txt new file mode 100644 index 0000000000..a71bebb50f --- /dev/null +++ b/extra/math/blas/syntax/summary.txt @@ -0,0 +1 @@ +Literal syntax for BLAS vectors and matrices diff --git a/extra/math/blas/syntax/syntax.factor b/extra/math/blas/syntax/syntax.factor new file mode 100644 index 0000000000..d161739d80 --- /dev/null +++ b/extra/math/blas/syntax/syntax.factor @@ -0,0 +1,12 @@ +USING: kernel math.blas.vectors parser prettyprint.backend ; +IN: math.blas.syntax + +: svector{ ( accum -- accum ) + \ } [ >float-blas-vector ] parse-literal ; parsing +: dvector{ ( accum -- accum ) + \ } [ >double-blas-vector ] parse-literal ; parsing +: cvector{ ( accum -- accum ) + \ } [ >float-complex-blas-vector ] parse-literal ; parsing +: zvector{ ( accum -- accum ) + \ } [ >double-complex-blas-vector ] parse-literal ; parsing + diff --git a/extra/math/blas/syntax/tags.txt b/extra/math/blas/syntax/tags.txt new file mode 100644 index 0000000000..ede10ab61b --- /dev/null +++ b/extra/math/blas/syntax/tags.txt @@ -0,0 +1 @@ +math diff --git a/extra/math/blas/vectors/authors.txt b/extra/math/blas/vectors/authors.txt new file mode 100644 index 0000000000..f13c9c1e77 --- /dev/null +++ b/extra/math/blas/vectors/authors.txt @@ -0,0 +1 @@ +Joe Groff diff --git a/extra/math/blas/vectors/summary.txt b/extra/math/blas/vectors/summary.txt new file mode 100644 index 0000000000..91653e0938 --- /dev/null +++ b/extra/math/blas/vectors/summary.txt @@ -0,0 +1 @@ +Basic Linear Algebra words for accelerated vector and matrix math diff --git a/extra/math/blas/vectors/tags.txt b/extra/math/blas/vectors/tags.txt new file mode 100644 index 0000000000..ede10ab61b --- /dev/null +++ b/extra/math/blas/vectors/tags.txt @@ -0,0 +1 @@ +math diff --git a/extra/math/blas/vectors/vectors-tests.factor b/extra/math/blas/vectors/vectors-tests.factor new file mode 100644 index 0000000000..e059d2943d --- /dev/null +++ b/extra/math/blas/vectors/vectors-tests.factor @@ -0,0 +1,173 @@ +USING: kernel math.blas.vectors math.blas.syntax sequences tools.test ; +IN: math.blas.vectors.tests + +! clone + +[ svector{ 1.0 2.0 3.0 } ] [ svector{ 1.0 2.0 3.0 } clone ] unit-test +[ f ] [ svector{ 1.0 2.0 3.0 } dup clone eq? ] unit-test +[ dvector{ 1.0 2.0 3.0 } ] [ dvector{ 1.0 2.0 3.0 } clone ] unit-test +[ f ] [ dvector{ 1.0 2.0 3.0 } dup clone eq? ] unit-test +[ cvector{ 1.0 C{ 2.0 3.0 } 4.0 } ] [ cvector{ 1.0 C{ 2.0 3.0 } 4.0 } clone ] unit-test +[ f ] [ cvector{ 1.0 C{ 2.0 3.0 } 4.0 } dup clone eq? ] unit-test +[ zvector{ 1.0 C{ 2.0 3.0 } 4.0 } ] [ zvector{ 1.0 C{ 2.0 3.0 } 4.0 } clone ] unit-test +[ f ] [ zvector{ 1.0 C{ 2.0 3.0 } 4.0 } dup clone eq? ] unit-test + +! nth + +[ 1.0 ] [ 2 svector{ 3.0 2.0 1.0 } nth ] unit-test +[ 1.0 ] [ 2 dvector{ 3.0 2.0 1.0 } nth ] unit-test + +[ C{ 1.0 2.0 } ] +[ 2 cvector{ C{ -3.0 -2.0 } C{ -1.0 0.0 } C{ 1.0 2.0 } } nth ] unit-test + +[ C{ 1.0 2.0 } ] +[ 2 zvector{ C{ -3.0 -2.0 } C{ -1.0 0.0 } C{ 1.0 2.0 } } nth ] unit-test + +! set-nth + +[ svector{ 3.0 2.0 0.0 } ] [ 0.0 2 svector{ 3.0 2.0 1.0 } [ set-nth ] keep ] unit-test +[ dvector{ 3.0 2.0 0.0 } ] [ 0.0 2 dvector{ 3.0 2.0 1.0 } [ set-nth ] keep ] unit-test + +[ cvector{ C{ -3.0 -2.0 } C{ -1.0 0.0 } C{ 3.0 4.0 } } ] [ + C{ 3.0 4.0 } 2 + cvector{ C{ -3.0 -2.0 } C{ -1.0 0.0 } C{ 1.0 2.0 } } + [ set-nth ] keep +] unit-test +[ zvector{ C{ -3.0 -2.0 } C{ -1.0 0.0 } C{ 3.0 4.0 } } ] [ + C{ 3.0 4.0 } 2 + zvector{ C{ -3.0 -2.0 } C{ -1.0 0.0 } C{ 1.0 2.0 } } + [ set-nth ] keep +] unit-test + +! V+ + +[ svector{ 11.0 22.0 } ] [ svector{ 1.0 2.0 } svector{ 10.0 20.0 } V+ ] unit-test +[ dvector{ 11.0 22.0 } ] [ dvector{ 1.0 2.0 } dvector{ 10.0 20.0 } V+ ] unit-test + +[ cvector{ 11.0 C{ 22.0 33.0 } } ] +[ cvector{ 1.0 C{ 2.0 3.0 } } cvector{ 10.0 C{ 20.0 30.0 } } V+ ] +unit-test + +[ zvector{ 11.0 C{ 22.0 33.0 } } ] +[ zvector{ 1.0 C{ 2.0 3.0 } } zvector{ 10.0 C{ 20.0 30.0 } } V+ ] +unit-test + +! V- + +[ svector{ 9.0 18.0 } ] [ svector{ 10.0 20.0 } svector{ 1.0 2.0 } V- ] unit-test +[ dvector{ 9.0 18.0 } ] [ dvector{ 10.0 20.0 } dvector{ 1.0 2.0 } V- ] unit-test + +[ cvector{ 9.0 C{ 18.0 27.0 } } ] +[ cvector{ 10.0 C{ 20.0 30.0 } } cvector{ 1.0 C{ 2.0 3.0 } } V- ] +unit-test + +[ zvector{ 9.0 C{ 18.0 27.0 } } ] +[ zvector{ 10.0 C{ 20.0 30.0 } } zvector{ 1.0 C{ 2.0 3.0 } } V- ] +unit-test + +! Vneg + +[ svector{ 1.0 -2.0 } ] [ svector{ -1.0 2.0 } Vneg ] unit-test +[ dvector{ 1.0 -2.0 } ] [ dvector{ -1.0 2.0 } Vneg ] unit-test + +[ cvector{ 1.0 C{ -2.0 3.0 } } ] [ cvector{ -1.0 C{ 2.0 -3.0 } } Vneg ] unit-test +[ zvector{ 1.0 C{ -2.0 3.0 } } ] [ zvector{ -1.0 C{ 2.0 -3.0 } } Vneg ] unit-test + +! n*V + +[ svector{ 100.0 200.0 } ] [ 10.0 svector{ 10.0 20.0 } n*V ] unit-test +[ dvector{ 100.0 200.0 } ] [ 10.0 dvector{ 10.0 20.0 } n*V ] unit-test + +[ cvector{ C{ 20.0 4.0 } C{ 8.0 12.0 } } ] +[ C{ 10.0 2.0 } cvector{ 2.0 C{ 1.0 1.0 } } n*V ] +unit-test + +[ zvector{ C{ 20.0 4.0 } C{ 8.0 12.0 } } ] +[ C{ 10.0 2.0 } zvector{ 2.0 C{ 1.0 1.0 } } n*V ] +unit-test + +! V*n + +[ svector{ 100.0 200.0 } ] [ svector{ 10.0 20.0 } 10.0 V*n ] unit-test +[ dvector{ 100.0 200.0 } ] [ dvector{ 10.0 20.0 } 10.0 V*n ] unit-test + +[ cvector{ C{ 20.0 4.0 } C{ 8.0 12.0 } } ] +[ cvector{ 2.0 C{ 1.0 1.0 } } C{ 10.0 2.0 } V*n ] +unit-test + +[ zvector{ C{ 20.0 4.0 } C{ 8.0 12.0 } } ] +[ zvector{ 2.0 C{ 1.0 1.0 } } C{ 10.0 2.0 } V*n ] +unit-test + +! V/n + +[ svector{ 1.0 2.0 } ] [ svector{ 4.0 8.0 } 4.0 V/n ] unit-test +[ dvector{ 1.0 2.0 } ] [ dvector{ 4.0 8.0 } 4.0 V/n ] unit-test + +[ cvector{ 2.0 1.0 } ] +[ cvector{ C{ 16.0 4.0 } C{ 8.0 2.0 } } C{ 8.0 2.0 } V/n ] +unit-test + +[ cvector{ 2.0 1.0 } ] +[ cvector{ C{ 16.0 4.0 } C{ 8.0 2.0 } } C{ 8.0 2.0 } V/n ] +unit-test + +! V. + +[ 7.0 ] [ svector{ 1.0 2.5 } svector{ 2.0 2.0 } V. ] unit-test +[ 7.0 ] [ dvector{ 1.0 2.5 } dvector{ 2.0 2.0 } V. ] unit-test +[ C{ 7.0 7.0 } ] [ cvector{ C{ 1.0 1.0 } 2.5 } cvector{ 2.0 C{ 2.0 2.0 } } V. ] unit-test +[ C{ 7.0 7.0 } ] [ zvector{ C{ 1.0 1.0 } 2.5 } zvector{ 2.0 C{ 2.0 2.0 } } V. ] unit-test + +! V.conj + +[ C{ 7.0 3.0 } ] [ cvector{ C{ 1.0 1.0 } 2.5 } cvector{ 2.0 C{ 2.0 2.0 } } V.conj ] unit-test +[ C{ 7.0 3.0 } ] [ zvector{ C{ 1.0 1.0 } 2.5 } zvector{ 2.0 C{ 2.0 2.0 } } V.conj ] unit-test + +! Vnorm + +[ 5.0 ] [ svector{ 3.0 4.0 } Vnorm ] unit-test +[ 5.0 ] [ dvector{ 3.0 4.0 } Vnorm ] unit-test + +[ 13.0 ] [ cvector{ C{ 3.0 4.0 } 12.0 } Vnorm ] unit-test +[ 13.0 ] [ zvector{ C{ 3.0 4.0 } 12.0 } Vnorm ] unit-test + +! Vasum + +[ 6.0 ] [ svector{ 1.0 2.0 -3.0 } Vasum ] unit-test +[ 6.0 ] [ dvector{ 1.0 2.0 -3.0 } Vasum ] unit-test + +[ 15.0 ] [ cvector{ 1.0 C{ -2.0 3.0 } C{ 4.0 -5.0 } } Vasum ] unit-test +[ 15.0 ] [ zvector{ 1.0 C{ -2.0 3.0 } C{ 4.0 -5.0 } } Vasum ] unit-test + +! Vswap + +[ svector{ 2.0 2.0 } svector{ 1.0 1.0 } ] +[ svector{ 1.0 1.0 } svector{ 2.0 2.0 } Vswap ] +unit-test + +[ dvector{ 2.0 2.0 } dvector{ 1.0 1.0 } ] +[ dvector{ 1.0 1.0 } dvector{ 2.0 2.0 } Vswap ] +unit-test + +[ cvector{ 2.0 C{ 2.0 2.0 } } cvector{ C{ 1.0 1.0 } 1.0 } ] +[ cvector{ C{ 1.0 1.0 } 1.0 } cvector{ 2.0 C{ 2.0 2.0 } } Vswap ] +unit-test + +[ zvector{ 2.0 C{ 2.0 2.0 } } zvector{ C{ 1.0 1.0 } 1.0 } ] +[ zvector{ C{ 1.0 1.0 } 1.0 } zvector{ 2.0 C{ 2.0 2.0 } } Vswap ] +unit-test + +! Viamax + +[ 3 ] [ svector{ 1.0 -5.0 4.0 -6.0 -1.0 } Viamax ] unit-test +[ 3 ] [ dvector{ 1.0 -5.0 4.0 -6.0 -1.0 } Viamax ] unit-test +[ 0 ] [ cvector{ C{ 2.0 -5.0 } 4.0 -6.0 -1.0 } Viamax ] unit-test +[ 0 ] [ zvector{ C{ 2.0 -5.0 } 4.0 -6.0 -1.0 } Viamax ] unit-test + +! Vamax + +[ -6.0 ] [ svector{ 1.0 -5.0 4.0 -6.0 -1.0 } Vamax ] unit-test +[ -6.0 ] [ dvector{ 1.0 -5.0 4.0 -6.0 -1.0 } Vamax ] unit-test +[ C{ 2.0 -5.0 } ] [ cvector{ C{ 2.0 -5.0 } 4.0 -6.0 -1.0 } Vamax ] unit-test +[ C{ 2.0 -5.0 } ] [ zvector{ C{ 2.0 -5.0 } 4.0 -6.0 -1.0 } Vamax ] unit-test diff --git a/extra/math/blas/vectors/vectors.factor b/extra/math/blas/vectors/vectors.factor new file mode 100644 index 0000000000..bec1daa855 --- /dev/null +++ b/extra/math/blas/vectors/vectors.factor @@ -0,0 +1,273 @@ +USING: accessors alien alien.c-types arrays byte-arrays fry +kernel macros math math.blas.cblas math.complex math.functions +math.order multi-methods qualified sequences sequences.private +shuffle ; +QUALIFIED: syntax +IN: math.blas.vectors + +TUPLE: blas-vector-base data length inc ; +TUPLE: float-blas-vector < blas-vector-base ; +TUPLE: double-blas-vector < blas-vector-base ; +TUPLE: float-complex-blas-vector < blas-vector-base ; +TUPLE: double-complex-blas-vector < blas-vector-base ; + +INSTANCE: float-blas-vector sequence +INSTANCE: double-blas-vector sequence +INSTANCE: float-complex-blas-vector sequence +INSTANCE: double-complex-blas-vector sequence + +C: float-blas-vector +C: double-blas-vector +C: float-complex-blas-vector +C: double-complex-blas-vector + +GENERIC: zero-vector ( v -- zero ) + +GENERIC: n*V+V-in-place ( n v1 v2 -- v2=n*v1+v2 ) +GENERIC: n*V-in-place ( n v -- v=n*v ) + +GENERIC: V. ( v1 v2 -- v1.v2 ) +GENERIC: V.conj ( v1 v2 -- v1^H.v2 ) +GENERIC: Vnorm ( v -- norm ) +GENERIC: Vasum ( v -- abs-sum ) +GENERIC: Vswap ( v1 v2 -- v1=v2 v2=v1 ) + +GENERIC: Viamax ( v -- abs-max-index ) + +> ] [ data>> ] [ inc>> ] tri ] dip + 4 npick * + 1 ; + +MACRO: (do-copy) ( copy make-vector -- ) + '[ over 6 npick , 2dip 1 @ ] ; + +: (prepare-swap) ( v1 v2 -- length v1-data v1-inc v2-data v2-inc v1 v2 ) + [ + [ [ length>> ] bi@ min ] + [ [ [ data>> ] [ inc>> ] bi ] bi@ ] 2bi + ] 2keep ; + +: (prepare-axpy) ( n v1 v2 -- length n v1-data v1-inc v2-data v2-inc v2 ) + [ + [ [ length>> ] bi@ min swap ] + [ [ [ data>> ] [ inc>> ] bi ] bi@ ] 2bi + ] keep ; + +: (prepare-scal) ( n v -- length n v-data v-inc v ) + [ [ length>> swap ] [ data>> ] [ inc>> ] tri ] keep ; + +: (prepare-dot) ( v1 v2 -- length v1-data v1-inc v2-data v2-inc ) + [ [ length>> ] bi@ min ] + [ [ [ data>> ] [ inc>> ] bi ] bi@ ] 2bi ; + +: (prepare-nrm2) ( v -- length v1-data v1-inc ) + [ length>> ] [ data>> ] [ inc>> ] tri ; + +: (flatten-complex-sequence) ( seq -- seq' ) + [ [ real-part ] [ imaginary-part ] bi 2array ] map concat ; + +: (>c-complex) ( complex -- alien ) + [ real-part ] [ imaginary-part ] bi 2array >c-float-array ; +: (>z-complex) ( complex -- alien ) + [ real-part ] [ imaginary-part ] bi 2array >c-double-array ; + +: (c-complex>) ( alien -- complex ) + 2 c-float-array> first2 rect> ; +: (z-complex>) ( alien -- complex ) + 2 c-double-array> first2 rect> ; + +: (prepare-nth) ( n v -- n*inc v-data ) + [ inc>> ] [ data>> ] bi [ * ] dip ; + +MACRO: (complex-nth) ( nth-quot -- ) + '[ + [ 2 * dup 1+ ] dip + , curry bi@ rect> + ] ; + +: (c-complex-nth) ( n alien -- complex ) + [ float-nth ] (complex-nth) ; +: (z-complex-nth) ( n alien -- complex ) + [ double-nth ] (complex-nth) ; + +MACRO: (set-complex-nth) ( set-nth-quot -- ) + '[ + [ + [ [ real-part ] [ imaginary-part ] bi ] + [ 2 * dup 1+ ] bi* + swapd + ] dip + , curry 2bi@ + ] ; + +: (set-c-complex-nth) ( complex n alien -- ) + [ set-float-nth ] (set-complex-nth) ; +: (set-z-complex-nth) ( complex n alien -- ) + [ set-double-nth ] (set-complex-nth) ; + +PRIVATE> + +METHOD: zero-vector { float-blas-vector } + length>> 0.0 swap 0 ; +METHOD: zero-vector { double-blas-vector } + length>> 0.0 swap 0 ; +METHOD: zero-vector { float-complex-blas-vector } + length>> "CBLAS_C" swap 0 ; +METHOD: zero-vector { double-complex-blas-vector } + length>> "CBLAS_Z" swap 0 ; + +syntax:M: blas-vector-base length + length>> ; + +syntax:M: float-blas-vector nth-unsafe + (prepare-nth) float-nth ; +syntax:M: float-blas-vector set-nth-unsafe + (prepare-nth) set-float-nth ; + +syntax:M: double-blas-vector nth-unsafe + (prepare-nth) double-nth ; +syntax:M: double-blas-vector set-nth-unsafe + (prepare-nth) set-double-nth ; + +syntax:M: float-complex-blas-vector nth-unsafe + (prepare-nth) (c-complex-nth) ; +syntax:M: float-complex-blas-vector set-nth-unsafe + (prepare-nth) (set-c-complex-nth) ; + +syntax:M: double-complex-blas-vector nth-unsafe + (prepare-nth) (z-complex-nth) ; +syntax:M: double-complex-blas-vector set-nth-unsafe + (prepare-nth) (set-z-complex-nth) ; + +: >float-blas-vector ( seq -- v ) + [ >c-float-array ] [ length ] bi 1 ; +: >double-blas-vector ( seq -- v ) + [ >c-double-array ] [ length ] bi 1 ; +: >float-complex-blas-vector ( seq -- v ) + [ (flatten-complex-sequence) >c-float-array ] [ length ] bi 1 ; +: >double-complex-blas-vector ( seq -- v ) + [ (flatten-complex-sequence) >c-double-array ] [ length ] bi 1 ; + +syntax:M: float-blas-vector clone + "float" heap-size (prepare-copy) + [ cblas_scopy ] [ ] (do-copy) ; +syntax:M: double-blas-vector clone + "double" heap-size (prepare-copy) + [ cblas_dcopy ] [ ] (do-copy) ; +syntax:M: float-complex-blas-vector clone + "CBLAS_C" heap-size (prepare-copy) + [ cblas_ccopy ] [ ] (do-copy) ; +syntax:M: double-complex-blas-vector clone + "CBLAS_Z" heap-size (prepare-copy) + [ cblas_zcopy ] [ ] (do-copy) ; + +METHOD: Vswap { float-blas-vector float-blas-vector } + (prepare-swap) [ cblas_sswap ] 2dip ; +METHOD: Vswap { double-blas-vector double-blas-vector } + (prepare-swap) [ cblas_dswap ] 2dip ; +METHOD: Vswap { float-complex-blas-vector float-complex-blas-vector } + (prepare-swap) [ cblas_cswap ] 2dip ; +METHOD: Vswap { double-complex-blas-vector double-complex-blas-vector } + (prepare-swap) [ cblas_zswap ] 2dip ; + +METHOD: n*V+V-in-place { real float-blas-vector float-blas-vector } + (prepare-axpy) [ cblas_saxpy ] dip ; +METHOD: n*V+V-in-place { real double-blas-vector double-blas-vector } + (prepare-axpy) [ cblas_daxpy ] dip ; +METHOD: n*V+V-in-place { number float-complex-blas-vector float-complex-blas-vector } + [ (>c-complex) ] 2dip + (prepare-axpy) [ cblas_caxpy ] dip ; +METHOD: n*V+V-in-place { number double-complex-blas-vector double-complex-blas-vector } + [ (>z-complex) ] 2dip + (prepare-axpy) [ cblas_zaxpy ] dip ; + +METHOD: n*V-in-place { real float-blas-vector } + (prepare-scal) [ cblas_sscal ] dip ; +METHOD: n*V-in-place { real double-blas-vector } + (prepare-scal) [ cblas_dscal ] dip ; +METHOD: n*V-in-place { number float-complex-blas-vector } + [ (>c-complex) ] dip + (prepare-scal) [ cblas_cscal ] dip ; +METHOD: n*V-in-place { number double-complex-blas-vector } + [ (>z-complex) ] dip + (prepare-scal) [ cblas_zscal ] dip ; + + + +: n*V+V ( n v1 v2 -- n*v1+v2 ) clone n*V+V-in-place ; +: n*V ( n v1 -- n*v1 ) clone n*V-in-place ; + +: V+ ( v1 v2 -- v1+v2 ) + 1.0 -rot n*V+V ; +: V- ( v1 v2 -- v1+v2 ) + -1.0 spin n*V+V ; + +: Vneg ( v1 -- -v1 ) + [ zero-vector ] keep V- ; + +: V*n ( v n -- v*n ) + swap n*V ; +: V/n ( v n -- v*n ) + recip swap n*V ; + +METHOD: V. { float-blas-vector float-blas-vector } + (prepare-dot) cblas_sdot ; +METHOD: V. { double-blas-vector double-blas-vector } + (prepare-dot) cblas_ddot ; +METHOD: V. { float-complex-blas-vector float-complex-blas-vector } + (prepare-dot) + "CBLAS_C" [ cblas_cdotu_sub ] keep (c-complex>) ; +METHOD: V. { double-complex-blas-vector double-complex-blas-vector } + (prepare-dot) + "CBLAS_Z" [ cblas_zdotu_sub ] keep (z-complex>) ; + +METHOD: V.conj { float-complex-blas-vector float-complex-blas-vector } + (prepare-dot) + "CBLAS_C" [ cblas_cdotc_sub ] keep (c-complex>) ; +METHOD: V.conj { double-complex-blas-vector double-complex-blas-vector } + (prepare-dot) + "CBLAS_Z" [ cblas_zdotc_sub ] keep (z-complex>) ; + +METHOD: Vnorm { float-blas-vector } + (prepare-nrm2) cblas_snrm2 ; +METHOD: Vnorm { double-blas-vector } + (prepare-nrm2) cblas_dnrm2 ; +METHOD: Vnorm { float-complex-blas-vector } + (prepare-nrm2) cblas_scnrm2 ; +METHOD: Vnorm { double-complex-blas-vector } + (prepare-nrm2) cblas_dznrm2 ; + +METHOD: Vasum { float-blas-vector } + (prepare-nrm2) cblas_sasum ; +METHOD: Vasum { double-blas-vector } + (prepare-nrm2) cblas_dasum ; +METHOD: Vasum { float-complex-blas-vector } + (prepare-nrm2) cblas_scasum ; +METHOD: Vasum { double-complex-blas-vector } + (prepare-nrm2) cblas_dzasum ; + +METHOD: Viamax { float-blas-vector } + (prepare-nrm2) cblas_isamax ; +METHOD: Viamax { double-blas-vector } + (prepare-nrm2) cblas_idamax ; +METHOD: Viamax { float-complex-blas-vector } + (prepare-nrm2) cblas_icamax ; +METHOD: Viamax { double-complex-blas-vector } + (prepare-nrm2) cblas_izamax ; + +: Vamax ( v -- max ) + [ Viamax ] keep nth ;