blas1.lisp 8.75 KB
Newer Older
1 2
;; BLAS level 1, Vector operations
;; Liam Healy, Wed Apr 26 2006 - 15:23
3
;; Time-stamp: <2010-07-07 14:25:01EDT blas1.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
(defmfun dot ((vec1 vector) (vec2 vector))
lhealy's avatar
lhealy committed
26
  ("gsl_blas_" :type "dot" :suffix)
27
  (((mpointer vec1) :pointer) ((mpointer vec2) :pointer)
28
   (result (:pointer :element-c-type)))
29 30
  :definition :generic
  :element-types :float-complex
31
  :inputs (vec1 vec2)
lhealy's avatar
lhealy committed
32 33
  :documentation			; FDL
  "Dot, or inner, product between vectors.")
34 35 36 37 38 39

;;; gsl_blas_sdsdot doesn't make much sense, but here it is.
(defmfun sdot (result alpha vec1 vec2)
  "gsl_blas_sdsdot"
  ((alpha :float) ((mpointer vec1) :pointer) ((mpointer vec2) :pointer)
   (result :pointer))
40
  :inputs (vec1 vec2)
41
  :outputs (result)
lhealy's avatar
lhealy committed
42 43
  :documentation			; FDL
  "Sum of a scalar and a dot product for single-floats.")
44 45 46

;;; Not porting gsl_blas_dsdot, stupid.

47
(defmfun cdot ((x vector) (y vector))
lhealy's avatar
lhealy committed
48
  ("gsl_blas_" :type "dotc")
49
  (((mpointer x) :pointer) ((mpointer y) :pointer)
50
   (result (:pointer :element-c-type)))
51 52 53
  :definition :generic
  :element-types :complex
  :inputs (x y)
lhealy's avatar
lhealy committed
54
  :documentation			; FDL
55 56
  "The complex conjugate scalar product x^H y for the vectors.")

57
(defmfun euclidean-norm ((vec vector))
lhealy's avatar
lhealy committed
58
  ("gsl_blas_" :component-float-type :type "nrm2")
59
  (((mpointer vec) :pointer))
60 61
  :definition :generic
  :element-types :float-complex
62 63
  :inputs (vec)
  :c-return :component-float-type
lhealy's avatar
lhealy committed
64
  :documentation			; FDL
65
  "The Euclidean norm ||x||_2 = \sqrt {\sum x_i^2} of the vector x.")
66

lhealy's avatar
lhealy committed
67 68
(defmfun absolute-sum ((x vector))
  ("gsl_blas_" :component-float-type :type "asum")
Liam Healy's avatar
Liam Healy committed
69
  (((mpointer x) :pointer))
lhealy's avatar
lhealy committed
70 71
  :definition :generic
  :element-types :float-complex
Liam Healy's avatar
Liam Healy committed
72
  :c-return :component-float-type
lhealy's avatar
lhealy committed
73 74 75 76 77 78
  :inputs (x)
  :documentation			; FDL
  "The absolute sum \sum |x_i| of the elements of the vector x.")

