blas3.lisp 7.41 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: <2010-07-07 14:25:00EDT blas3.lisp>
4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
;;
;; Copyright 2006, 2007, 2008, 2009 Liam M. Healy
;; Distributed under the terms of the GNU General Public License
;;
;; This program is free software: you can redistribute it and/or modify
;; it under the terms of the GNU General Public License as published by
;; the Free Software Foundation, either version 3 of the License, or
;; (at your option) any later version.
;;
;; This program is distributed in the hope that it will be useful,
;; but WITHOUT ANY WARRANTY; without even the implied warranty of
;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
;; GNU General Public License for more details.
;;
;; You should have received a copy of the GNU General Public License
;; along with this program.  If not, see <http://www.gnu.org/licenses/>.
20 21 22

(in-package :gsl)

23 24
;;; /usr/include/gsl/gsl_blas.h

25 26 27 28
;;;;****************************************************************************
;;;; Options
;;;;****************************************************************************

Liam Healy's avatar
Liam Healy committed
29
#+fsbv
30 31
(fsbv:defcenum-aux cblas-side)

32
;;;;****************************************************************************
33
;;;; Functions
34 35
;;;;****************************************************************************

36
(defmfun matrix-product
37 38
    ((A matrix) (B matrix)
     &optional
39 40 41 42 43
     C
     (alpha 1) (beta 1) (TransA :notrans) (TransB :notrans)
     &aux
     (Carr
      (or C
44
	  (grid:make-foreign-array element-type :dimensions (matrix-product-dimensions A B)
45
		       :initial-element 0))))
46 47 48
  ("gsl_blas_" :type "gemm")
  ((TransA cblas-transpose) (TransB cblas-transpose)
   (alpha :element-c-type) ((mpointer A) :pointer)
49
   ((mpointer B) :pointer) (beta :element-c-type) ((mpointer Carr) :pointer))
50
  :definition :methods
Liam Healy's avatar
Liam Healy committed
51
  :element-types #+fsbv :float-complex #-fsbv :float
52 53
  :inputs (A B Carr)
  :outputs (Carr))
54 55

(defmfun matrix-product-symmetric
56
    ((A matrix) (B matrix)
57 58 59 60
     &optional C (alpha 1) (beta 1) (uplo :upper) (side :left)
     &aux
     (Carr
      (or C
61
	  (grid:make-foreign-array element-type :dimensions (matrix-product-dimensions A B)
62
		       :initial-element 0))))
63 64
  ("gsl_blas_" :type "symm")
  ((side cblas-side) (uplo cblas-uplo) (alpha :element-c-type)
Liam Healy's avatar
Liam Healy committed
65
   ((mpointer A) :pointer) ((mpointer B) :pointer)
66
   (beta :element-c-type) ((mpointer Carr) :pointer)) 
Liam Healy's avatar
Liam Healy committed
67
  :definition :methods
Liam Healy's avatar
Liam Healy committed
68
  :element-types #+fsbv :float-complex #-fsbv :float
69 70
  :inputs (A B Carr)
  :outputs (Carr))
71

Liam Healy's avatar
Liam Healy committed
72
#+fsbv
73
(defmfun matrix-product-hermitian
74 75
    ((A matrix) (B matrix)
     &optional
76
     (C (grid:make-foreign-array element-type :dimensions (matrix-product-dimensions A B)
77
		     :initial-element 0))
78
     (alpha 1) (beta 1) (uplo :upper) (side :left))
Liam Healy's avatar
Liam Healy committed
79 80
  ;; This always signals an error because you can't pass a
  ;; struct in CFFI yet.
81 82 83 84 85 86 87
  ("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)
88
  :outputs (C))
89 90 91 92 93 94 95 96 97 98

(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
Liam Healy's avatar
Liam Healy committed
99
  :element-types #+fsbv :float-complex #-fsbv :float
100
  :inputs (A B)
101
  :outputs (B))
102 103 104 105 106

(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
107 108
  ;; This signals an error for complex arguments because you can't pass a
  ;; struct in CFFI yet.
109
  ("gsl_blas_" :type "trsm")
Liam Healy's avatar
Liam Healy committed
110 111 112
  ((side cblas-side) (uplo cblas-uplo)
   (TransA cblas-transpose) (diag cblas-diag)
   (alpha :element-c-type) ((mpointer A) :pointer) ((mpointer B) :pointer))
113
  :definition :methods
Liam Healy's avatar
Liam Healy committed
114
  :element-types #+fsbv :float-complex #-fsbv :float
115
  :inputs (A B)
116
  :outputs (B))
117 118 119 120 121 122 123 124 125

(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
Liam Healy's avatar
Liam Healy committed
126
  :element-types #+fsbv :float-complex #-fsbv :float
127
  :inputs (A C)
128
  :outputs (C))
129

Liam Healy's avatar
Liam Healy committed
130
#+fsbv
131 132 133 134 135 136 137 138 139 140
(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)
141
  :outputs (C))
142 143 144 145 146 147 148 149 150

(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))
151 152
  :definition :methods 
  :element-types #+fsbv :float-complex #-fsbv :float
153
  :inputs (A B C)
154
  :outputs (C))
155

Liam Healy's avatar
Liam Healy committed
156
#+fsbv
157 158 159 160 161 162 163 164 165 166 167
(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)
168
  :outputs (C))
169 170 171 172 173

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

Liam Healy's avatar
Liam Healy committed
174
(generate-all-array-tests matrix-product #+fsbv :float-complex #-fsbv :float
175
 (let ((m1 (array-default '(3 3)))
176 177 178 179
       (m2 (array-default '(3 3)))
       (m3 (array-default '(3 3)))
       (s1 (scalar-default))
       (s2 (scalar-default)))
180
   (grid:copy-to (matrix-product m1 m2 m3 s1 s2))))
Liam Healy's avatar
Liam Healy committed
181

182 183 184
(generate-all-array-tests matrix-product-nonsquare :float-complex
 (let ((m1 (array-default '(2 3)))
       (m2 (array-default '(3 2))))
185
   (grid:copy-to (matrix-product m1 m2))))
186

Liam Healy's avatar
Liam Healy committed
187 188
(generate-all-array-tests matrix-product-triangular
			  #+fsbv :float-complex #-fsbv :float
189
 (let ((m1 (array-default '(3 3)))
190 191
       (m2 (array-default '(3 3)))
       (s1 (scalar-default)))
192
   (grid:copy-to (matrix-product-triangular m1 m2 s1))))
Liam Healy's avatar
Liam Healy committed
193

Liam Healy's avatar
Liam Healy committed
194 195
(generate-all-array-tests inverse-matrix-product 
			  #+fsbv :float-complex #-fsbv :float
196
 (let ((m1 (array-default '(3 3)))
197 198
       (m2 (array-default '(3 3)))
       (s1 (scalar-default)))
199
   (grid:copy-to (inverse-matrix-product m1 m2 s1))))
Liam Healy's avatar
Liam Healy committed
200 201

(generate-all-array-tests matrix-product-symmetric :float
202
 (let ((m1 (array-default '(3 3)))
203 204 205 206
       (m2 (array-default '(3 3)))
       (m3 (array-default '(3 3)))
       (s1 (scalar-default))
       (s2 (scalar-default)))
207
   (grid:copy-to (matrix-product-symmetric m1 m2 m3 s1 s2))))
Liam Healy's avatar
Liam Healy committed
208

Liam Healy's avatar
Liam Healy committed
209
#+fsbv
Liam Healy's avatar
Liam Healy committed
210
(generate-all-array-tests matrix-product-hermitian :complex
211
 (let ((m1 (array-default '(3 3)))
212 213 214 215
       (m2 (array-default '(3 3)))
       (m3 (array-default '(3 3)))
       (s1 (scalar-default))
       (s2 (scalar-default)))
216
   (grid:copy-to (matrix-product-hermitian m1 m2 m3 s1 s2))))