From 2557a24bf7203a0d26c0c29231d21ab27d0828f7 Mon Sep 17 00:00:00 2001 From: Liam Healy <liam@thinkpad.local> Date: Sat, 9 Aug 2008 19:55:07 -0400 Subject: [PATCH] Port BLAS3 to ffa, with changes to BLAS2 Port of all BLAS3 to ffa, which also means changing BLAS2 because common generic functions are used. These compile but are not tested, and there are no regression tests/examples. --- gsll.asd | 4 +- linear-algebra/blas2.lisp | 211 ++++++++++++++++++++------- linear-algebra/blas3.lisp | 298 +++++++++++++++----------------------- 3 files changed, 279 insertions(+), 234 deletions(-) diff --git a/gsll.asd b/gsll.asd index 3c2b9032..1fc4d032 100644 --- a/gsll.asd +++ b/gsll.asd @@ -1,6 +1,6 @@ ;; Definition of GSLL system ;; Liam Healy -;; Time-stamp: <2008-08-06 22:41:57EDT gsll.asd> +;; Time-stamp: <2008-08-10 17:01:08EDT gsll.asd> ;; $Id$ (asdf:defsystem "gsll" @@ -81,7 +81,7 @@ :components ((:file "blas1") (:file "blas2") - ;(:file "blas3" :depends-on (blas2)) + (:file "blas3" :depends-on (blas2)) ;(:file "lu") ;(:file "qr") ;(:file "qrpt") diff --git a/linear-algebra/blas2.lisp b/linear-algebra/blas2.lisp index cf604f48..642ac4ed 100644 --- a/linear-algebra/blas2.lisp +++ b/linear-algebra/blas2.lisp @@ -1,6 +1,6 @@ ;; BLAS level 2, Matrix-vector operations ;; Liam Healy, Wed Apr 26 2006 - 21:08 -;; Time-stamp: <2008-08-06 22:49:49EDT blas2.lisp> +;; Time-stamp: <2008-08-09 19:43:20EDT blas2.lisp> ;; $Id$ (in-package :gsl) @@ -25,8 +25,9 @@ ;;;; Functions ;;;;**************************************************************************** -(defmfun matrix-vector-product - (TransA alpha (A matrix) (x vector) beta (y vector)) +(defmfun matrix-product + ((A matrix) (x vector) (y vector) + &optional (alpha 1) (beta 1) (TransA :notrans) TransB) ("gsl_blas_" :type "gemv") ((transa cblas-transpose) (alpha :element-c-type) ((mpointer A) :pointer) ((mpointer x) :pointer) (beta :element-c-type) ((mpointer y) :pointer)) @@ -34,32 +35,55 @@ :element-types :float-complex :inputs (A x y) :outputs (y) + :return (y) :documentation ; FDL - "The matrix-vector product and sum + "If the second and third arguments are vectors, compute + the matrix-vector product and sum y = alpha op(A) x + beta y, where op(A) = A, A^T, A^H - for TransA = :notrans, :trans, :conjtrans.") + for TransA = :notrans, :trans, :conjtrans. + If the second and third arguments are matrices, compute + the matrix-matrix product and sum C = alpha + op(A) op(B) + beta C where op(A) = A, A^T, A^H for TransA = + :notrans, :trans, :conjtrans and similarly for the + parameter TransB.") -(defmfun matrix-vector-product-triangular - (uplo TransA diag (A matrix) (x vector)) +(defmfun matrix-product-triangular + ((A matrix) (x vector) + &optional (alpha 1) (uplo :upper) (TransA :notrans) + (diag :nonunit) (side :left)) ("gsl_blas_" :type "trmv") ((uplo cblas-uplo) (TransA cblas-transpose) (diag cblas-diag) ((mpointer A) :pointer) ((mpointer x) :pointer)) :definition :generic :element-types :float-complex - :inputs (TransA A x) + :inputs (A x) :outputs (x) + :return (x) :documentation ; FDL - "The matrix-vector product x =\alpha op(A) x + "If the second argument is a vector, compute + the matrix-vector product x =alpha op(A) x for the triangular matrix A, where op(A) = A, A^T, A^H for TransA = :NoTrans, :Trans, :ConjTrans. When Uplo is :Upper then the upper triangle of A is used, and when Uplo is :Lower then the lower triangle of A is used. If Diag is :NonUnit then the diagonal of the matrix is used, but if Diag is :Unit then the diagonal elements of the - matrix A are taken as unity and are not referenced.") + matrix A are taken as unity and are not referenced. + If the second argument is a matrix, compute + the matrix-matrix product B = alpha op(A) B + if Side is :Left and B = \alpha B op(A) if Side is + :Right. The matrix A is triangular and op(A) = A, A^T, A^H + for TransA = :NoTrans, :Trans, :ConjTrans When Uplo + is :Upper then the upper triangle of A is used, and when + Uplo is :Lower then the lower triangle of A is used. If + Diag is :NonUnit then the diagonal of A is used, but if + Diag is :Unit then the diagonal elements of the matrix A + are taken as unity and are not referenced.") -(defmfun inverse-matrix-vector-product - (uplo TransA diag (A matrix) (x vector)) +(defmfun inverse-matrix-product + ((A matrix) (x vector) + &optional (alpha 1) (uplo :upper) (TransA :notrans) + (diag :nonunit) (side :left)) ("gsl_blas_" :type "trsv") ((uplo cblas-uplo) (TransA cblas-transpose) (diag cblas-diag) ((mpointer A) :pointer) ((mpointer x) :pointer)) @@ -67,17 +91,30 @@ :element-types :float-complex :inputs (A x) :outputs (x) + :return (x) :documentation ; FDL - "Compute inv(op(A)) x for x, where op(A) = A, A^T, A^H for - TransA = :NoTrans, :Trans, :ConjTrans. When Uplo is - :Upper then the upper triangle of A is used, and when Uplo is - :Lower then the lower triangle of A is used. If Diag is - :NonUnit then the diagonal of the matrix is used, but if Diag - is :Unit then the diagonal elements of the matrix A are taken - as unity and are not referenced.") - -(defmfun matrix-vector-product-symmetric - (uplo alpha (A matrix) (x vector) beta (y vector)) + "If the second argument is a vector, compute + inv(op(A)) x for x, where op(A) = A, A^T, A^H for + TransA = :NoTrans, :Trans, :ConjTrans. When Uplo is + :Upper then the upper triangle of A is used, and when Uplo is + :Lower then the lower triangle of A is used. If Diag is + :NonUnit then the diagonal of the matrix is used, but if Diag + is :Unit then the diagonal elements of the matrix A are taken + as unity and are not referenced. + If the second argument is a matrix, compute + the inverse-matrix matrix product B = alpha op(inv(A))B if + Side is :Left and B = alpha B op(inv(A)) if Side is + :Right. The matrix A is triangular and op(A) = A, A^T, A^H + for TransA = :NoTrans, :Trans, :ConjTrans When + Uplo is :Upper then the upper triangle of A is used, and + when Uplo is :Lower then the lower triangle of A is + used. If Diag is :NonUnit then the diagonal of A is used, + but if Diag is :Unit then the diagonal elements of the + matrix A are taken as unity and are not referenced.") + +(defmfun matrix-product-symmetric + ((A matrix) (x vector) (y vector) + &optional (alpha 1) (beta 1) (uplo :upper) (side :left)) ("gsl_blas_" :type "symv") ((uplo cblas-uplo) (alpha :element-c-type) ((mpointer A) :pointer) ((mpointer x) :pointer) (beta :element-c-type) ((mpointer y) :pointer)) @@ -85,16 +122,26 @@ :element-types :float :inputs (A x y) :outputs (y) + :return (y) :documentation ; FDL - "The matrix-vector product and sum y = alpha A + "If the second and third arguments are vectors, compute + the matrix-vector product and sum y = alpha A x + beta y for the symmetric matrix A. Since the matrix A is symmetric only its upper half or lower half need to be stored. When Uplo is :Upper then the upper triangle and diagonal of A are used, and when Uplo is :Lower then the - lower triangle and diagonal of A are used.") + lower triangle and diagonal of A are used. + If the second and third arguments are matrices, compute + the matrix-matrix product and sum C = alpha A + B + beta C for Side is :Left and C = alpha B A + beta C + for Side is :Right, where the matrix A is symmetric. When + Uplo is :Upper then the upper triangle and diagonal of A + are used, and when Uplo is :Lower then the lower triangle + and diagonal of A are used.") -(defmfun matrix-vector-product-hermitian - (uplo alpha (A matrix) (x vector) beta (y vector)) +(defmfun matrix-product-hermitian + ((A matrix) (x vector) (y vector) + &optional (alpha 1) (beta 1) (uplo :upper) (side :left)) ("gsl_blas_" :type "hemv") ((uplo cblas-uplo) (alpha :element-c-type) ((mpointer A) :pointer) ((mpointer x) :pointer) (beta :element-c-type) ((mpointer y) :pointer)) @@ -102,14 +149,23 @@ :element-types :complex :inputs (A x y) :outputs (y) + :return (y) :documentation ; FDL - "The matrix-vector product and sum y = alpha A x + beta y for the + "If the second and third arguments are vectors, compute the + matrix-vector product and sum y = alpha A x + beta y for the hermitian matrix A. Since the matrix A is hermitian only its upper - half or lower half need to be stored. When Uplo is CblasUpper then + half or lower half need to be stored. When Uplo is :upper then the upper triangle and diagonal of A are used, and when Uplo is - CblasLower then the lower triangle and diagonal of A are used. The + :lower then the lower triangle and diagonal of A are used. The imaginary elements of the diagonal are automatically assumed to be - zero and are not referenced.") + zero and are not referenced. If the second and third arguments are + matrices, compute the matrix-matrix product and sum C = alpha A B + + beta C if Side is :left and C = \alpha B A + \beta C if Side + is :right, where the matrix A is hermitian. When Uplo is + :upper then the upper triangle and diagonal of A are used, and + when Uplo is :lower then the lower triangle and diagonal of A + are used. The imaginary elements of the diagonal are automatically + set to zero.") (defmfun rank-1-update (alpha (x vector) (y vector) (A matrix)) ("gsl_blas_" :type "ger" :suffix) @@ -119,6 +175,7 @@ :element-types :float-complex :inputs (x y A) :outputs (A) + :return (A) :documentation ; FDL "The rank-1 update A = alpha x y^T + A of the matrix A.") @@ -130,10 +187,13 @@ :element-types :complex :inputs (x y A) :outputs (A) + :return (A) :documentation ; FDL - "The conjugate rank-1 update A = alpha x y^H + A of the matrix A.") + "The conjugate rank-1 update A = alpha x y^H + A of the matrix A.") -(defmfun symmetric-rank-1-update (uplo alpha (x vector) (A matrix)) +(defmfun symmetric-rank-1-update + ((x vector) (A matrix) + &optional (alpha 1) (beta 1) (uplo :upper) (trans :notrans)) ("gsl_blas_" :type "syr") ((uplo cblas-uplo) (alpha :element-c-type) ((mpointer x) :pointer) ((mpointer A) :pointer)) @@ -141,15 +201,25 @@ :element-types :float :inputs (x A) :outputs (A) + :return (A) :documentation ; FDL - "The symmetric rank-1 update A = \alpha x x^T + - A of the symmetric matrix A. Since the matrix A is symmetric - only its upper half or lower half need to be stored. When Uplo - is :Upper then the upper triangle and diagonal of A are - used, and when Uplo is :Lower then the lower triangle and - diagonal of A are used.") - -(defmfun hermitian-rank-1-update (uplo alpha (x vector) (A matrix)) + "If the first argument is a vector, + the symmetric rank-1 update A = \alpha x x^T + A of the symmetric + matrix A. Since the matrix A is symmetric only its upper half or + lower half need to be stored. When Uplo is :Upper then the upper + triangle and diagonal of A are used, and when Uplo is :Lower then + the lower triangle and diagonal of A are used. If the first + argument is a matrix, a rank-k update of the symmetric matrix C, C = + \alpha A A^T + \beta C when Trans is CblasNoTrans and C = \alpha A^T + A + \beta C when Trans is CblasTrans. Since the matrix C is + symmetric only its upper half or lower half need to be stored. When + Uplo is CblasUpper then the upper triangle and diagonal of C are + used, and when Uplo is CblasLower then the lower triangle and + diagonal of C are used.") + +(defmfun hermitian-rank-1-update + ((x vector) (A matrix) + &optional (alpha 1) (beta 1) (uplo :upper) (trans :notrans)) ("gsl_blas_" :type "her") ((uplo cblas-uplo) (alpha :element-c-type) ((mpointer x) :pointer) ((mpointer A) :pointer)) @@ -157,15 +227,27 @@ :element-types :complex :inputs (x A) :outputs (A) + :return (A) :documentation ; FDL - "Compute the hermitian rank-1 update A = alpha x x^H + A of the + "If the first argument is a vector, + compute the hermitian rank-1 update A = alpha x x^H + A of the hermitian matrix A. Since the matrix A is hermitian only its upper - half or lower half need to be stored. When Uplo is :upper then - the upper triangle and diagonal of A are used, and when Uplo is + half or lower half need to be stored. When Uplo is :upper then the + upper triangle and diagonal of A are used, and when Uplo is :lower then the lower triangle and diagonal of A are used. The - imaginary elements of the diagonal are automatically set to zero.") + imaginary elements of the diagonal are automatically set to zero. + If the first argument is a matrix, compute a rank-k update of the + hermitian matrix C, C = \alpha A A^H + \beta C when Trans is + :notrans and C = \alpha A^H A + \beta C when Trans is + :trans. Since the matrix C is hermitian only its upper half or + lower half need to be stored. When Uplo is :upper then the upper + triangle and diagonal of C are used, and when Uplo is :lower + then the lower triangle and diagonal of C are used. The imaginary + elements of the diagonal are automatically set to zero.") -(defmfun symmetric-rank-2-update (uplo alpha (x vector) (y vector) (A matrix)) +(defmfun symmetric-rank-2-update + ((x vector) (y vector) (A matrix) + &optional (alpha 1) (beta 1) (uplo :upper) (trans :notrans)) ("gsl_blas_" :type "syr2") ((uplo cblas-uplo) (alpha :element-c-type) ((mpointer x) :pointer) ((mpointer y) :pointer) @@ -174,15 +256,25 @@ :element-types :float :inputs (x y A) :outputs (A) + :return (A) :documentation ; FDL - "The symmetric rank-2 update A = alpha x y^T + - alpha y x^T + A of the symmetric matrix A. Since the matrix A - is symmetric only its upper half or lower half need to be - stored. When Uplo is :Upper then the upper triangle and - diagonal of A are used, and when Uplo is :Lower then the - lower triangle and diagonal of A are used.") + "If the first two arguments are vectors, compute + the symmetric rank-2 update A = alpha x y^T + alpha y x^T + A of the + symmetric matrix A. Since the matrix A is symmetric only its upper + half or lower half need to be stored. When Uplo is :upper then the + upper triangle and diagonal of A are used, and when Uplo is :lower + then the lower triangle and diagonal of A are used. If the first + two arguments are matrices, compute a rank-2k update of the + symmetric matrix C, C = \alpha A B^T + \alpha B A^T + \beta C when + Trans is :notrans and C = \alpha A^T B + \alpha B^T A + \beta C + when Trans is :trans. Since the matrix C is symmetric only its + upper half or lower half need to be stored. When Uplo is :upper + then the upper triangle and diagonal of C are used, and when Uplo is + :lower then the lower triangle and diagonal of C are used.") -(defmfun hermitian-rank-2-update (uplo alpha (x vector) (y vector) (A matrix)) +(defmfun hermitian-rank-2-update + ((x vector) (y vector) (A matrix) + &optional (alpha 1) (beta 1) (uplo :upper) (trans :notrans)) ("gsl_blas_" :type "her2") ((uplo cblas-uplo) (alpha :element-c-type) ((mpointer x) :pointer) ((mpointer A) :pointer)) @@ -190,13 +282,24 @@ :element-types :complex :inputs (x A) :outputs (A) + :return (A) :documentation ; FDL - "The hermitian rank-2 update A = alpha x y^H + alpha^* y x^H A of + "If the first two arguments are vectors, compute the + hermitian rank-2 update A = alpha x y^H + alpha^* y x^H A of the hermitian matrix A. Since the matrix A is hermitian only its upper half or lower half need to be stored. When uplo is :upper then the upper triangle and diagonal of A are used, and when uplo is :lower then the lower triangle and diagonal of A are used. The - imaginary elements of the diagonal are automatically set to zero.") + imaginary elements of the diagonal are automatically set to zero. + If the first two arguments are matrices, compute a rank-2k update of + the hermitian matrix C, C = \alpha A B^H + \alpha^* B A^H + \beta C + when Trans is :notrans and C = \alpha A^H B + \alpha^* B^H A + + \beta C when Trans is :conjtrans. Since the matrix C is + hermitian only its upper half or lower half need to be stored. When + Uplo is :upper then the upper triangle and diagonal of C are + used, and when Uplo is :lower then the lower triangle and + diagonal of C are used. The imaginary elements of the diagonal are + automatically set to zero.") ;;;;**************************************************************************** ;;;; Examples and unit test diff --git a/linear-algebra/blas3.lisp b/linear-algebra/blas3.lisp index 32aeb549..4d733a8d 100644 --- a/linear-algebra/blas3.lisp +++ b/linear-algebra/blas3.lisp @@ -1,6 +1,6 @@ ;; BLAS level 3, Matrix-matrix operations ;; Liam Healy, Wed Apr 26 2006 - 21:08 -;; Time-stamp: <2008-03-09 14:13:42EDT blas3.lisp> +;; Time-stamp: <2008-08-09 19:46:23EDT blas3.lisp> ;; $Id$ (in-package :gsl) @@ -14,185 +14,127 @@ (:left 141) :right) ;;;;**************************************************************************** -;;;; Generic +;;;; Functions ;;;;**************************************************************************** -(export '(gemm symm trmm trsm syrk syr2k)) - -(defgeneric gemm (TransA TransB alpha A B beta C) - (:method :after (TransA TransB alpha A B beta C) - (cl-invalidate C)) - (:documentation ; FDL - "The matrix-matrix product and sum C = alpha - op(A) op(B) + beta C where op(A) = A, A^T, A^H for TransA = - :NoTrans, :Trans, :ConjTrans and similarly for the - parameter TransB.")) - -(defgeneric symm (side uplo alpha A B beta C) - (:method :after (side uplo alpha A B beta C) - (cl-invalidate C)) - (:documentation ; FDL - "The matrix-matrix product and sum C = alpha A - B + beta C for Side is :Left and C = alpha B A + beta C - for Side is :Right, where the matrix A is symmetric. When - Uplo is :Upper then the upper triangle and diagonal of A - are used, and when Uplo is :Lower then the lower triangle - and diagonal of A are used.")) - -(defgeneric trmm (side uplo TransA diag alpha A B) - (:method :after (side uplo TransA diag alpha A B) - (cl-invalidate B)) - (:documentation ; FDL - "The matrix-matrix product B = \alpha op(A) B - for Side is :Left and B = \alpha B op(A) for Side is - :Right. The matrix A is triangular and op(A) = A, A^T, A^H - for TransA = :NoTrans, :Trans, :ConjTrans When Uplo - is :Upper then the upper triangle of A is used, and when - Uplo is :Lower then the lower triangle of A is used. If - Diag is :NonUnit then the diagonal of A is used, but if - Diag is :Unit then the diagonal elements of the matrix A - are taken as unity and are not referenced.")) - -(defgeneric trsm (side uplo TransA diag alpha A B) - (:method :after (side uplo TransA diag alpha A B) - (cl-invalidate B)) - (:documentation - "The inverse-matrix matrix product B = \alpha op(inv(A))B for - Side is :Left and B = \alpha B op(inv(A)) for Side is - :Right. The matrix A is triangular and op(A) = A, A^T, A^H - for TransA = :NoTrans, :Trans, :ConjTrans When - Uplo is :Upper then the upper triangle of A is used, and - when Uplo is :Lower then the lower triangle of A is - used. If Diag is :NonUnit then the diagonal of A is used, - but if Diag is :Unit then the diagonal elements of the - matrix A are taken as unity and are not referenced.")) - -(defgeneric syrk (uplo trans alpha A beta C) - (:method :after (uplo trans alpha A beta C) - (cl-invalidate C)) - (:documentation "A rank-k update of the symmetric matrix C, C = - alpha A A^T + beta C when Trans is :NoTrans and C = - alpha A^T A + beta C when Trans is :Trans. Since the - matrix C is symmetric only its upper half or lower half need to - be stored. When Uplo is :Upper then the upper triangle and - diagonal of C are used, and when Uplo is :Lower then the - lower triangle and diagonal of C are used.")) - -(defgeneric syr2k (uplo trans alpha A B beta C) - (:method :after (uplo trans alpha A B beta C) - (cl-invalidate C)) - (:documentation "A rank-2k update of the symmetric matrix C, C - = alpha A B^T + alpha B A^T + beta C when Trans is - :NoTrans and C = alpha A^T B + alpha B^T A + \beta C when - Trans is :Trans. Since the matrix C is symmetric only its - upper half or lower half need to be stored. When Uplo is - :Upper then the upper triangle and diagonal of C are used, - and when Uplo is :Lower then the lower triangle and - diagonal of C are used.")) - -;;;;**************************************************************************** -;;;; Single -;;;;**************************************************************************** - -(defmfun gemm - (TransA TransB alpha (A matrix-single-float) (B matrix-single-float) - beta (C matrix-single-float)) - "gsl_blas_sgemm" - ((transa cblas-transpose) (transb cblas-transpose) - (alpha :float) ((pointer A) gsl-matrix-c) ((pointer B) gsl-matrix-c) - (beta :float) ((pointer C) gsl-matrix-c)) - :type :method) - -(defmfun symm - (side uplo alpha (A matrix-single-float) (B matrix-single-float) - beta (C matrix-single-float)) - "gsl_blas_ssymm" - ((side cblas-side) (uplo cblas-uplo) (alpha :float) - ((pointer A) gsl-matrix-c) ((pointer B) gsl-matrix-c) - (beta :float) ((pointer C) gsl-matrix-c)) - :type :method) - -(defmfun trmm - (side uplo transa diag alpha (A matrix-single-float) (B matrix-single-float)) - "gsl_blas_strmm" - ((side cblas-side) (uplo cblas-uplo) (transa cblas-transpose) (diag cblas-diag) - (alpha :float) ((pointer A) gsl-matrix-c) ((pointer B) gsl-matrix-c)) - :type :method) - -(defmfun trsm - (side uplo transa diag alpha (A matrix-single-float) (B matrix-single-float)) - "gsl_blas_strsm" - ((side cblas-side) (uplo cblas-uplo) (transa cblas-transpose) (diag cblas-diag) - (alpha :float) ((pointer A) gsl-matrix-c) ((pointer B) gsl-matrix-c)) - :type :method) - -(defmfun syrk - (uplo trans alpha (A matrix-single-float) beta (C matrix-single-float)) - "gsl_blas_ssyrk" - ((uplo cblas-uplo) (trans cblas-transpose) (alpha :float) - ((pointer A) gsl-matrix-c) (beta :float) ((pointer C) gsl-matrix-c)) - :type :method) - -(defmfun syr2k - (uplo trans alpha (A matrix-single-float) (B matrix-single-float) - beta (C matrix-single-float)) - "gsl_blas_ssyr2k" - ((uplo cblas-uplo) (trans cblas-transpose) (alpha :float) - ((pointer A) gsl-matrix-c) ((pointer B) gsl-matrix-c) - (beta :float) ((pointer C) gsl-matrix-c)) - :type :method) - -;;;;**************************************************************************** -;;;; Double -;;;;**************************************************************************** - -(defmfun gemm - (TransA TransB alpha (A matrix-double-float) (B matrix-double-float) - beta (C matrix-double-float)) - "gsl_blas_dgemm" - ((transa cblas-transpose) (transb cblas-transpose) - (alpha :double) (A gsl-matrix-c) (B gsl-matrix-c) - (beta :double) (C gsl-matrix-c)) - :type :method) - -(defmfun symm - (side uplo alpha (A matrix-double-float) (B matrix-double-float) - beta (C matrix-double-float)) - "gsl_blas_dsymm" - ((side cblas-side) (uplo cblas-uplo) (alpha :double) - ((pointer A) gsl-matrix-c) ((pointer B) gsl-matrix-c) - (beta :double) ((pointer C) gsl-matrix-c)) - :type :method) - -(defmfun trmm - (side uplo transa diag alpha (A matrix-double-float) (B matrix-double-float)) - "gsl_blas_dtrmm" - ((side cblas-side) (uplo cblas-uplo) (transa cblas-transpose) (diag cblas-diag) - (alpha :double) ((pointer A) gsl-matrix-c) ((pointer B) gsl-matrix-c)) - :type :method) - -(defmfun trsm - (side uplo transa diag alpha (A matrix-double-float) (B matrix-double-float)) - "gsl_blas_dtrsm" - ((side cblas-side) (uplo cblas-uplo) (transa cblas-transpose) (diag cblas-diag) - (alpha :double) ((pointer A) gsl-matrix-c) ((pointer B) gsl-matrix-c)) - :type :method) - -(defmfun syrk - (uplo trans alpha (A matrix-double-float) beta (C matrix-double-float)) - "gsl_blas_dsyrk" - ((uplo cblas-uplo) (trans cblas-transpose) (alpha :double) - ((pointer A) gsl-matrix-c) (alpha :double) ((pointer C) gsl-matrix-c)) - :type :method) - -(defmfun syr2k - (uplo trans (alpha float) (A matrix-double-float) (B matrix-double-float) - beta (C matrix-double-float)) - "gsl_blas_dsyr2k" - ((uplo cblas-uplo) (trans cblas-transpose) (alpha :double) - ((pointer A) gsl-matrix-c) ((pointer B) gsl-matrix-c) - (beta :double) ((pointer C) gsl-matrix-c)) - :type :method) +(defmfun matrix-product + ((A matrix) (B matrix) (C matrix) + &optional (alpha 1) (beta 1) (TransA :notrans) (TransB :notrans)) + ("gsl_blas_" :type "gemm") + ((TransA cblas-transpose) (TransB cblas-transpose) + (alpha :element-c-type) ((mpointer A) :pointer) + ((mpointer B) :pointer) (beta :element-c-type) ((mpointer C) :pointer)) + :definition :methods + :element-types :float-complex + :inputs (A B C) + :outputs (C) + :return (C)) + +(defmfun matrix-product-symmetric + ((A matrix) (B matrix) (C matrix) + &optional (alpha 1) (beta 1) (uplo :upper) (side :left)) + ("gsl_blas_" :type "symm") + ((side cblas-side) (uplo cblas-uplo) (alpha :element-c-type) + (A :pointer) (B :pointer) (beta :element-c-type) (C :pointer)) + :definition :generic + :element-types :float-complex + :inputs (A B C) + :outputs (C) + :return (C)) + +(defmfun matrix-product-hermitian + ((A matrix) (B matrix) (C matrix) + &optional (alpha 1) (beta 1) (uplo :upper) (side :left)) + ("gsl_blas_" :type "hemm") + ((side cblas-side) (uplo cblas-uplo) (alpha :element-c-type) + ((mpointer A) :pointer) ((mpointer B) :pointer) + (beta :element-c-type) ((mpointer C) :pointer)) + :definition :methods + :element-types :complex + :inputs (A B C) + :outputs (C) + :return (C)) + +(defmfun matrix-product-triangular + ((A matrix) (B matrix) + &optional (alpha 1) (uplo :upper) (TransA :notrans) + (diag :nonunit) (side :left)) + ("gsl_blas_" :type "trmm") + ((side cblas-side) (uplo cblas-uplo) (TransA cblas-transpose) + (diag cblas-diag) + (alpha :element-c-type) ((mpointer A) :pointer) ((mpointer B) :pointer)) + :definition :methods + :element-types :float-complex + :inputs (A B) + :outputs (B) + :return (B)) + +(defmfun inverse-matrix-product + ((A matrix) (B matrix) + &optional (alpha 1) (uplo :upper) (TransA :notrans) + (diag :nonunit) (side :left)) + ("gsl_blas_" :type "trsm") + ((uplo cblas-uplo) (TransA cblas-transpose) (diag cblas-diag) + ((mpointer A) :pointer) ((mpointer B) :pointer)) + :definition :methods + :element-types :float-complex + :inputs (A B) + :outputs (B) + :return (B)) + +(defmfun symmetric-rank-1-update + ((A matrix) (C matrix) + &optional (alpha 1) (beta 1) (uplo :upper) (trans :notrans)) + ("gsl_blas_" :type "syrk") + ((uplo cblas-uplo) (trans cblas-transpose) + (alpha :element-c-type) ((mpointer A) :pointer) + (beta :element-c-type) ((mpointer C) :pointer)) + :definition :methods + :element-types :float-complex + :inputs (A C) + :outputs (C) + :return (C)) + +(defmfun hermitian-rank-1-update + ((A matrix) (C matrix) + &optional (alpha 1) (beta 1) (uplo :upper) (trans :notrans)) + ("gsl_blas_" :type "herk") + ((uplo cblas-uplo) (trans cblas-transpose) + (alpha :element-c-type) ((mpointer A) :pointer) + (beta :element-c-type) ((mpointer C) :pointer)) + :definition :methods + :element-types :complex + :inputs (A C) + :outputs (C) + :return (C)) + +(defmfun symmetric-rank-2-update + ((A matrix) (B matrix) (C matrix) + &optional (alpha 1) (beta 1) (uplo :upper) (trans :notrans)) + ("gsl_blas_" :type "syr2k") + ((uplo cblas-uplo) (trans cblas-transpose) + (alpha :element-c-type) ((mpointer A) :pointer) + ((mpointer B) :pointer) (beta :element-c-type) + ((mpointer C) :pointer)) + :definition :methods + :element-types :float + :inputs (A B C) + :outputs (C) + :return (C)) + +(defmfun hermitian-rank-2-update + ((A matrix) (B matrix) (C matrix) + &optional (alpha 1) (beta 1) (uplo :upper) (trans :notrans)) + ("gsl_blas_" :type "her2k") + ((uplo cblas-uplo) (trans cblas-transpose) + (alpha :element-c-type) ((mpointer A) :pointer) + ((mpointer B) :pointer) (beta :element-c-type) + ((mpointer C) :pointer)) + :definition :methods + :element-types :complex + :inputs (A B C) + :outputs (C) + :return (C)) ;;;;**************************************************************************** ;;;; Examples and unit test -- GitLab