-
Notifications
You must be signed in to change notification settings - Fork 1
/
cblaslib.l
152 lines (136 loc) · 5.5 KB
/
cblaslib.l
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
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
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
;;; This file was automatically generated by SWIG (http://www.swig.org).
;;; Version 2.0.9
;;;
;;; Do not make changes to this file unless you know what you are doing--modify
;;; the SWIG interface file instead.
;; (setq *cblaslib* (load-foreign (format nil "/usr/lib/libcblas.so")))
(setq *cblaslib* (load-foreign (format nil "/usr/lib/openblas-base/libblas.so")))
(defun make-float-vector (n &optional (v 0.0))
(make-sequence float-vector n :initial-element v))
(defun matrix-reshape (m row &optional (col (/ (array-total-size m) row)))
(setq (m . dim0) row (m . dim1) col)
m)
(defun float-vector-arange (n &aux (v (instantiate float-vector n)))
(dotimes (i n) (setf (elt v i) i)) v)
(defun matrix-arange (n &aux (m (make-matrix 1 n)))
(dotimes (i n) (setf (elt (m . entity) i) i))
m)
(defun array-arange (n &aux (m (make-array (list 1 n))))
(dotimes (i n) (setf (elt (m . entity) i) i))
m)
(defun array-shape (a &aux s)
(dotimes (i (array-rank a)) (push (slot a array (+ i 5)) s))
(nreverse s))
(defun array-reshape (a shape)
(setq (a . rank) (length shape))
(dotimes (i (length shape)) (setslot a array (+ 5 i) (pop shape)))
a)
(defforeign cblas_dnrm2 *cblaslib* "cblas_dnrm2" () :float)
;;double cblas_dnrm2(const int N, const double *X, const int incX);
;; (sqrt (x . x))
(defun cblas-dnrm2 (x &key (incx 1)) (cblas_dnrm2 (length x) x incx))
(defforeign cblas_ddot *cblaslib* "cblas_ddot" () :float)
;;double cblas_ddot(const int N, const double *X, const int incX,
;; const double *Y, const int incY);
(defun cblas-ddot (x y &key (incx 1) (incy 1))
(cblas_ddot (length x) x incx y incy))
(defforeign cblas_dcopy *cblaslib* "cblas_dcopy" () :integer)
;;void cblas_dcopy(const int N, const double *X, const int incX,
;; double *Y, const int incY);
(defun cblas-dcopy (x y &key (incx 1) (incy 1))
(cblas_dcopy (length y) y incy x incx)
x)
(defforeign cblas_dger *cblaslib* "cblas_dger" () :integer)
;;void cblas_dger(CBLAS_LAYOUT layout, const int M, const int N,
;; const double alpha, const double *X, const int incX,
;; const double *Y, const int incY, double *A, const int lda);
;; A = alpha* x * y^t + A
(defun cblas-dger (xv yv &key (incx 1) (incy 1) (alpha 1.0))
(let* ((major 101) ;; 101: rowmajor, 102: colmajor
(m (length xv))
(n (length yv))
(mat (make-matrix m n)))
(cblas_dger major m n alpha xv incx yv incy (mat . entity) (max m n))
mat))
(defforeign cblas_daxpy *cblaslib* "cblas_daxpy" () :integer)
;;void cblas_daxpy(const int N, const double alpha, const double *X,
;; const int incX, double *Y, const int incY);
;; Y = alpha * X + Y
(defun cblas-daxpy (x y &key (alpha 1.0) (incx 1) (incy 1))
(cblas_daxpy (length x) alpha x incx y incy)
y)
(defun cblas-mat- (x &optional y)
(if y
(cblas-daxpy (y . entity) (x . entity) :alpha -1.0)
(cblas-daxpy (x . entity) (x . entity) :alpha -2.0)) ;; modified
x)
(defun cblas-mat+ (x y)
(cblas-daxpy (y . entity) (x . entity))
x)
(defforeign cblas_dgemv *cblaslib* "cblas_dgemv" () :integer)
;;void cblas_dgemv(CBLAS_LAYOUT layout,
;; CBLAS_TRANSPOSE TransA, const int M, const int N,
;; const double alpha, const double *A, const int lda,
;; const double *X, const int incX, const double beta,
;; double *Y, const int incY);
(defun cblas-dgemv (a x y &key (alpha 1.0) (beta 0.0) (incx 1) (incy 1))
(let* ((major 101) ;; 101: rowmajor, 102: colmajor
(trans 111) ;; 111: notrans, 112: trans, 113: conjtrans
(m (array-dimension a 0))
(n (array-dimension a 1))
)
(cblas_dgemv major trans
m n alpha (a . entity) (max m n)
x incx beta y incy)
y))
(defforeign cblas_dgemm *cblaslib* "cblas_dgemm" () :integer)
;; ( CBLAS_LAYOUT CBLAS_TRANSPOSE CBLAS_TRANSPOSE :integer :integer :integer :float :integer :integer :integer :integer :float :integer :integer)
(defun cblas-dgemm(a b c &optional (alpha 1.0) (beta 1.0))
(let ((major 101) ;; 101: rowmajor, 102: colmajor
(trans 111) ;; 111: notrans, 112: trans, 113: conjtrans
(m (array-dimension a 0))
(n (array-dimension b 1))
(k (array-dimension a 1))
)
(cblas_dgemm major trans trans
m n k alpha
(a . entity) k
(b . entity) n
beta
(c . entity) n)
c)
)
(defun test1 nil
(setq m1 (make-matrix 2 4 '((1.0 2.0 3.0 4.0) (1.0 2.0 3.0 4.0))))
(setq m2 (make-matrix 4 3 '((1.0 2.0 3.0) (1.0 2.0 3.0) (1.0 2.0 3.0) (1.0 2.0 3.0))))
(setq m3 (make-matrix 2 3))
(setq v4 (float-vector 1 2 3 4))
(setq v2 (float-vector 1 1))
(cblas-dgemm m1 m2 m3 1.0 0.0)
(cblas-dgemv m1 v4 v2)
(cblas-dnrm2 v4)
)
(defun test2 nil
(setq m1 (unit-matrix 5000))
(setq m2 (unit-matrix 5000))
(setq m3 (unit-matrix 5000))
(timing 3 (cblas-mat- m1 m2)) ;; 70.000 milisec
(timing 3 (m- m1 m2 m3)) ;; 53.333 milisec
(timing 3 (cblas-mat+ m1 m2)) ;; 73.333 milisec
(timing 3 (m+ m1 m2 m3)) ;; 53.333 milisec
)
(defun test3 nil
(setq m1 (make-matrix 1000 10000))
(setq m2 (make-matrix 10000 500))
(setq m3 (make-matrix 1000 500))
(setq v4 (make-float-vector 10000 1.0))
(setq v1 (make-float-vector 1000 1.0))
(setq v2 (make-float-vector 1000 1.0))
(setq v10 (make-float-vector 10000 1.0))
(timing 10 (cblas-dgemv m1 v4 v1)) ;; 5.000 milisec
(timing 10 (transform m1 v4 v2)) ;; 15.000 milisec
(timing 1000 (cblas-dnrm2 v10)) ;; 0.010 milisec
(timing 1000 (cblas-ddot v10 v10)) ;; 0.010 milisec
(timing 10 (cblas-dgemm m1 m2 m3)) ;; 1990.000 milisec
(timing 1 (m* m1 m2 m3)) ;; 25305.000 milisec
)