(defmfun index-max ((vec vector))
  ("gsl_blas_i" :type "amax")
Liam Healy's avatar
Liam Healy committed
79
  (((mpointer vec) :pointer))
lhealy's avatar
lhealy committed
80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111
  :definition :generic
  :element-types :float-complex
  :c-return sizet
  :inputs (vec)
  :documentation			; FDL
  "The index of the largest element of the vector
   x. The largest element is determined by its absolute magnitude for
   real vectors and by the sum of the magnitudes of the real and
   imaginary parts |\Re(x_i)| + |\Im(x_i)| for complex vectors. If the
   largest value occurs several times then the index of the first
   occurrence is returned.")

(defmfun blas-swap ((vec1 vector) (vec2 vector))
  ("gsl_blas_" :type "swap")
  (((mpointer vec1) :pointer) ((mpointer vec2) :pointer))
  :definition :generic
  :element-types :float-complex
  :inputs (vec1 vec2)
  :outputs (vec1 vec2)
  :documentation			; FDL
  "Exchange the elements of the vectors.")

(defmfun blas-copy ((x vector) (y vector))
  ("gsl_blas_" :type "copy")
  (((mpointer x) :pointer) ((mpointer y) :pointer))
  :definition :generic
  :element-types :float-complex
  :inputs (x)
  :outputs (y)
  :documentation			; FDL
  "Copy the elements of the vector x into the vector y.")

112 113 114
(defmfun axpy
    (alpha (x vector)
	   &optional
115
	   (y (grid:make-foreign-array element-type :dimensions (dimensions x))))
Liam Healy's avatar
Liam Healy committed
116 117
  ;; This gets an error for complex types because you can't pass a
  ;; struct in CFFI yet.
lhealy's avatar
lhealy committed
118 119 120
  ("gsl_blas_" :type "axpy")
  ((alpha  :element-c-type) ((mpointer x) :pointer) ((mpointer y) :pointer))
  :definition :generic
Liam Healy's avatar
Liam Healy committed
121
  :element-types #+fsbv :float-complex #-fsbv :float
lhealy's avatar
lhealy committed
122 123 124
  :inputs (x y)
  :outputs (y)
  :documentation			; FDL
Liam Healy's avatar
Liam Healy committed
125
  "Compute the sum y = \alpha x + y for the vectors x and y.")
lhealy's avatar
lhealy committed
126

Liam Healy's avatar
Liam Healy committed
127
(defmfun scale ((alpha :element-type) (x vector))
128 129
  ;; Alpha is the same type as the elements of vector, so for complex
  ;; vectors it must be complex.
lhealy's avatar
lhealy committed
130
  ("gsl_blas_" :type "scal")
131
  ((alpha :element-c-type) ((mpointer x) :pointer))
lhealy's avatar
lhealy committed
132
  :definition :generic
Liam Healy's avatar
Liam Healy committed
133
  :element-types #+fsbv :float-complex #-fsbv :float
Liam Healy's avatar
Liam Healy committed
134
  :c-return :void
lhealy's avatar
lhealy committed
135 136 137 138 139
  :inputs (x)
  :outputs (x)
  :documentation			; FDL
  "Rescale the vector x by the multiplicative factor alpha.")

140 141
(defmfun scale ((alpha :component-float-type) (x vector))
  ;; Alpha is a float and the vector is complex
lhealy's avatar
lhealy committed
142 143 144 145
  ("gsl_blas_" :type :component-float-type "scal")
  ((alpha :component-float-type) ((mpointer x) :pointer))
  :definition :methods
  :element-types :complex
Liam Healy's avatar
Liam Healy committed
146
  :c-return :void
lhealy's avatar
lhealy committed
147
  :inputs (x)
148
  :outputs (x))
lhealy's avatar
lhealy committed
149

Liam Healy's avatar
Liam Healy committed
150 151 152
;;; The Givens rotations come in two forms, those that work on bare C
;;; arrays, and those that work on GSL vectors.  Ports of both are
;;; present.
lhealy's avatar
lhealy committed
153

Liam Healy's avatar
Liam Healy committed
154 155
(defmfun givens-rotation ((x vector) (y vector) (c vector) (s vector))
  ("gsl_blas_" :type "rotg")
156 157
  (((foreign-pointer x) :pointer) ((foreign-pointer y) :pointer)
   ((foreign-pointer c) :pointer) ((foreign-pointer s) :pointer))
Liam Healy's avatar
Liam Healy committed
158 159 160 161 162 163 164 165 166
  :definition :generic
  :element-types :float
  :inputs (x y c s)
  :outputs (x y)
  :documentation			; FDL
  "These functions compute a Givens rotation (c,s) to the vector (x,y),
          [  c  s ] [ x ] = [ r ]
          [ -s  c ] [ y ]   [ 0 ]
   The variables x and y are overwritten by the routine.")
167

Liam Healy's avatar
Liam Healy committed
168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184
(defmfun givens-rotation-m ((x vector) (y vector) (c vector) (s vector))
  ("gsl_blas_" :type "rot")
  (((mpointer x) :pointer) ((mpointer y) :pointer)
   ((mpointer c) :pointer) ((mpointer s) :pointer))
  :definition :generic
  :element-types :float
  :inputs (x y c s)
  :outputs (x y)
  :documentation			; FDL
  "These functions compute a Givens rotation (c,s) to the vector (x,y),
          [  c  s ] [ x ] = [ r ]
          [ -s  c ] [ y ]   [ 0 ]
   The variables x and y are overwritten by the routine.")

(defmfun modified-givens-rotation
    ((d1 vector) (d2 vector) (b1 vector) b2 (P vector))
  ("gsl_blas_" :type "rotmg")
185 186 187
  (((foreign-pointer d1) :pointer) ((foreign-pointer d2) :pointer)
   ((foreign-pointer b1) :pointer) (b2 :element-c-type)
   ((foreign-pointer P) :pointer))
Liam Healy's avatar
Liam Healy committed
188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204
  :definition :generic
  :element-types :float
  :inputs (d1 d2 b1 P)
  :outputs ()				; I have no idea
  :documentation
  "Not explained")

(defmfun modified-givens-rotation-m ((x vector) (y vector) (P vector))
  ("gsl_blas_" :type "rotm")
  (((mpointer x) :pointer) ((mpointer y) :pointer)
   ((mpointer P) :pointer))
  :definition :generic
  :element-types :float
  :inputs (x y P)
  :outputs (x y)			;?????
  :documentation	
  "Not explained")
205 206 207 208 209

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

Liam Healy's avatar
Liam Healy committed
210
(generate-all-array-tests dot :float-complex
211
 (let ((v1 (array-default 8))
Liam Healy's avatar
Liam Healy committed
212
       (v2 (array-default 8)))
Liam Healy's avatar
Liam Healy committed
213 214 215
   (dot v1 v2)))

(generate-all-array-tests cdot :complex
216
 (let ((v1 (array-default 8))
Liam Healy's avatar
Liam Healy committed
217
       (v2 (array-default 8)))
Liam Healy's avatar
Liam Healy committed
218 219 220
   (cdot v1 v2)))

(generate-all-array-tests euclidean-norm :float-complex
221
 (let ((v1 (array-default 8)))
Liam Healy's avatar
Liam Healy committed
222 223 224
   (euclidean-norm v1)))

(generate-all-array-tests absolute-sum :float-complex
225
 (let ((v1 (array-default 8)))
Liam Healy's avatar
Liam Healy committed
226 227 228
   (absolute-sum v1)))

(generate-all-array-tests index-max :float-complex
229
 (let ((v1 (array-default 8)))
Liam Healy's avatar
Liam Healy committed
230 231 232
   (index-max v1)))

(generate-all-array-tests blas-swap :float-complex
233
 (let ((v1 (array-default 3))
Liam Healy's avatar
Liam Healy committed
234 235
	(v2 (array-default 3)))
   (blas-swap v2 v1)
236
   (list (grid:copy-to v1) (grid:copy-to v2))))
Liam Healy's avatar
Liam Healy committed
237 238

(generate-all-array-tests blas-copy :float-complex
239
 (let ((v1 (array-default 3))
Liam Healy's avatar
Liam Healy committed
240 241
	(v2 (array-default 3 t)))
   (blas-copy v1 v2)
242
   (grid:copy-to v2)))
Liam Healy's avatar
Liam Healy committed
243 244

(generate-all-array-tests axpy :float-complex
245
 (let ((v1 (array-default 8))
Liam Healy's avatar
Liam Healy committed
246 247
	(v2 (array-default 8))
	(scalar (scalar-default)))
248
   (grid:copy-to (axpy scalar v1 v2))))
Liam Healy's avatar
Liam Healy committed
249

Liam Healy's avatar
Liam Healy committed
250
(generate-all-array-tests scale #+fsbv :float-complex #-fsbv :float
251
 (let ((v1 (array-default 8))
Liam Healy's avatar
Liam Healy committed
252
	(scalar (scalar-default)))
253
   (grid:copy-to (scale scalar v1))))
Liam Healy's avatar
Liam Healy committed
254 255

(generate-all-array-tests scale :complex
256 257
 (let ((v1 (array-default 8))
       (scalar (scalar-default t)))
258
   (grid:copy-to (scale scalar v1))))
Liam Healy's avatar
Liam Healy committed
259 260

(generate-all-array-tests givens :float
261 262 263 264 265
 (let ((v1 (array-default 8))
       (v2 (array-default 8))
       (angles (array-default 8))
       (sines (array-default 8 t))
       (cosines (array-default 8 t)))
Liam Healy's avatar
Liam Healy committed
266
   (loop for i below 8 do
Liam Healy's avatar
Liam Healy committed
267 268
	(setf (grid:gref sines i) (sin (grid:gref angles i)))
	(setf (grid:gref cosines i) (cos (grid:gref angles i))))
Liam Healy's avatar
Liam Healy committed
269
   (givens-rotation v1 v2 cosines sines)
270
   (list (grid:copy-to v1) (grid:copy-to v2))))