blas3.lisp 5.78 KB
Newer Older
liam's avatar
liam committed
1
2
;; BLAS level 3, Matrix-matrix operations
;; Liam Healy, Wed Apr 26 2006 - 21:08
3
;; Time-stamp: <2008-12-07 18:27:01EST blas3.lisp>
lhealy's avatar
lhealy committed
4
;; $Id$
5
6
7
8
9
10
11
12
13
14
15
16

(in-package :gsl)

;;;;****************************************************************************
;;;; Options
;;;;****************************************************************************

(cffi:defcenum cblas-side
  "/usr/include/gsl/gsl_cblas.h."
  (:left 141) :right)

;;;;****************************************************************************
17
;;;; Functions
18
19
;;;;****************************************************************************

20
21
22
23
24
25
26
27
28
29
(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)
30
  :outputs (C))
31
32
33
34
35
36

(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)
Liam Healy's avatar
Liam Healy committed
37
38
   ((mpointer A) :pointer) ((mpointer B) :pointer)
   (beta :element-c-type) ((mpointer C) :pointer)) 
Liam Healy's avatar
Liam Healy committed
39
  :definition :methods
40
41
  :element-types :float-complex
  :inputs (A B C)
42
  :outputs (C))
43
44
45
46

(defmfun matrix-product-hermitian
    ((A matrix) (B matrix) (C matrix)
     &optional (alpha 1) (beta 1) (uplo :upper) (side :left))
Liam Healy's avatar
Liam Healy committed
47
48
  ;; This always signals an error because you can't pass a
  ;; struct in CFFI yet.
49
50
51
52
53
54
55
  ("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)
56
  :outputs (C))
57
58
59
60
61
62
63
64
65
66
67
68

(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)
69
  :outputs (B))
70
71
72
73
74

(defmfun inverse-matrix-product
    ((A matrix) (B matrix)
     &optional (alpha 1) (uplo :upper) (TransA :notrans)
     (diag :nonunit) (side :left))
Liam Healy's avatar
Liam Healy committed
75
76
  ;; This signals an error for complex arguments because you can't pass a
  ;; struct in CFFI yet.
77
  ("gsl_blas_" :type "trsm")
Liam Healy's avatar
Liam Healy committed
78
79
80
  ((side cblas-side) (uplo cblas-uplo)
   (TransA cblas-transpose) (diag cblas-diag)
   (alpha :element-c-type) ((mpointer A) :pointer) ((mpointer B) :pointer))
81
82
83
  :definition :methods
  :element-types :float-complex
  :inputs (A B)
84
  :outputs (B))
85
86
87
88
89
90
91
92
93
94
95

(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)
96
  :outputs (C))
97
98
99
100
101
102
103
104
105
106
107

(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)
108
  :outputs (C))
109
110
111
112
113
114
115
116
117
118
119
120

(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)
121
  :outputs (C))
122
123
124
125
126
127
128
129
130
131
132
133

(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)
134
  :outputs (C))
135
136
137
138
139

;;;;****************************************************************************
;;;; Examples and unit test
;;;;****************************************************************************

Liam Healy's avatar
Liam Healy committed
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
(generate-all-array-tests matrix-product :float-complex
 (letm ((m1 (array-default '(3 3)))
	(m2 (array-default '(3 3)))
	(answer (array-default '(3 3) t))
	(s1 (scalar-default))
	(s2 (scalar-default)))
   (cl-array (matrix-product m1 m2 answer s1 s2))))

(generate-all-array-tests matrix-product-triangular :float-complex
 (letm ((m1 (array-default '(3 3)))
	(m2 (array-default '(3 3)))
	(s1 (scalar-default)))
   (cl-array (matrix-product-triangular m1 m2 s1))))

(generate-all-array-tests inverse-matrix-product :float-complex
 (letm ((m1 (array-default '(3 3)))
	(m2 (array-default '(3 3)))
	(s1 (scalar-default)))
   (cl-array (inverse-matrix-product m1 m2 s1))))

(generate-all-array-tests matrix-product-symmetric :float
 (letm ((m1 (array-default '(3 3)))
	(m2 (array-default '(3 3)))
	(answer (array-default '(3 3) t))
	(s1 (scalar-default))
	(s2 (scalar-default)))
   (cl-array (matrix-product-symmetric m1 m2 answer s1 s2))))

(generate-all-array-tests matrix-product-hermitian :complex
 (letm ((m1 (array-default '(3 3)))
	(m2 (array-default '(3 3)))
	(answer (array-default '(3 3) t))
	(s1 (scalar-default))
	(s2 (scalar-default)))
   (cl-array (matrix-product-hermitian m1 m2 answer s1 s2))))