Commit 3f52993b authored by liam's avatar liam

Add blas2 and blas3, untested, but they do compile.

git-svn-id: svn+ssh://pop/opt/space/mathematics/gsl/trunk@3050 a3d8a0fb-c1db-0310-ace7-a616afeb9e30
parent a82143c7
;********************************************************
; file: blas1.lisp
; description: BLAS level 1
; description: BLAS level 1, Vector operations
; date: Wed Apr 26 2006 - 15:23
; author: Liam Healy
; modified: Wed Apr 26 2006 - 17:36
; modified: Wed Apr 26 2006 - 21:54
;********************************************************
;;; $Id: $
......@@ -52,7 +52,7 @@
"Apply a Givens rotation (x', y') = (c x + s y, -s x + c y) to the vectors x, y."))
;;;;****************************************************************************
;;;; Double vectors
;;;; Double
;;;;****************************************************************************
(defun-gsl dot (((pointer vec1) :pointer) ((pointer vec2) :pointer))
......
;********************************************************
; file: blas2.lisp
; description: BLAS level 2, Matrix-vector operations
; date: Wed Apr 26 2006 - 21:08
; author: Liam M. Healy
; modified: Wed Apr 26 2006 - 22:43
;********************************************************
;;; $Id: $
(in-package :gsl)
;;;;****************************************************************************
;;;; Options
;;;;****************************************************************************
(cffi:defcenum cblas-transpose
"CBLAS_TRANSPOSE from /usr/include/gsl/gsl_cblas.h."
(:notrans 111) :trans :conjtrans)
(cffi:defcenum cblas-uplo
"/usr/include/gsl/gsl_cblas.h."
(:upper 121) :lower)
(cffi:defcenum cblas-diag
"/usr/include/gsl/gsl_cblas.h."
(:nonunit 131) :unit)
;;;;****************************************************************************
;;;; Generic
;;;;****************************************************************************
(export '(gemv trmv trsv symv ger syr syr2))
(defgeneric gemv (TransA alpha A x beta y)
(:documentation "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.")
(:method :after (TransA alpha A x beta y)
(cl-invalidate y)))
(defgeneric trmv (uplo TransA diag A x)
(:documentation "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.")
(:method :after (uplo TransA diag A x)
(cl-invalidate x)))
(defgeneric trsv (uplo TransA diag A x)
(:documentation
"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.")
(:method :after (uplo TransA diag A x)
(cl-invalidate x)))
(defgeneric symv (uplo alpha A x beta y)
(:documentation "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.")
(:method :after (uplo alpha A x beta y)
(cl-invalidate y)))
(defgeneric ger (alpha x y A)
(:documentation "The rank-1 update A = \alpha x y^T + A of the matrix A.")
(:method :after (alpha x y A)
(cl-invalidate A)))
(defgeneric syr (uplo alpha x A)
(:documentation "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.")
(:method :after (uplo alpha x A)
(cl-invalidate A)))
(defgeneric syr2 (uplo alpha x y A)
(:documentation "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.")
(:method :after (uplo alpha x y A)
(cl-invalidate A)))
;;;;****************************************************************************
;;;; Double
;;;;****************************************************************************
(defun-gsl gemv
((transa cblas-transpose) (alpha :double) (A gsl-matrix-c)
(x gsl-vector-c) (beta :double) (y gsl-vector-c))
"gsl_blas_dgemv"
:method
(TransA alpha (A gsl-matrix) (x gsl-vector) beta (y gsl-vector)))
(defun-gsl trmv
((uplo cblas-uplo) (transa cblas-transpose) (diag cblas-diag)
(A gsl-matrix-c) (x gsl-vector-c))
"gsl_blas_dtrmv"
:method
(uplo transa diag (A gsl-matrix) (x gsl-vector)))
(defun-gsl trsv
((uplo cblas-uplo) (transa cblas-transpose) (diag cblas-diag)
(A gsl-matrix-c) (x gsl-vector-c))
"gsl_blas_dtrsv"
:method
(uplo transa diag (A gsl-matrix) (x gsl-vector)))
(defun-gsl symv
((uplo cblas-uplo) (alpha :double) (A gsl-matrix-c)
(x gsl-vector-c) (beta :double) (y gsl-vector-c))
"gsl_blas_dsymv"
:method
(uplo alpha (A gsl-matrix) (x gsl-vector) beta (y gsl-vector)))
(defun-gsl ger
((alpha :double) (x gsl-vector-c) (y gsl-vector-c) (A gsl-matrix-c))
"gsl_blas_dger"
:method
(alpha (x gsl-vector) (y gsl-vector) (A gsl-matrix)))
(defun-gsl syr
((uplo cblas-uplo) (alpha :double) (x gsl-vector-c) (A gsl-matrix-c))
"gsl_blas_dsyr"
:method
(uplo alpha (x gsl-vector) (A gsl-matrix)))
(defun-gsl syr2
((uplo cblas-uplo) (alpha :double)
(x gsl-vector-c) (y gsl-vector-c) (A gsl-matrix-c))
"gsl_blas_dsyr2"
:method
(uplo alpha (x gsl-vector) (y gsl-vector) (A gsl-matrix)))
;;;;****************************************************************************
;;;; Examples and unit test
;;;;****************************************************************************
;********************************************************
; file: blas3.lisp
; description: BLAS level 3, Matrix-matrix operations
; date: Wed Apr 26 2006 - 21:08
; author: Liam M. Healy
; modified: Wed Apr 26 2006 - 22:47
;********************************************************
;;; $Id: $
(in-package :gsl)
;;;;****************************************************************************
;;;; Options
;;;;****************************************************************************
(cffi:defcenum cblas-side
"/usr/include/gsl/gsl_cblas.h."
(:left 141) :right)
;;;;****************************************************************************
;;;; Generic
;;;;****************************************************************************
(export '(gemm symm trmm trsm syrk syr2k))
(defgeneric gemm (TransA TransB alpha A B beta C)
(:documentation "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.")
(:method :after (TransA TransB alpha A B beta C)
(cl-invalidate C)))
(defgeneric symm (side uplo alpha A B beta C)
(:documentation "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.")
(:method :after (side uplo alpha A B beta C)
(cl-invalidate C)))
(defgeneric trmm (side uplo TransA diag alpha A B)
(:documentation "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.")
(:method :after (side uplo TransA diag alpha A B)
(cl-invalidate B)))
(defgeneric trsm (side uplo TransA diag alpha A 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.")
(:method :after (side uplo TransA diag alpha A B)
(cl-invalidate B)))
(defgeneric syrk (uplo trans alpha A beta 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.")
(:method :after (uplo trans alpha A beta C)
(cl-invalidate C)))
(defgeneric syr2k (uplo trans alpha A B beta 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.")
(:method :after (uplo trans alpha A B beta C)
(cl-invalidate C)))
;;;;****************************************************************************
;;;; Double
;;;;****************************************************************************
(defun-gsl gemm
((transa cblas-transpose) (transb cblas-transpose)
(alpha :double) (A gsl-matrix-c) (B gsl-matrix-c)
(beta :double) (C gsl-matrix-c))
"gsl_blas_dgemm"
:method
(TransA TransB alpha (A gsl-matrix) (B gsl-matrix) beta (C gsl-matrix)))
(defun-gsl symm
((side cblas-side) (uplo cblas-uplo) (alpha :double)
(A gsl-matrix-c) (B gsl-matrix-c) (beta :double) (C gsl-matrix-c))
"gsl_blas_dsymm"
:method
(side uplo alpha (A gsl-matrix) (B gsl-matrix) beta (C gsl-matrix)))
(defun-gsl trmm
((side cblas-side) (uplo cblas-uplo) (transa cblas-transpose) (diag cblas-diag)
(alpha :double) (A gsl-matrix-c) (B gsl-matrix-c))
"gsl_blas_dtrmm"
:method
(side uplo transa diag alpha (A gsl-matrix) (B gsl-matrix)))
(defun-gsl trsm
((side cblas-side) (uplo cblas-uplo) (transa cblas-transpose) (diag cblas-diag)
(alpha :double) (A gsl-matrix-c) (B gsl-matrix-c))
"gsl_blas_dtrsm"
:method
(side uplo transa diag alpha (A gsl-matrix) (B gsl-matrix)))
(defun-gsl syrk
((uplo cblas-uplo) (trans cblas-transpose) (alpha :double)
(A gsl-matrix-c) (alpha :double) (C gsl-matrix-c))
"gsl_blas_dsyrk"
:method
(uplo trans alpha (A gsl-matrix) beta (C gsl-matrix)))
(defun-gsl syr2k
((uplo cblas-uplo) (trans cblas-transpose) (alpha :double)
(A gsl-matrix-c) (B gsl-matrix-c) (beta :double) (C gsl-matrix-c))
"gsl_blas_dsyr2k"
:method
(uplo trans alpha (A gsl-matrix) (B gsl-matrix) beta (C gsl-matrix)))
;;;;****************************************************************************
;;;; Examples and unit test
;;;;****************************************************************************
......@@ -3,7 +3,7 @@
; description: Definition of GSLL system
; date:
; author: Liam Healy
; modified: Wed Apr 26 2006 - 16:20
; modified: Wed Apr 26 2006 - 22:31
;********************************************************
;;; $Id: $
......@@ -13,7 +13,7 @@
:version "0"
:author "Liam M. Healy"
:licence "GPL"
:depends-on (cffi cffi-unix) ; http://www.cliki.net/cffi-unix
:depends-on (cffi cffi-unix) ; http://www.cliki.net/cffi-unix
:components
((:module init
:components
......@@ -21,7 +21,7 @@
(:file "interface" :depends-on (init))
;; http://www.cs.northwestern.edu/academics/courses/325/readings/lisp-unit.html
(:file "lisp-unit")
(:file "tests")))
(:file "tests" :depends-on (lisp-unit))))
(:module general
:depends-on (init)
:components
......@@ -29,7 +29,7 @@
(:file "mathematical")))
;; complex numbers not necessary? Just make a struct.
(:file "cffi-array")
(:file "polynomial" :depends-on (init cffi-array)) ; see file
(:file "polynomial" :depends-on (init cffi-array)) ; see file
(:module data
:depends-on (init)
:components
......@@ -39,7 +39,9 @@
(:file "matrix" :depends-on (data))
(:file "permutation" :depends-on (data))
(:file "combination" :depends-on (data))
(:file "blas1")))
(:file "blas1" :depends-on (vector))
(:file "blas2" :depends-on (vector matrix))
(:file "blas3" :depends-on (blas2 matrix))))
(:module special-functions
:depends-on (init)
:components
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment