Menu
▾
▴
matlisp-commit — Commit log
You can subscribe to this list here.
| 2012 |
Jan
|
Feb
|
Mar
(34) |
Apr
(4) |
May
(2) |
Jun
(11) |
Jul
(22) |
Aug
(9) |
Sep
|
Oct
|
Nov
|
Dec
(4) |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 2013 |
Jan
(15) |
Feb
(17) |
Mar
(3) |
Apr
|
May
|
Jun
(3) |
Jul
(1) |
Aug
(5) |
Sep
(5) |
Oct
(2) |
Nov
(1) |
Dec
(1) |
| 2014 |
Jan
|
Feb
(1) |
Mar
(1) |
Apr
|
May
(1) |
Jun
|
Jul
|
Aug
|
Sep
|
Oct
|
Nov
|
Dec
|
| 2016 |
Jan
|
Feb
|
Mar
|
Apr
|
May
|
Jun
|
Jul
|
Aug
(1) |
Sep
|
Oct
|
Nov
|
Dec
|
|
From: Akshay S. <ak...@us...> - 2012-03-17 06:05:54
|
This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "matlisp".
The branch, matlisp-cffi has been updated
via 54341c25f149263190e4ffad1c516d93a79ad3ed (commit)
via 146847f922138620cb4b1ad064cc2ad8f80bb304 (commit)
from 517949e31b41b303dab91670e80b207bc45d3256 (commit)
Those revisions listed above that are new to this repository have
not appeared on any other notification email; so we list those
revisions in full, below.
- Log -----------------------------------------------------------------
commit 54341c25f149263190e4ffad1c516d93a79ad3ed
Author: Akshay Srinivasan <aks...@gm...>
Date: Sat Mar 17 11:32:32 2012 +0530
-> Gemv works!
diff --git a/packages.lisp b/packages.lisp
index 566dbd6..74385fa 100644
--- a/packages.lisp
+++ b/packages.lisp
@@ -154,8 +154,10 @@
;;; Define the packages and symbols for Matlisp.
(defpackage :utilities
- (:use :common-lisp)
- (:export #:lcase
+ (:use :common-lisp)
+ (:export #:ensure-list
+ #:zip
+ #:zip-eq
#:cut-cons-chain!
#:when-let
#:if-ret
diff --git a/src/gemv.lisp b/src/gemv.lisp
index 0e49fac..c58f746 100644
--- a/src/gemv.lisp
+++ b/src/gemv.lisp
@@ -39,8 +39,6 @@
(* beta (matrix-ref-2d y i 0)))))))))
y))
-
-
;;
(defgeneric gemv! (alpha A x beta y &optional job)
(:documentation
@@ -100,6 +98,7 @@
(defmethod gemv! ((alpha cl:real) (A real-matrix) (x real-matrix)
(beta cl:real) (y real-matrix) &optional (job :n))
+ ;; y <- \beta . y + \alpha . A o x
(real-double-gemv!-typed (coerce alpha 'double-float) A x
(coerce beta 'double-float) y job))
@@ -110,57 +109,126 @@
(defmethod gemv! ((alpha number) (A complex-matrix) (x complex-matrix)
(beta number) (y complex-matrix) &optional (job :n))
+ ;; y <- \beta . y + \alpha . A o x
(complex-double-gemv!-typed (complex-coerce alpha) A x
(complex-coerce beta) y job))
;
-(defmethod gemv! ((alpha cl:real) (A real-matrix) (x real-matrix)
- (beta cl:real) (y complex-matrix) &optional (job :n))
- (real-double-gemv!-typed (coerce alpha 'double-float) A x
- (coerce beta 'double-float) (realpart! y) job))
-
-(defmethod gemv! ((alpha cl:real) (A real-matrix) (x real-matrix)
+(defmethod gemv! ((alpha number) (A real-matrix) (x real-matrix)
(beta complex) (y complex-matrix) &optional (job :n))
+ ;; y <- \beta * y
(scal! (complex-coerce beta) y)
- (real-double-gemv!-typed (coerce alpha 'double-float) A x
- 1d0 (realpart! y) job))
+ ;; y <- y + \alpha * A o x
+ (gemv! alpha A x 1d0 y job))
-;
+(defmethod gemv! ((alpha cl:real) (A real-matrix) (x real-matrix)
+ (beta cl:real) (y complex-matrix) &optional (job :n))
+ (let ((r-be (coerce beta 'double-float))
+ (r-al (coerce alpha 'double-float))
+ (r-y (realpart! y)))
+ (declare (type double-float r-be r-al)
+ (type real-matrix r-y))
+ ;; y <- \beta * y
+ (scal! r-be y)
+ ;; (realpart! y) <- (realpart! y) + \alpha * A o x
+ (real-double-gemv!-typed r-al A x 1d0 r-y job))
+ y)
(defmethod gemv! ((alpha complex) (A real-matrix) (x real-matrix)
(beta cl:real) (y complex-matrix) &optional (job :n))
- (real-double-gemv!-typed (coerce (realpart alpha) 'double-float) A x
- (coerce beta 'double-float) (realpart! y) job)
- (real-double-gemv!-typed (coerce (imagpart alpha) 'double-float) A x
- (coerce beta 'double-float) (imagpart! y) job))
+ (let ((r-al (coerce (realpart alpha) 'double-float))
+ (i-al (coerce (imagpart alpha) 'double-float))
+ (r-be (coerce beta 'double-float))
+ (r-y (realpart! y))
+ (i-y (imagpart! y)))
+ (declare (type double-float r-al i-al r-be)
+ (type real-matrix r-y i-y))
+ ;; (realpart! y) <- \beta * (realpart! y) + (realpart \alpha) . A o x
+ (real-double-gemv!-typed r-al A x r-be r-y job)
+ ;; (imagpart! y) <- \beta * (imagpart! y) + (imagpart \alpha) . A o x
+ (real-double-gemv!-typed i-al A x r-be i-y job))
+ y)
-(defmethod gemv! ((alpha complex) (A real-matrix) (x real-matrix)
- (beta complex) (y complex-matrix) &optional (job :n))
- (scal! (complex-coerce beta) y)
- (real-double-gemv!-typed (coerce (realpart alpha) 'double-float) A x
- 1d0 (realpart! y) job)
- (real-double-gemv!-typed (coerce (imagpart alpha) 'double-float) A x
- 1d0 (imagpart! y) job))
;
-
(defmethod gemv! ((alpha number) (A real-matrix) (x complex-matrix)
- (beta number) (y complex-matrix) &optional (job :n))
- (gemv! alpha A (realpart! x)
- beta y job)
- (gemv! (* #c(0d0 1d0) alpha) A (imagpart! x)
- beta y job))
-;
-(defmethod gemv! ((alpha number) (A complex-matrix) (x real-matrix)
- (beta number) (y complex-matrix) &optional (job :n))
- (gemv! alpha (realpart! A) x
- beta y job)
- (gemv! (* #c(0d0 1d0) alpha) (imagpart! A) x
- beta y job))
+ (beta complex) (y complex-matrix) &optional (job :n))
+ ;; y <- \beta y
+ (scal! beta y)
+ ;; y <- y + \alpha . A o x
+ (gemv! alpha A x 1d0 y job))
+
+(defmethod gemv! ((alpha cl:real) (A real-matrix) (x complex-matrix)
+ (beta cl:real) (y complex-matrix) &optional (job :n))
+ (let ((r-x (realpart! x))
+ (i-x (imagpart! x))
+ (r-y (realpart! y))
+ (i-y (imagpart! y))
+ (r-al (coerce (realpart alpha) 'double-float))
+ (r-be (coerce beta 'double-float)))
+ (declare (type double-float r-al r-be)
+ (type real-matrix r-x i-x r-y i-y))
+ ;; (realpart! y) <- \beta * (realpart! y) + \alpha . A o (realpart! x)
+ (real-double-gemv!-typed r-al A r-x r-be r-y job)
+ ;; (imagpart! y) <- \beta * (imagpart! y) + \alpha . A o (realpart! x)
+ (real-double-gemv!-typed r-al A i-x r-be i-y job))
+ y)
+
+(defmethod gemv! ((alpha complex) (A real-matrix) (x complex-matrix)
+ (beta cl:real) (y complex-matrix) &optional (job :n))
+ (let ((r-x (realpart! x))
+ (i-x (imagpart! x))
+ (r-y (realpart! y))
+ (i-y (imagpart! y))
+ (r-al (coerce (realpart alpha) 'double-float))
+ (i-al (coerce (imagpart alpha) 'double-float))
+ (r-be (coerce beta 'double-float)))
+ (declare (type double-float r-al r-be i-al)
+ (type real-matrix r-x i-x r-y i-y))
+ (real-double-gemv!-typed r-al A r-x r-be r-y job)
+ (real-double-gemv!-typed (- i-al) A i-x 1d0 r-y job)
+ ;;
+ (real-double-gemv!-typed i-al A r-x r-be i-y job)
+ (real-double-gemv!-typed r-al A i-x 1d0 i-y job))
+ y)
-
;
-(defun gemv! ((alpha number) (A complex-matrix) (x real-matrix)
- (beta number) (y complex-matrix))
- )
+(defmethod gemv! ((alpha number) (A complex-matrix) (x real-matrix)
+ (beta complex) (y complex-matrix) &optional (job :n))
+ ;; y <- \beta y
+ (scal! beta y)
+ ;; y <- y + \alpha . A o x
+ (gemv! alpha A x 1d0 y job))
-;;
\ No newline at end of file
+(defmethod gemv! ((alpha cl:real) (A complex-matrix) (x real-matrix)
+ (beta cl:real) (y complex-matrix) &optional (job :n))
+ (let ((r-A (realpart! A))
+ (i-A (imagpart! A))
+ (r-y (realpart! y))
+ (i-y (imagpart! y))
+ (r-al (coerce (realpart alpha) 'double-float))
+ (r-be (coerce beta 'double-float)))
+ (declare (type double-float r-al r-be)
+ (type real-matrix r-A i-A r-y i-y))
+ ;; (realpart! y) <- \beta * (realpart! y) + \alpha . A o (realpart! x)
+ (real-double-gemv!-typed r-al r-A x r-be r-y job)
+ ;; (imagpart! y) <- \beta * (imagpart! y) + \alpha . A o (realpart! x)
+ (real-double-gemv!-typed r-al i-A x r-be i-y job))
+ y)
+
+(defmethod gemv! ((alpha complex) (A complex-matrix) (x real-matrix)
+ (beta cl:real) (y complex-matrix) &optional (job :n))
+ (let ((r-A (realpart! A))
+ (i-A (imagpart! A))
+ (r-y (realpart! y))
+ (i-y (imagpart! y))
+ (r-al (coerce (realpart alpha) 'double-float))
+ (i-al (coerce (imagpart alpha) 'double-float))
+ (r-be (coerce beta 'double-float)))
+ (declare (type double-float r-al r-be i-al)
+ (type real-matrix r-A i-A r-y i-y))
+ (real-double-gemv!-typed r-al r-A x r-be r-y job)
+ (real-double-gemv!-typed (- i-al) i-A x 1d0 r-y job)
+ ;;
+ (real-double-gemv!-typed i-al r-A x r-be i-y job)
+ (real-double-gemv!-typed r-al i-A x 1d0 i-y job))
+ y)
\ No newline at end of file
diff --git a/src/utilities.lisp b/src/utilities.lisp
index 76caf8c..6cbf2bc 100644
--- a/src/utilities.lisp
+++ b/src/utilities.lisp
@@ -145,15 +145,23 @@ else run else-body"
(cut-cons-chain-tin lst test lst)))
;;
-(defmacro lcase (keyform &rest body)
- (let ((key-eval (gensym)))
- (labels ((app-equal (lst)
- (if (null lst)
- nil
- `(((equal ,key-eval ,(caar lst)) ,@(cdar lst))
- ,@(app-equal (cdr lst))))))
- `(let ((,key-eval ,keyform))
- (cond
- ,@(app-equal body))))))
+(defun zip (&rest args)
+ (apply #'map 'list #'list args))
+
+;;
+(defmacro mcase (keyform &rest body)
+ (labels ((app-equal (lst)
+ (if (null lst)
+ nil
+ `(((and ,@(mapcar (lambda (pair) (cons 'eq pair))
+ (zip keyform (caar lst))))
+ ,@(cdar lst))
+ ,@(app-equal (cdr lst))))))
+ `(cond
+ ,@(app-equal body))))
+
+(defmacro zip-eq (a b)
+ `(and ,@(mapcar (lambda (pair) (cons 'eq pair))
+ (zip (ensure-list a) (ensure-list b)))))
;;
commit 146847f922138620cb4b1ad064cc2ad8f80bb304
Author: Akshay Srinivasan <aks...@gm...>
Date: Sat Mar 17 08:53:49 2012 +0530
Gemv core function works reasonably well.
diff --git a/src/axpy.lisp b/src/axpy.lisp
index d6cee37..f001fd8 100644
--- a/src/axpy.lisp
+++ b/src/axpy.lisp
@@ -85,20 +85,18 @@
((hd-b st-b) (slot-values mat-b '(head store)) :type (fixnum (,store-type *))))
(if (and cp-a cp-b)
(,blas-func sz alpha st-a inc-a st-b inc-b :head-x hd-a :head-y hd-b)
- (symbol-macrolet
- ((common-code
- (mlet* (((nr-a nc-a rs-a cs-a) (slot-values mat-a '(number-of-rows number-of-cols row-stride col-stride))
- :type (fixnum fixnum fixnum fixnum))
- ((rs-b cs-b) (slot-values mat-b '(row-stride col-stride))
- :type (fixnum fixnum)))
- (loop for i from 0 below nr-a
- do (,blas-func nc-a alpha st-a cs-a st-b cs-b :head-x (+ hd-a (* i rs-a)) :head-y (+ hd-b (* i rs-b)))))))
- ;;Choose the smaller of the loops
- (if (> (nrows mat-a) (ncols mat-a))
- (with-transpose! (mat-a mat-b)
- common-code)
- common-code)))
- mat-b)))
+ (mlet* (((nr-a nc-a rs-a cs-a) (slot-values mat-a '(number-of-rows number-of-cols row-stride col-stride))
+ :type (fixnum fixnum fixnum fixnum))
+ ((rs-b cs-b) (slot-values mat-b '(row-stride col-stride))
+ :type (fixnum fixnum)))
+ ;;Choose the smaller of the loops
+ (when (> nr-a nc-a)
+ (rotatef nr-a nc-a)
+ (rotatef rs-a cs-a)
+ (rotatef rs-b cs-b))
+ (loop for i from 0 below nr-a
+ do (,blas-func nc-a alpha st-a cs-a st-b cs-b :head-x (+ hd-a (* i rs-a)) :head-y (+ hd-b (* i rs-b)))))))
+ mat-b))
;;
(defgeneric axpy! (alpha x y)
diff --git a/src/blas.lisp b/src/blas.lisp
index 901b2fc..e7f76d9 100644
--- a/src/blas.lisp
+++ b/src/blas.lisp
@@ -498,19 +498,16 @@
(def-fortran-routine mzdotu :void
(result (* :complex-double-float) :output)
(n :integer :input)
- (zx (* :complex-double-float) :input)
+ (zx (* :complex-double-float :inc head-x) :input)
(incx :integer :input)
- (zy (* :complex-double-float) :input)
+ (zy (* :complex-double-float :inc head-y) :input)
(incy :integer :input)
)
(let ((result (make-array 2 :element-type 'double-float)))
- (defun zdotu (n zx incx zy incy)
- (with-vector-data-addresses ((addr-result result)
- (addr-zx zx)
- (addr-zy zy))
- (mzdotu addr-result n addr-zx incx addr-zy incy)
- (complex (aref result 0) (aref result 1)))))
+ (defun zdotu (n zx incx zy incy &key (head-x 0) (head-y 0))
+ (mzdotu result n zx incx zy incy :head-x head-x :head-y head-y)
+ (complex (aref result 0) (aref result 1))))
#-(and)
(def-fortran-routine zdotc :complex-double-float
@@ -559,19 +556,16 @@
(def-fortran-routine mzdotc :void
(result (* :complex-double-float) :output)
(n :integer :input)
- (zx (* :complex-double-float) :input)
+ (zx (* :complex-double-float :inc head-x) :input)
(incx :integer :input)
- (zy (* :complex-double-float) :input)
+ (zy (* :complex-double-float :inc head-y) :input)
(incy :integer :input)
)
(let ((result (make-array 2 :element-type 'double-float)))
- (defun zdotc (n zx incx zy incy)
- (with-vector-data-addresses ((addr-result result)
- (addr-zx zx)
- (addr-zy zy))
- (mzdotc addr-result n addr-zx incx addr-zy incy)
- (complex (aref result 0) (aref result 1)))))
+ (defun zdotc (n zx incx zy incy &key (head-x 0) (head-y 0))
+ (mzdotc result n zx incx zy incy :head-x head-x :head-y head-y)
+ (complex (aref result 0) (aref result 1))))
(def-fortran-routine idamax :integer
@@ -2143,12 +2137,12 @@
(m :integer )
(n :integer )
(alpha :complex-double-float )
- (a (* :complex-double-float) )
+ (a (* :complex-double-float :inc head-a) )
(lda :integer )
- (x (* :complex-double-float) )
+ (x (* :complex-double-float :inc head-x) )
(incx :integer )
(beta :complex-double-float )
- (y (* :complex-double-float) :output)
+ (y (* :complex-double-float :inc head-y) :output)
(incy :integer )
)
diff --git a/src/complex-matrix.lisp b/src/complex-matrix.lisp
index 0e02c9b..704f823 100644
--- a/src/complex-matrix.lisp
+++ b/src/complex-matrix.lisp
@@ -97,8 +97,8 @@
(store (allocate-complex-store size)))
(multiple-value-bind (row-stride col-stride)
(ecase order
- (:row-major (values n 1))
- (:col-major (values 1 m)))
+ (:row-major (values m 1))
+ (:col-major (values 1 n)))
(let ((matrix
(make-instance 'complex-matrix
:nrows n :ncols m
@@ -306,3 +306,17 @@
:nrows (nrows mat) :ncols (ncols mat)
:row-stride (* 2 (row-stride mat)) :col-stride (* 2 (col-stride mat))
:head (+ 1 (* 2 (head mat)))))))
+
+
+(defun conjugate! (mat)
+ (cond
+ ((typep mat 'real-matrix) nil)
+ ((typep mat 'complex-matrix) (progn
+ (transpose! (scal! -1d0 (imagpart! mat)))
+ mat))))
+
+(defmacro with-conjugate! (matlst &rest body)
+ `(progn
+ ,@(mapcar #'(lambda (mat) `(conjugate! ,mat)) matlst)
+ ,@body
+ ,@(mapcar #'(lambda (mat) `(conjugate! ,mat)) matlst)))
\ No newline at end of file
diff --git a/src/copy.lisp b/src/copy.lisp
index 03c6e8e..8ce687d 100644
--- a/src/copy.lisp
+++ b/src/copy.lisp
@@ -82,7 +82,6 @@
(defmacro generate-typed-copy!-func (func store-type matrix-type blas-func)
;;Be very careful when using functions generated by this macro.
;;Indexes can be tricky and this has no safety net
- ;;(you don't see a matrix-ref do you ?)
;;Use only after checking the arguments for compatibility.
`(defun ,func (mat-a mat-b)
(declare (type ,matrix-type mat-a mat-b)
diff --git a/src/gemv.lisp b/src/gemv.lisp
index cffa7d6..0e49fac 100644
--- a/src/gemv.lisp
+++ b/src/gemv.lisp
@@ -1,32 +1,166 @@
(in-package :matlisp)
-(defun real-double-gemv!-typed (alpha A x beta y job)
- ;;Be very careful when you use this function!
+;;There's no support for ":c", because there is no
+;;equivalent of ":n" with complex conjugation.
+(defmacro generate-typed-gemv!-func (func element-type store-type matrix-type blas-gemv-func blas-axpy-func blas-dot-func)
+ ;;Be very careful when you use functions generated by this macro!
;;Indexes can be tricky and this has no safety net
- ;;(you don't see a matrix-ref do you ?)
;;Use only after checking the arguments for compatibility.
- (declare (type double-float alpha beta)
- (type real-matrix A x y)
- (type symbol job))
- (mlet* ((fort-op (ecase job (:n "N") (:t "T")) :type ((string 1)))
- ((st-a hd-a nr-a nc-a rs-a cs-a) (slot-values A '(store head number-of-rows number-of-cols row-stride col-stride))
- :type ((real-matrix-store-type *) fixnum fixnum fixnum fixnum fixnum))
- ((st-x hd-x rs-x) (slot-values x '(store head row-stride))
- :type ((real-matrix-store-type *) fixnum fixnum))
- ((st-y hd-y rs-y) (slot-values y '(store head row-stride))
- :type ((real-matrix-store-type *) fixnum fixnum))
- ((sym lda tf-op) (blas-matrix-compatible-p A fort-op) :type (symbol fixnum (string 1))))
- (if (not (string= tf-op "?"))
- (progn
- (when (eq sym :row-major)
- (rotatef nr-a nc-a)
- (rotatef rs-a cs-a))
- (blas:dgemv tf-op nr-a nc-a alpha st-a lda st-x rs-x beta st-y rs-y :head-a hd-a :head-x hd-x :head-y hd-y))
- (progn
- (when (string= fort-op "T")
- (rotatef nr-a nc-a)
- (rotatef rs-a cs-a))
- (dotimes (i nr-a)
- (setf (matrix-ref-2d y i 0) (+ (* alpha (blas:ddot nc-a st-a cs-a st-x rs-x :head-x (+ hd-a (* i rs-a)) :head-y hd-x))
- (* beta (matrix-ref-2d y i 0))))))))
- y)
\ No newline at end of file
+ `(defun ,func (alpha A x beta y job)
+ (declare (type ,element-type alpha beta)
+ (type ,matrix-type A x y)
+ (type symbol job))
+ (mlet* ((fort-op (ecase job (:n "N") (:t "T")) :type ((string 1)))
+ ((st-a hd-a nr-a nc-a rs-a cs-a) (slot-values A '(store head number-of-rows number-of-cols row-stride col-stride))
+ :type ((,store-type *) fixnum fixnum fixnum fixnum fixnum))
+ ((st-x hd-x rs-x) (slot-values x '(store head row-stride))
+ :type ((,store-type *) fixnum fixnum))
+ ((st-y hd-y rs-y) (slot-values y '(store head row-stride))
+ :type ((,store-type *) fixnum fixnum))
+ ((sym lda tf-op) (blas-matrix-compatible-p A fort-op) :type (symbol fixnum (string 1))))
+ (if (not (string= tf-op "?"))
+ (progn
+ (when (eq sym :row-major)
+ (rotatef nr-a nc-a)
+ (rotatef rs-a cs-a))
+ (,blas-gemv-func tf-op nr-a nc-a alpha st-a lda st-x rs-x beta st-y rs-y :head-a hd-a :head-x hd-x :head-y hd-y))
+ (progn
+ (when (string= fort-op "T")
+ (rotatef nr-a nc-a)
+ (rotatef rs-a cs-a))
+ ;;Use the smaller of the loops.
+ (if (> nr-a nc-a)
+ (progn
+ (scal! beta y)
+ (dotimes (i nc-a)
+ (,blas-axpy-func nr-a (* alpha (matrix-ref-2d x i 0)) st-a rs-a st-y rs-y :head-x (+ hd-a (* i cs-a)) :head-y hd-y)))
+ (dotimes (i nr-a)
+ (setf (matrix-ref-2d y i 0) (+ (* alpha (,blas-dot-func nc-a st-a cs-a st-x rs-x :head-x (+ hd-a (* i rs-a)) :head-y hd-x))
+ (* beta (matrix-ref-2d y i 0)))))))))
+ y))
+
+
+
+;;
+(defgeneric gemv! (alpha A x beta y &optional job)
+ (:documentation
+"
+ Syntax
+ ======
+ (GEMV! alpha A x beta y [job])
+
+ Purpose
+ =======
+ Performs the GEneral Matrix Vector operation given by
+ -- - -
+
+ Y <- alpha * op(A) * x + beta * y
+
+ and returns y.
+
+ alpha,beta are scalars,
+ A is a matrix, and x,y are vectors.
+
+ op(A) means either A or A'.
+
+ JOB Operation
+ ---------------------------------------------------
+ :N (default) alpha * A * x + beta * y
+ :T alpha * A'* x + beta * y
+
+ Note
+ ====
+ Take caution when using GEMM! as follows:
+
+ (GEMV! alpha a x beta x)
+
+ The results may be unpredictable depending
+ on the underlying DGEMM, ZGEMM routines
+ from BLAS, ATLAS or LIBCRUFT.
+"))
+
+(defmethod gemv! :before ((alpha number) (A standard-matrix) (x standard-matrix)
+ (beta number) (y standard-matrix)
+ &optional (job :n))
+ (mlet* (((nr-a nc-a) (slot-values A '(number-of-rows number-of-cols)) :type (fixnum fixnum))
+ ((nr-x nc-x) (slot-values x '(number-of-rows number-of-cols)) :type (fixnum fixnum))
+ ((nr-y nc-y) (slot-values y '(number-of-rows number-of-cols)) :type (fixnum fixnum)))
+ (unless (member job '(:n :t))
+ (error "Argument JOB to GEMV! is not recognized"))
+ (when (eq job :t)
+ (rotatef nr-a nc-a))
+ (unless (and (= nc-x 1) (= nc-y 1)
+ (= nc-a nr-x) (= nr-a nr-y))
+ (error "Dimensions of A,x,y given to GEMV! do not match"))))
+
+;;
+(generate-typed-gemv!-func real-double-gemv!-typed
+ double-float real-matrix-store-type real-matrix
+ blas:dgemv blas:daxpy blas:ddot)
+
+(defmethod gemv! ((alpha cl:real) (A real-matrix) (x real-matrix)
+ (beta cl:real) (y real-matrix) &optional (job :n))
+ (real-double-gemv!-typed (coerce alpha 'double-float) A x
+ (coerce beta 'double-float) y job))
+
+;;
+(generate-typed-gemv!-func complex-double-gemv!-typed
+ complex-double-float complex-matrix-store-type complex-matrix
+ blas:zgemv blas:zaxpy blas:zdotu)
+
+(defmethod gemv! ((alpha number) (A complex-matrix) (x complex-matrix)
+ (beta number) (y complex-matrix) &optional (job :n))
+ (complex-double-gemv!-typed (complex-coerce alpha) A x
+ (complex-coerce beta) y job))
+
+;
+(defmethod gemv! ((alpha cl:real) (A real-matrix) (x real-matrix)
+ (beta cl:real) (y complex-matrix) &optional (job :n))
+ (real-double-gemv!-typed (coerce alpha 'double-float) A x
+ (coerce beta 'double-float) (realpart! y) job))
+
+(defmethod gemv! ((alpha cl:real) (A real-matrix) (x real-matrix)
+ (beta complex) (y complex-matrix) &optional (job :n))
+ (scal! (complex-coerce beta) y)
+ (real-double-gemv!-typed (coerce alpha 'double-float) A x
+ 1d0 (realpart! y) job))
+
+;
+
+(defmethod gemv! ((alpha complex) (A real-matrix) (x real-matrix)
+ (beta cl:real) (y complex-matrix) &optional (job :n))
+ (real-double-gemv!-typed (coerce (realpart alpha) 'double-float) A x
+ (coerce beta 'double-float) (realpart! y) job)
+ (real-double-gemv!-typed (coerce (imagpart alpha) 'double-float) A x
+ (coerce beta 'double-float) (imagpart! y) job))
+
+(defmethod gemv! ((alpha complex) (A real-matrix) (x real-matrix)
+ (beta complex) (y complex-matrix) &optional (job :n))
+ (scal! (complex-coerce beta) y)
+ (real-double-gemv!-typed (coerce (realpart alpha) 'double-float) A x
+ 1d0 (realpart! y) job)
+ (real-double-gemv!-typed (coerce (imagpart alpha) 'double-float) A x
+ 1d0 (imagpart! y) job))
+;
+
+(defmethod gemv! ((alpha number) (A real-matrix) (x complex-matrix)
+ (beta number) (y complex-matrix) &optional (job :n))
+ (gemv! alpha A (realpart! x)
+ beta y job)
+ (gemv! (* #c(0d0 1d0) alpha) A (imagpart! x)
+ beta y job))
+;
+(defmethod gemv! ((alpha number) (A complex-matrix) (x real-matrix)
+ (beta number) (y complex-matrix) &optional (job :n))
+ (gemv! alpha (realpart! A) x
+ beta y job)
+ (gemv! (* #c(0d0 1d0) alpha) (imagpart! A) x
+ beta y job))
+
+
+;
+(defun gemv! ((alpha number) (A complex-matrix) (x real-matrix)
+ (beta number) (y complex-matrix))
+ )
+
+;;
\ No newline at end of file
diff --git a/src/standard-matrix.lisp b/src/standard-matrix.lisp
index 0f1ec71..089ef9b 100644
--- a/src/standard-matrix.lisp
+++ b/src/standard-matrix.lisp
@@ -228,25 +228,6 @@ Matrix only has ~A elements." idx (number-of-elements matrix))))
(list (nrows matrix) (ncols matrix)))
;;
-(defgeneric transpose! (matrix)
- (:documentation
- "
- Syntax
- ======
- (transpose! matrix)
-
- Purpose
- =======
- Exchange row and column strides so that effectively
- the matrix is transposed in place (without much effort).
-"))
-
-(defmethod transpose! ((matrix standard-matrix))
- (rotatef (nrows matrix) (ncols matrix))
- (rotatef (row-stride matrix) (col-stride matrix))
- matrix)
-
-;;
(defgeneric fill-matrix (matrix fill-element)
(:documentation
"
@@ -277,7 +258,58 @@ matrix and a number"))
(format stream "~%")))
;;
-(defun make-sub-matrix (matrix i j nrows ncols)
+(defun transpose-i! (matrix)
+"
+ Syntax
+ ======
+ (transpose-i! matrix)
+
+ Purpose
+ =======
+ Exchange row and column strides so that effectively
+ the matrix is transposed in place (without much effort).
+"
+ (cond
+ ((typep matrix 'standard-matrix)
+ (progn
+ (rotatef (nrows matrix) (ncols matrix))
+ (rotatef (row-stride matrix) (col-stride matrix))
+ matrix))
+ ((typep matrix 'number) matrix)
+ (t (error "Don't know how to take the transpose of ~A." matrix))))
+
+;;
+(defun transpose! (matrix)
+"
+ Syntax
+ ======
+ (transpose! matrix)
+
+ Purpose
+ =======
+ Create a new matrix object which represents the transpose of the
+ the given matrix, but shares the store with matrix.
+"
+ (cond
+ ((typep matrix 'standard-matrix)
+ (mlet* (((hd nr nc rs cs) (slot-values matrix '(head number-of-rows number-of-cols row-stride col-stride)) :type (fixnum fixnum fixnum fixnum)))
+ (make-instance (class-of matrix)
+ :nrows nc :ncols nr
+ :store (store matrix)
+ :head hd
+ :row-stride cs :col-stride rs)))
+ ((typep matrix 'number) matrix)
+ (t (error "Don't know how to take the transpose of ~A." matrix))))
+
+;;
+(defmacro with-transpose! (matlst &rest body)
+ `(progn
+ ,@(mapcar #'(lambda (mat) `(transpose-i! ,mat)) matlst)
+ ,@body
+ ,@(mapcar #'(lambda (mat) `(transpose-i! ,mat)) matlst)))
+
+;;
+(defun sub! (matrix i j nrows ncols)
(declare (type standard-matrix matrix)
(type fixnum i j nrows ncols))
(let ((hd (head matrix))
@@ -286,20 +318,49 @@ matrix and a number"))
(rs (row-stride matrix))
(cs (col-stride matrix)))
(declare (type fixnum hd nr nc rs cs))
- (when (or (> i (- nr nrows)) (> j (- nc ncols)))
- (error "Parent matrix is not big enough to create sub-matrix."))
+ (unless (and (< -1 i (+ i nrows) nr) (< -1 j (+ j ncols) nc))
+ (error "Bad index and/or size.
+Cannot create a sub-matrix of size (~a ~a) starting at (~a ~a)" nrows ncols i j))
(make-instance (class-of matrix)
:nrows nrows :ncols ncols
:store (store matrix)
:head (store-indexing i j hd rs cs)
:row-stride rs :col-stride cs)))
+(defun (setf sub!) (mat-b mat-a i j nrows ncols)
+ (copy! mat-b (sub! mat-a i j nrows ncols)))
+
;;
-(defmacro with-transpose! (matlst &rest body)
- `(progn
- ,@(mapcar #'(lambda (mat) `(transpose! ,mat)) matlst)
- ,@body
- ,@(mapcar #'(lambda (mat) `(transpose! ,mat)) matlst)))
+(defun row! (matrix i)
+ (declare (type standard-matrix matrix)
+ (type fixnum i))
+ (mlet* (((hd nr nc rs cs) (slot-values matrix '(head number-of-rows number-of-cols row-stride col-stride)) :type (fixnum fixnum fixnum fixnum)))
+ (unless (< -1 i nr)
+ (error "Index ~a is outside the valid range for the given matrix." i))
+ (make-instance (class-of matrix)
+ :nrows 1 :ncols nc
+ :store (store matrix)
+ :head (store-indexing i 0 hd rs cs)
+ :row-stride rs :col-stride cs)))
+
+(defun (setf row!) (mat-b mat-a i)
+ (copy! mat-b (row! mat-a i)))
+
+;;
+(defun col! (matrix j)
+ (declare (type standard-matrix matrix)
+ (type fixnum j))
+ (mlet* (((hd nr nc rs cs) (slot-values matrix '(head number-of-rows number-of-cols row-stride col-stride)) :type (fixnum fixnum fixnum fixnum)))
+ (unless (< -1 j nc)
+ (error "Index ~a is outside the valid range for the given matrix." j))
+ (make-instance (class-of matrix)
+ :nrows nr :ncols 1
+ :store (store matrix)
+ :head (store-indexing 0 j hd rs cs)
+ :row-stride rs :col-stride cs)))
+
+(defun (setf col!) (mat-b mat-a j)
+ (copy! mat-b (col! mat-a j)))
;;
(defun blas-copyable-p (matrix)
@@ -325,6 +386,8 @@ matrix and a number"))
((= cs 1) (values :row-major rs (cond
((string= fortran-op "N" ) "T")
((string= fortran-op "T" ) "N"))))
- ((= rs 1) (values :col-major cs fortran-op))
+ ((= rs 1) (values :col-major cs (cond
+ ((string= fortran-op "N" ) "N")
+ ((string= fortran-op "N" ) "T"))))
;;Lets not confound lisp's type declaration.
(t (values nil -1 "?")))))
\ No newline at end of file
-----------------------------------------------------------------------
Summary of changes:
packages.lisp | 6 +-
src/axpy.lisp | 26 ++---
src/blas.lisp | 32 +++----
src/complex-matrix.lisp | 18 +++-
src/copy.lisp | 1 -
src/gemv.lisp | 258 +++++++++++++++++++++++++++++++++++++++++-----
src/standard-matrix.lisp | 119 ++++++++++++++++-----
src/utilities.lisp | 28 +++--
8 files changed, 384 insertions(+), 104 deletions(-)
hooks/post-receive
--
matlisp
|
|
From: Raymond T. <rt...@us...> - 2012-03-15 05:57:54
|
This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "matlisp".
The branch, master has been updated
via 55a0678924a585263937fdff69c5d522f1589d7a (commit)
via ced92c7a0cbf00be10fd0dff010a4768130da9c1 (commit)
via 74e3f184bf5e3b0974d145a823b4f86e33e00a1a (commit)
via 45962c22d9931aa0961dca11931b50d884e3ec01 (commit)
from 8a2722a7452cac1de3edd5fc8cc65d97c0d1a0a8 (commit)
Those revisions listed above that are new to this repository have
not appeared on any other notification email; so we list those
revisions in full, below.
- Log -----------------------------------------------------------------
commit 55a0678924a585263937fdff69c5d522f1589d7a
Author: Raymond Toy <toy...@gm...>
Date: Wed Mar 14 22:57:39 2012 -0700
Regenerated.
diff --git a/configure b/configure
index 35ee372..2c91c83 100755
--- a/configure
+++ b/configure
@@ -15236,6 +15236,82 @@ else
fi
+# Check to see if the ATLAS libraries are compatible with matlisp's
+# ffi. Basically the same test as above that checks to see if -ff2c
+# is needed. We call zdotu which is a Fortran function returning a
+# complex number. Matlisp assumes such functions return the result by
+# storing the answer at address given by a hidden first parameter to
+# the function.
+
+if test x"$atlas" = xtrue; then
+ { $as_echo "$as_me:${as_lineno-$LINENO}: checking if ATLAS is compatible with f2c calling conventions" >&5
+$as_echo_n "checking if ATLAS is compatible with f2c calling conventions... " >&6; }
+ # From the value of f77_name, figure out the actual name for
+ # Fortran's zdotu.
+ case $f77_name in
+ F77*) case $f77_name in
+ F77_NAME) ZDOTU="ZDOTU" ;;
+ F77_NAME_) ZDOTU="ZDOTU_" ;;
+ F77_NAME__) ZDOTU="ZDOTU_" ;;
+ esac
+ ;;
+ f77*) case $f77_name in
+ f77_name) ZDOTU="zdotu" ;;
+ f77_name_) ZDOTU="zdotu_" ;;
+ f77_name__) ZDOTU="zdotu_" ;;
+ esac
+ ;;
+ esac
+
+ cat > conftest.c <<EOF
+
+#include <stdio.h>
+#include <math.h>
+
+extern void ${ZDOTU}(double *, int *, double *, int *, double *, int *);
+
+int main()
+{
+ int rc;
+
+ int n = 2;
+ int incx = 1;
+ double x[4];
+ double out[2];
+ x[0] = 1;
+ x[1] = 0;
+ x[2] = 2;
+ x[3] = 0;
+ out[0] = 0;
+ out[1] = 0;
+
+ ${ZDOTU}(out, &n, x, &incx, x, &incx);
+
+ if (fabs(out[0] - 5) < 1e-10) {
+ rc = 0;
+ } else {
+ printf("Actual output = %lg, %lg, instead of 5, 0\n", out[0], out[1]);
+ rc = 1;
+ }
+
+ return rc;
+}
+
+EOF
+ $CC $CFLAGS -c conftest.c
+ $F77 $FFLAGS -o a.out conftest.o -L${ATLAS_DIR} -latlas -lcblas -lf77blas -llapack
+ if a.out; then
+ { $as_echo "$as_me:${as_lineno-$LINENO}: result: yes" >&5
+$as_echo "yes" >&6; }
+ F2C=-ff2c
+ else
+ { $as_echo "$as_me:${as_lineno-$LINENO}: result: no" >&5
+$as_echo "no" >&6; }
+ F2C=""
+ fi
+fi
+
+# The following variables will be substituted into the .in files
commit ced92c7a0cbf00be10fd0dff010a4768130da9c1
Author: Raymond Toy <toy...@gm...>
Date: Wed Mar 14 22:57:28 2012 -0700
Put back the check for ATLAS compatible with f2c conventions. We need
to know this to compile the libmatisp compatibility function correctly
because it calls zdotu/c.
diff --git a/configure.ac b/configure.ac
index ef543fb..9bbb440 100644
--- a/configure.ac
+++ b/configure.ac
@@ -314,78 +314,79 @@ AC_HELP_STRING([--with-atlas=libpath], [Location of the ATLAS libraries]),
])
AM_CONDITIONAL([ATLAS], [test x$atlas = xtrue])
-dnl # Check to see if the ATLAS libraries are compatible with matlisp's
-dnl # ffi. Basically the same test as above that checks to see if -ff2c
-dnl # is needed. We call zdotu which is a Fortran function returning a
-dnl # complex number. Matlisp assumes such functions return the result by
-dnl # storing the answer at address given by a hidden first parameter to
-dnl # the function.
-dnl
-dnl if test x"$atlas" = xtrue; then
-dnl AC_MSG_CHECKING([if ATLAS is compatible with matlisp's FFI])
-dnl # From the value of f77_name, figure out the actual name for
-dnl # Fortran's zdotu.
-dnl case $f77_name in
-dnl F77*) case $f77_name in
-dnl F77_NAME) ZDOTU="ZDOTU" ;;
-dnl F77_NAME_) ZDOTU="ZDOTU_" ;;
-dnl F77_NAME__) ZDOTU="ZDOTU_" ;;
-dnl esac
-dnl ;;
-dnl f77*) case $f77_name in
-dnl f77_name) ZDOTU="zdotu" ;;
-dnl f77_name_) ZDOTU="zdotu_" ;;
-dnl f77_name__) ZDOTU="zdotu_" ;;
-dnl esac
-dnl ;;
-dnl esac
-dnl
-dnl cat > conftest.c <<EOF
-dnl [
-dnl #include <stdio.h>
-dnl #include <math.h>
-dnl
-dnl extern void ${ZDOTU}(double *, int *, double *, int *, double *, int *);
-dnl
-dnl int main()
-dnl {
-dnl int rc;
-dnl
-dnl int n = 2;
-dnl int incx = 1;
-dnl double x[4];
-dnl double out[2];
-dnl x[0] = 1;
-dnl x[1] = 0;
-dnl x[2] = 2;
-dnl x[3] = 0;
-dnl out[0] = 0;
-dnl out[1] = 0;
-dnl
-dnl ${ZDOTU}(out, &n, x, &incx, x, &incx);
-dnl
-dnl if (fabs(out[0] - 5) < 1e-10) {
-dnl rc = 0;
-dnl } else {
-dnl printf("Actual output = %lg, %lg, instead of 5, 0\n", out[0], out[1]);
-dnl rc = 1;
-dnl }
-dnl
-dnl return rc;
-dnl }
-dnl ]
-dnl EOF
-dnl $CC $CFLAGS -c conftest.c
-dnl $F77 $FFLAGS -o a.out conftest.o -L${ATLAS_DIR} -latlas -lcblas -lf77blas -llapack
-dnl if a.out; then
-dnl AC_MSG_RESULT([yes])
-dnl else
-dnl AC_MSG_RESULT([no])
-dnl AC_MSG_ERROR([ATLAS libraries are not compatible with matlisp.])
-dnl fi
-dnl fi
-dnl
-dnl The following variables will be substituted into the .in files
+# Check to see if the ATLAS libraries are compatible with matlisp's
+# ffi. Basically the same test as above that checks to see if -ff2c
+# is needed. We call zdotu which is a Fortran function returning a
+# complex number. Matlisp assumes such functions return the result by
+# storing the answer at address given by a hidden first parameter to
+# the function.
+
+if test x"$atlas" = xtrue; then
+ AC_MSG_CHECKING([if ATLAS is compatible with f2c calling conventions])
+ # From the value of f77_name, figure out the actual name for
+ # Fortran's zdotu.
+ case $f77_name in
+ F77*) case $f77_name in
+ F77_NAME) ZDOTU="ZDOTU" ;;
+ F77_NAME_) ZDOTU="ZDOTU_" ;;
+ F77_NAME__) ZDOTU="ZDOTU_" ;;
+ esac
+ ;;
+ f77*) case $f77_name in
+ f77_name) ZDOTU="zdotu" ;;
+ f77_name_) ZDOTU="zdotu_" ;;
+ f77_name__) ZDOTU="zdotu_" ;;
+ esac
+ ;;
+ esac
+
+ cat > conftest.c <<EOF
+[
+#include <stdio.h>
+#include <math.h>
+
+extern void ${ZDOTU}(double *, int *, double *, int *, double *, int *);
+
+int main()
+{
+ int rc;
+
+ int n = 2;
+ int incx = 1;
+ double x[4];
+ double out[2];
+ x[0] = 1;
+ x[1] = 0;
+ x[2] = 2;
+ x[3] = 0;
+ out[0] = 0;
+ out[1] = 0;
+
+ ${ZDOTU}(out, &n, x, &incx, x, &incx);
+
+ if (fabs(out[0] - 5) < 1e-10) {
+ rc = 0;
+ } else {
+ printf("Actual output = %lg, %lg, instead of 5, 0\n", out[0], out[1]);
+ rc = 1;
+ }
+
+ return rc;
+}
+]
+EOF
+ $CC $CFLAGS -c conftest.c
+ $F77 $FFLAGS -o a.out conftest.o -L${ATLAS_DIR} -latlas -lcblas -lf77blas -llapack
+ if a.out; then
+ AC_MSG_RESULT([yes])
+ F2C=-ff2c
+ else
+ AC_MSG_RESULT([no])
+ F2C=""
+ fi
+fi
+
+# The following variables will be substituted into the .in files
AC_SUBST(LISPEXEC)
AC_SUBST(BLAS_OBJS)
AC_SUBST(NO_ATLAS_LAPACK_OBJS)
@@ -404,7 +405,7 @@ AC_SUBST(FLIBS_OPTS)
AC_SUBST(share_ext)
AC_SUBST(GIT_VERSION)
AC_SUBST(HAVE_QL)
-dnl AC_SUBST(F2C)
+AC_SUBST(F2C)
echo host = $host
commit 74e3f184bf5e3b0974d145a823b4f86e33e00a1a
Author: Raymond Toy <toy...@gm...>
Date: Wed Mar 14 22:56:23 2012 -0700
Load matlisp after loading all the rest of the libraries.
diff --git a/lib/lazy-loader.lisp.in b/lib/lazy-loader.lisp.in
index d558183..ba1296c 100644
--- a/lib/lazy-loader.lisp.in
+++ b/lib/lazy-loader.lisp.in
@@ -180,8 +180,8 @@
)
(t
(cffi:use-foreign-library blas)
- (cffi:use-foreign-library matlisp)
- (cffi:use-foreign-library lapack))))
+ (cffi:use-foreign-library lapack)
+ (cffi:use-foreign-library matlisp))))
#+:allegro
(defun load-blas-&-lapack-libraries ()
commit 45962c22d9931aa0961dca11931b50d884e3ec01
Author: Raymond Toy <toy...@gm...>
Date: Wed Mar 14 22:55:54 2012 -0700
Don't need the dependency on libblas.
diff --git a/lib-src/compat/Makefile.am b/lib-src/compat/Makefile.am
index afa1180..95a71cd 100644
--- a/lib-src/compat/Makefile.am
+++ b/lib-src/compat/Makefile.am
@@ -5,8 +5,8 @@ if LIB32
AM_FFLAGS += -m32
endif
-libmatlisp_la_LDFLAGS = ../../LAPACK/BLAS/SRC/libblas.la $(FLIBS_OPTS)
-libmatlisp_la_LIBADD = $(FLIBS)
+libmatlisp_la_LDFLAGS = $(FLIBS_OPTS)
+libmatlisp_la_LIBADD = $(FLIBS_LIBS)
libmatlisp_la_SOURCES = \
compat.f
-----------------------------------------------------------------------
Summary of changes:
configure | 76 +++++++++++++++++++++++
configure.ac | 147 ++++++++++++++++++++++----------------------
lib-src/compat/Makefile.am | 4 +-
lib/lazy-loader.lisp.in | 4 +-
4 files changed, 154 insertions(+), 77 deletions(-)
hooks/post-receive
--
matlisp
|
|
From: Akshay S. <ak...@us...> - 2012-03-14 05:58:28
|
This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "matlisp".
The branch, matlisp-cffi has been updated
via 517949e31b41b303dab91670e80b207bc45d3256 (commit)
via 452bd980d066101fad780815ff5ddfccdcf5e683 (commit)
via f4b94df81eab2264c071a5a6f592e6becf76f770 (commit)
via 8a2722a7452cac1de3edd5fc8cc65d97c0d1a0a8 (commit)
via ecab7fc4f75f615eabed198e0e7a325b2c055fd0 (commit)
via 14707c10ced9ee08431313323cc1905ac77e2a21 (commit)
via 7f45afca17c2501c0b2a777809e1998683e5a17f (commit)
via 8708b5e4f149f2ba16b879bfa87aa3085097cb34 (commit)
via 4d23b62fbced07e732614e290a08e90a731942ab (commit)
via 579b17087f685d198095ef8616a6374ed59d6848 (commit)
via ec16280b56855cff043ec1fc6bbddc11cf2cebca (commit)
via 6695b66525322e6817de899b467fb3505ec22775 (commit)
via b4da181a45007d71d24371851089844b257e2fa8 (commit)
from 2cbbcb64c997e7ed95c2e9a344347ca9606fce66 (commit)
Those revisions listed above that are new to this repository have
not appeared on any other notification email; so we list those
revisions in full, below.
- Log -----------------------------------------------------------------
commit 517949e31b41b303dab91670e80b207bc45d3256
Author: Akshay Srinivasan <aks...@gm...>
Date: Wed Mar 14 11:24:44 2012 +0530
-> gemv core function works.
diff --git a/src/copy.lisp b/src/copy.lisp
index cf987d9..03c6e8e 100644
--- a/src/copy.lisp
+++ b/src/copy.lisp
@@ -80,6 +80,10 @@
;;
(defmacro generate-typed-copy!-func (func store-type matrix-type blas-func)
+ ;;Be very careful when using functions generated by this macro.
+ ;;Indexes can be tricky and this has no safety net
+ ;;(you don't see a matrix-ref do you ?)
+ ;;Use only after checking the arguments for compatibility.
`(defun ,func (mat-a mat-b)
(declare (type ,matrix-type mat-a mat-b)
(optimize (safety 0) (speed 3)))
@@ -89,24 +93,26 @@
((hd-b st-b) (slot-values mat-b '(head store)) :type (fixnum (,store-type *))))
(if (and cp-a cp-b)
(,blas-func sz st-a inc-a st-b inc-b :head-x hd-a :head-y hd-b)
- (symbol-macrolet
- ((common-code
- (mlet* (((nr-a nc-a rs-a cs-a) (slot-values mat-a '(number-of-rows number-of-cols row-stride col-stride))
- :type (fixnum fixnum fixnum fixnum))
- ((rs-b cs-b) (slot-values mat-b '(row-stride col-stride))
- :type (fixnum fixnum)))
- (loop for i from 0 below nr-a
- do (,blas-func nc-a st-a cs-a st-b cs-b :head-x (+ hd-a (* i rs-a)) :head-y (+ hd-b (* i rs-b)))))))
- ;;Choose the smaller of the loops
- (if (> (nrows mat-a) (ncols mat-a))
- (with-transpose! (mat-a mat-b)
- common-code)
- common-code)))
- mat-b)))
+ (mlet* (((nr-a nc-a rs-a cs-a) (slot-values mat-a '(number-of-rows number-of-cols row-stride col-stride))
+ :type (fixnum fixnum fixnum fixnum))
+ ((rs-b cs-b) (slot-values mat-b '(row-stride col-stride))
+ :type (fixnum fixnum)))
+ ;;Choose the smaller of the loops
+ (when (> (nrows mat-a) (ncols mat-a))
+ (rotatef nr-a nc-a)
+ (rotatef rs-a cs-a)
+ (rotatef rs-b cs-b))
+ (loop for i from 0 below nr-a
+ do (,blas-func nc-a st-a cs-a st-b cs-b :head-x (+ hd-a (* i rs-a)) :head-y (+ hd-b (* i rs-b)))))))
+ mat-b))
;;
(defmacro generate-typed-num-copy!-func (func element-type store-type matrix-type blas-func
array-decl)
+ ;;Be very careful when using functions generated by this macro.
+ ;;Indexes can be tricky and this has no safety net
+ ;;(you don't see a matrix-ref do you ?)
+ ;;Use only after checking the arguments for compatibility.
(let ((num-arg (car array-decl)))
(destructuring-bind (var var-maker-form var-setf-form &key type) (cadr array-decl)
`(mlet* (((,var) (,@var-maker-form) ,@(if type
@@ -120,18 +126,15 @@
((hd-x st-x sz) (slot-values mat-x '(head store number-of-elements)) :type (fixnum (,store-type *) fixnum)))
(if cp-x
(,blas-func sz ,var 0 st-x inc-x :head-y hd-x)
- (symbol-macrolet
- ((common-code
- (mlet* (((nr-x nc-x rs-x cs-x) (slot-values mat-x '(number-of-rows number-of-cols row-stride col-stride))
- :type (fixnum fixnum fixnum fixnum)))
- (loop for i from 0 below nr-x
- do (,blas-func nc-x ,var 0 st-x cs-x :head-y (+ hd-x (* i rs-x)))))))
- ;;Choose the smaller of the loops
- (if (> (nrows mat-x) (ncols mat-x))
- (with-transpose! (mat-x)
- common-code)
- common-code)))
- mat-x))))))
+ (mlet* (((nr-x nc-x rs-x cs-x) (slot-values mat-x '(number-of-rows number-of-cols row-stride col-stride))
+ :type (fixnum fixnum fixnum fixnum)))
+ ;;Choose the smaller of the loops
+ (when (> (nrows mat-x) (ncols mat-x))
+ (rotatef nr-x nc-x)
+ (rotatef rs-x cs-x))
+ (loop for i from 0 below nr-x
+ do (,blas-func nc-x ,var 0 st-x cs-x :head-y (+ hd-x (* i rs-x)))))))
+ mat-x)))))
;;
(defgeneric copy! (matrix new-matrix)
diff --git a/src/gemv.lisp b/src/gemv.lisp
index 9adb6cc..cffa7d6 100644
--- a/src/gemv.lisp
+++ b/src/gemv.lisp
@@ -1,35 +1,32 @@
(in-package :matlisp)
(defun real-double-gemv!-typed (alpha A x beta y job)
- (declare ((type double-float alpha beta)
- (type real-matrix A x y)
- (type symbol job)))
- (symbol-macrolet
- ((common-block
- (mlet* (((st-a hd-a nr-a nc-a rs-a cs-a) (slot-values A '(store head number-of-rows number-of-cols row-stride col-stride))
- :type ((real-matrix-store-type *) fixnum fixnum fixnum fixnum fixnum))
- ((st-x hd-x rs-x) (slot-values x '(store head row-stride))
- :type ((real-matrix-store-type *) fixnum fixnum))
- ((st-y hd-y rs-y) (slot-values y '(store head row-stride))
- :type ((real-matrix-store-type *) fixnum fixnum))
- ((sym lda tf-job) (blas-matrix-compatible-p A) :type (nil fixnum (string 1))))
- (if (not (string= tf-job "?"))
- (blas:dgemv tf-job nr-a nc-a alpha st-a lda st-x rs-x beta st-y rs-y :head-a hd-a :head-x hd-x :head-y hd-y)
- (dotimes (i nr-a)
- (setf (matrix-ref-2d y i 0) (+ (* alpha (blas:ddot nc-a st-a cs-a st-x rs-x :head-x (+ hd-a (* i rs-a)) :head-y hd-x))
- (* beta (matrix-ref-2d y i 0))))
-
- (mlet* ((fortran-job (ecase job (:n "N") (:t "T")) :type ((string 1)))
+ ;;Be very careful when you use this function!
+ ;;Indexes can be tricky and this has no safety net
+ ;;(you don't see a matrix-ref do you ?)
+ ;;Use only after checking the arguments for compatibility.
+ (declare (type double-float alpha beta)
+ (type real-matrix A x y)
+ (type symbol job))
+ (mlet* ((fort-op (ecase job (:n "N") (:t "T")) :type ((string 1)))
((st-a hd-a nr-a nc-a rs-a cs-a) (slot-values A '(store head number-of-rows number-of-cols row-stride col-stride))
:type ((real-matrix-store-type *) fixnum fixnum fixnum fixnum fixnum))
((st-x hd-x rs-x) (slot-values x '(store head row-stride))
:type ((real-matrix-store-type *) fixnum fixnum))
((st-y hd-y rs-y) (slot-values y '(store head row-stride))
:type ((real-matrix-store-type *) fixnum fixnum))
- ((sym lda tf-job) (blas-matrix-compatible-p A fortran-job) :type (nil fixnum (string 1))))
- (if (not (string= tf-job "?"))
- (blas:dgemv tf-job nr-a nc-a alpha st-a lda st-x rs-x beta st-y rs-y :head-a hd-a :head-x hd-x :head-y hd-y)
- (dotimes (i nr-a)
- (setf (matrix-ref-2d y i 0) (+ (* alpha (blas:ddot nc-a st-a cs-a st-x rs-x :head-x (+ hd-a (* i rs-a)) :head-y hd-x))
- (* beta (matrix-ref-2d y i 0))))
-
\ No newline at end of file
+ ((sym lda tf-op) (blas-matrix-compatible-p A fort-op) :type (symbol fixnum (string 1))))
+ (if (not (string= tf-op "?"))
+ (progn
+ (when (eq sym :row-major)
+ (rotatef nr-a nc-a)
+ (rotatef rs-a cs-a))
+ (blas:dgemv tf-op nr-a nc-a alpha st-a lda st-x rs-x beta st-y rs-y :head-a hd-a :head-x hd-x :head-y hd-y))
+ (progn
+ (when (string= fort-op "T")
+ (rotatef nr-a nc-a)
+ (rotatef rs-a cs-a))
+ (dotimes (i nr-a)
+ (setf (matrix-ref-2d y i 0) (+ (* alpha (blas:ddot nc-a st-a cs-a st-x rs-x :head-x (+ hd-a (* i rs-a)) :head-y hd-x))
+ (* beta (matrix-ref-2d y i 0))))))))
+ y)
\ No newline at end of file
diff --git a/src/scal.lisp b/src/scal.lisp
index 88da573..a81af20 100644
--- a/src/scal.lisp
+++ b/src/scal.lisp
@@ -70,6 +70,10 @@
(in-package "MATLISP")
(defmacro generate-typed-scal!-func (func element-type store-type matrix-type blas-func)
+ ;;Be very careful when using functions generated by this macro.
+ ;;Indexes can be tricky and this has no safety net
+ ;;(you don't see a matrix-ref do you ?)
+ ;;Use only after checking the arguments for compatibility.
`(defun ,func (alpha mat-x)
(declare (type ,matrix-type mat-x)
(type ,element-type alpha)
@@ -80,17 +84,15 @@
:type (fixnum (,store-type *))))
(if cp-x
(,blas-func sz-x alpha st-x inc-x :head-x hd-x)
- (symbol-macrolet
- ((common-code
- (mlet* (((nr-x nc-x rs-x cs-x) (slot-values mat-x '(number-of-rows number-of-cols row-stride col-stride))
- :type (fixnum fixnum fixnum fixnum)))
- (loop for i from 0 below nr-x
- do (,blas-func nc-x alpha st-x cs-x :head-x (+ hd-x (* i rs-x)))))))
- (if (> (nrows mat-x) (ncols mat-x))
- (with-transpose! (mat-x)
- common-code)
- common-code)))
- mat-x)))
+ (mlet* (((nr-x nc-x rs-x cs-x) (slot-values mat-x '(number-of-rows number-of-cols row-stride col-stride))
+ :type (fixnum fixnum fixnum fixnum)))
+ ;;Choose the smaller of the loops.
+ (when (> (nrows mat-x) (ncols mat-x))
+ (rotatef nr-x nc-x)
+ (rotatef rs-x cs-x))
+ (loop for i from 0 below nr-x
+ do (,blas-func nc-x alpha st-x cs-x :head-x (+ hd-x (* i rs-x)))))))
+ mat-x))
;;
(defgeneric scal! (alpha x)
diff --git a/src/standard-matrix.lisp b/src/standard-matrix.lisp
index 3c2cfce..0f1ec71 100644
--- a/src/standard-matrix.lisp
+++ b/src/standard-matrix.lisp
@@ -318,8 +318,8 @@ matrix and a number"))
;;
(defun blas-matrix-compatible-p (matrix &optional (fortran-op "N"))
(declare (optimize (safety 0) (speed 3))
- (type (or real-matrix complex-matrix) mat))
- (mlet* (((rs cs) (slot-values mat '(row-stride col-stride))
+ (type (or real-matrix complex-matrix) matrix))
+ (mlet* (((rs cs) (slot-values matrix '(row-stride col-stride))
:type (fixnum fixnum)))
(cond
((= cs 1) (values :row-major rs (cond
commit 452bd980d066101fad780815ff5ddfccdcf5e683
Author: Akshay Srinivasan <aks...@gm...>
Date: Tue Mar 13 22:03:18 2012 +0530
-> Ported axpy to support new classes.
Started work on Level-2 stuff. Does not work yet.
diff --git a/src/axpy.lisp b/src/axpy.lisp
index 57104cf..d6cee37 100644
--- a/src/axpy.lisp
+++ b/src/axpy.lisp
@@ -74,23 +74,78 @@
(in-package "MATLISP")
-#+nil (use-package "LAPACK")
-#+nil (use-package "BLAS")
-#+nil (use-package "FORTRAN-FFI-ACCESSORS")
+(defmacro generate-typed-axpy!-func (func element-type store-type matrix-type blas-func)
+ `(defun ,func (alpha mat-a mat-b)
+ (declare (type ,element-type alpha)
+ (type ,matrix-type mat-a mat-b)
+ (optimize (safety 0) (speed 3)))
+ (mlet* (((cp-a inc-a sz-a) (blas-copyable-p mat-a) :type (boolean fixnum nil))
+ ((cp-b inc-b sz-b) (blas-copyable-p mat-b) :type (boolean fixnum nil))
+ ((hd-a st-a sz) (slot-values mat-a '(head store number-of-elements)) :type (fixnum (,store-type *) fixnum))
+ ((hd-b st-b) (slot-values mat-b '(head store)) :type (fixnum (,store-type *))))
+ (if (and cp-a cp-b)
+ (,blas-func sz alpha st-a inc-a st-b inc-b :head-x hd-a :head-y hd-b)
+ (symbol-macrolet
+ ((common-code
+ (mlet* (((nr-a nc-a rs-a cs-a) (slot-values mat-a '(number-of-rows number-of-cols row-stride col-stride))
+ :type (fixnum fixnum fixnum fixnum))
+ ((rs-b cs-b) (slot-values mat-b '(row-stride col-stride))
+ :type (fixnum fixnum)))
+ (loop for i from 0 below nr-a
+ do (,blas-func nc-a alpha st-a cs-a st-b cs-b :head-x (+ hd-a (* i rs-a)) :head-y (+ hd-b (* i rs-b)))))))
+ ;;Choose the smaller of the loops
+ (if (> (nrows mat-a) (ncols mat-a))
+ (with-transpose! (mat-a mat-b)
+ common-code)
+ common-code)))
+ mat-b)))
-#+nil (export '(axpy!
- axpy))
+;;
+(defgeneric axpy! (alpha x y)
+ (:documentation
+ "
+ Syntax
+ ======
+ (AXPY! alpha x y)
+
+ Y <- alpha * x + y
+
+ Purpose
+ =======
+ Same as AXPY except that the result
+ is stored in Y and Y is returned.
+"))
+
+(defmethod axpy! :before ((alpha number) (x standard-matrix) (y standard-matrix))
+ (mlet* (((nr-x nc-x) (slot-values x '(number-of-rows number-of-cols)) :type (fixnum fixnum))
+ ((nr-y nc-y) (slot-values y '(number-of-rows number-of-cols)) :type (fixnum fixnum)))
+ (unless (and (= nr-x nr-y) (= nc-x nc-y))
+ (error "Arguments X,Y to AXPY! are of different dimensions."))))
-;; note: we should optimize the calls to axpy from other axpy's since
-;; we know exactly which one we are calling.
;;
-;; also, the most common type of bug is with ! operators, e.g. when
-;; you say (axpy! s a a)
+(generate-typed-axpy!-func real-double-axpy!-typed double-float real-matrix-store-type real-matrix blas:daxpy)
+
+(defmethod axpy! ((alpha number) (x complex-matrix) (y real-matrix))
+ (error "cannot AXPY! a complex X to a real Y,
+don't know how to coerce COMPLEX to REAL"))
-(deftype complex-double-float ()
- '(cl:complex (double-float * *)))
+(defmethod axpy! ((alpha cl:real) (x real-matrix) (y real-matrix))
+ (real-double-axpy!-typed (coerce alpha 'double-float) x y))
;;
+(generate-typed-axpy!-func complex-double-axpy!-typed complex-double-float complex-matrix-store-type complex-matrix blas:zaxpy)
+
+(defmethod axpy! ((alpha cl:real) (x real-matrix) (y complex-matrix))
+ (real-double-axpy!-typed (coerce alpha 'double-float) x (realpart! y)))
+
+(defmethod axpy! ((alpha complex) (x real-matrix) (y complex-matrix))
+ (real-double-axpy!-typed (coerce (realpart alpha) 'double-float) x (realpart! y))
+ (real-double-axpy!-typed (coerce (imagpart alpha) 'double-float) x (imagpart! y)))
+
+(defmethod axpy! ((alpha number) (x complex-matrix) (y complex-matrix))
+ (complex-double-axpy!-typed (complex-coerce alpha) x y))
+
+;;;;
(defgeneric axpy (alpha x y)
(:documentation
"
@@ -118,187 +173,29 @@
"))
(defmethod axpy :before ((alpha number) (x standard-matrix) (y standard-matrix))
- (let ((nxm-x (number-of-elements x))
- (nxm-y (number-of-elements y)))
- (declare (type fixnum nxm-x nxm-y))
-
- (if (not (= nxm-x nxm-y))
- (error "arguments X and Y to AXPY not of the same size"))))
+ (mlet* (((nr-x nc-x) (slot-values x '(number-of-rows number-of-cols)) :type (fixnum fixnum))
+ ((nr-y nc-y) (slot-values y '(number-of-rows number-of-cols)) :type (fixnum fixnum)))
+ (unless (and (= nr-x nr-y) (= nc-x nc-y))
+ (error "Arguments X,Y to AXPY are of different dimensions."))))
;;
(defmethod axpy ((alpha cl:real) (x real-matrix) (y real-matrix))
- (axpy (coerce alpha 'real-matrix-element-type) x y))
-
-(defmethod axpy ((alpha #+(or cmu sbcl) double-float #-(or cmu sbcl) float) (x real-matrix) (y real-matrix))
- (let* ((nxm (number-of-elements y))
- (result (copy y)))
- (declare (type fixnum nxm))
+ (let ((result (copy y)))
+ (axpy! alpha x result)))
- (daxpy nxm alpha (store x) 1 (store result) 1)
- result))
-
-;;
-(defmethod axpy ((alpha cl:real) (x complex-matrix) (y real-matrix))
- (axpy (coerce alpha 'real-matrix-element-type) x y))
+(defmethod axpy ((alpha complex) (x real-matrix) (y real-matrix))
+ (let ((result (scal alpha x)))
+ (axpy! 1d0 y result)))
-(defmethod axpy ((alpha #+(or cmu sbcl) double-float #-(or cmu sbcl) float) (x complex-matrix) (y real-matrix))
- (let* ((nxm (number-of-elements y))
- (n (nrows y))
- (m (ncols y))
- (result (make-complex-matrix-dim n m))
- (store-x (store x))
- (store-y (store y))
- (store-result (store result)))
- (declare (type fixnum n m nxm)
- (type (real-matrix-store-type *) store-y)
- (type (complex-matrix-store-type *) store-x store-result))
-
- (zcopy nxm store-x 1 store-result 1) ;; same as (COPY! x result)
- (zdscal nxm alpha store-result 1) ;; same as (SCAL! alpha result)
- (daxpy nxm 1.0d0 store-y 1 store-result 2) ;; same as (AXPY! 1d0 y result)
- result))
-
-;;
-(defmethod axpy ((alpha cl:real) (x real-matrix) (y complex-matrix))
- (axpy (coerce alpha 'complex-matrix-element-type) x y))
-
-(defmethod axpy ((alpha #+(or cmu sbcl) double-float #-(or cmu sbcl) float) (x real-matrix) (y complex-matrix))
- (let* ((nxm (number-of-elements y))
- (result (copy y)))
- (declare (type fixnum nxm))
- (daxpy nxm alpha (store x) 1 (store result) 2)
- result))
-
-;;
-(defmethod axpy ((alpha cl:real) (x complex-matrix) (y complex-matrix))
- (axpy (coerce alpha 'complex-matrix-element-type) x y))
-
-(defmethod axpy ((alpha #+(or cmu sbcl) double-float #-(or cmu sbcl) float) (x complex-matrix) (y complex-matrix))
- (let ((nxm (number-of-elements y))
- (result (copy y)))
- (declare (type fixnum nxm))
- (daxpy (* 2 nxm) alpha (store x) 1 (store result) 1)
- result))
-
-;;
-(defmethod axpy ((alpha number) (x real-matrix) (y complex-matrix))
- (let* ((nxm (number-of-elements y))
- (n (nrows y))
- (m (ncols y))
- (result (make-complex-matrix-dim n m))
- (store-x (store x))
- (store-y (store y))
- (store-result (store result))
- (c-alpha (complex-coerce alpha)))
- (declare (type complex-double-float c-alpha)
- (type fixnum n m nxm)
- (type (real-matrix-store-type *) store-x)
- (type (complex-matrix-store-type *) store-y store-result))
-
- (dcopy nxm store-x 1 store-result 2)
- (zscal nxm c-alpha store-result 1)
- (daxpy (* 2 nxm) 1.0d0 store-y 1 store-result 1)
-
- result))
-
-;;
(defmethod axpy ((alpha number) (x complex-matrix) (y real-matrix))
- (let* ((nxm (number-of-elements y))
- (result (copy x))
- (store-result (store result))
- (c-alpha (complex-coerce alpha)))
- (declare (type complex-double-float c-alpha)
- (type fixnum nxm)
- (type (complex-matrix-store-type *) store-result))
-
- (zscal nxm c-alpha store-result 1)
- (daxpy nxm 1.0d0 (store y) 1 store-result 2)
-
- result))
-
-;;
-(defmethod axpy ((alpha number) (x complex-matrix) (y complex-matrix))
- (let ((nxm (number-of-elements y))
- (result (copy y))
- (c-alpha (complex-coerce alpha)))
- (declare (type complex-double-float c-alpha)
- (type fixnum nxm))
-
- (zaxpy nxm c-alpha (store x) 1 (store result) 1)
-
- result))
+ (let ((result (scal alpha x)))
+ (axpy! 1d0 y result)))
;;
-(defgeneric axpy! (alpha x y)
- (:documentation
- "
- Syntax
- ======
- (AXPY! alpha x y)
-
- Purpose
- =======
- Same as AXPY except that the result
- is stored in Y and Y is returned.
-"))
-
-(defmethod axpy! :before ((alpha number) (x standard-matrix) (y standard-matrix))
- (let ((nxm-x (number-of-elements x))
- (nxm-y (number-of-elements y)))
- (declare (type fixnum nxm-x nxm-y))
-
- (if (not (= nxm-x nxm-y))
- (error "arguments X and Y to AXPY! not of the same size"))))
-
-(defmethod axpy! ((alpha number) (x complex-matrix) (y real-matrix))
- (error "cannot AXPY! a complex X to a real Y,
-don't know how to coerce COMPLEX to REAL"))
-
-;;
-(defmethod axpy! ((alpha cl:real) (x real-matrix) (y real-matrix))
- (axpy! (coerce alpha 'real-matrix-element-type) x y))
-
-(defmethod axpy! ((alpha #+(or cmu sbcl) double-float #-(or cmu sbcl) float) (x real-matrix) (y real-matrix))
- (let* ((nxm (number-of-elements y)))
- (declare (type fixnum nxm))
-
- (daxpy nxm alpha (store x) 1 (store y) 1)
- y))
-
-;;
-(defmethod axpy! ((alpha number) (x complex-matrix) (y complex-matrix))
- (let ((nxm (number-of-elements y))
- (c-alpha (complex-coerce alpha)))
- (declare (type complex-double-float c-alpha)
- (type fixnum nxm))
-
- (daxpy (* 2 nxm) c-alpha (store x) 1 (store y) 1)
- y))
-
-;;
-(defmethod axpy! ((alpha number) (x real-matrix) (y complex-matrix))
- (let* ((nxm (number-of-elements y))
- (store-x (store x))
- (store-y (store y))
- (c-alpha (complex-coerce alpha)))
- (declare (type complex-double-float c-alpha)
- (type fixnum nxm)
- (type (real-matrix-store-type *) store-x)
- (type (complex-matrix-store-type *) store-y))
-
- (daxpy nxm (realpart c-alpha) store-x 1 store-y 2)
- (with-vector-data-addresses ((addr-y store-y)
- (addr-x store-x))
- (incf-sap :double-float addr-y)
- (daxpy nxm (imagpart c-alpha) addr-x 1 addr-y 2))
- y))
+(defmethod axpy ((alpha number) (x real-matrix) (y complex-matrix))
+ (let ((result (copy y)))
+ (axpy! alpha x result)))
-;;
-(defmethod axpy! ((alpha number) (x complex-matrix) (y complex-matrix))
- (let ((nxm (number-of-elements y))
- (c-alpha (complex-coerce alpha)))
- (declare (type complex-double-float c-alpha)
- (type fixnum nxm))
-
- (zaxpy nxm c-alpha (store x) 1 (store y) 1)
- y))
\ No newline at end of file
+(defmethod axpy ((alpha number) (x complex-matrix) (y complex-matrix))
+ (let ((result (copy y)))
+ (axpy! alpha x result)))
\ No newline at end of file
diff --git a/src/blas.lisp b/src/blas.lisp
index 2875a5f..901b2fc 100644
--- a/src/blas.lisp
+++ b/src/blas.lisp
@@ -110,9 +110,9 @@
"
(n :integer :input)
(da :double-float :input)
- (dx (* :double-float))
+ (dx (* :double-float :inc head-x))
(incx :integer :input)
- (dy (* :double-float) :output)
+ (dy (* :double-float :inc head-y) :output)
(incy :integer :input)
)
@@ -286,11 +286,11 @@
"
(n :integer :input)
(za :complex-double-float)
- (zx (* :complex-double-float))
+ (zx (* :complex-double-float :inc head-x))
(incx :integer :input)
- (zy (* :complex-double-float) :output)
+ (zy (* :complex-double-float :inc head-y) :output)
(incy :integer :input)
- )
+)
(def-fortran-routine zcopy :void
"
@@ -632,9 +632,9 @@
Y(0),Y(2*INCY), ... , Y(2*(N-1)*INCY)
"
(n :integer :input)
- (dx (* :double-float) :input)
+ (dx (* :double-float :inc head-x) :input)
(incx :integer :input)
- (dy (* :double-float) :input)
+ (dy (* :double-float :inc head-y) :input)
(incy :integer :input)
)
@@ -744,12 +744,12 @@
(m :integer )
(n :integer )
(alpha :double-float )
- (a (* :double-float) )
+ (a (* :double-float :inc head-a) )
(lda :integer )
- (x (* :double-float) )
+ (x (* :double-float :inc head-x) )
(incx :integer )
(beta :double-float )
- (y (* :double-float) :output)
+ (y (* :double-float :inc head-y) :output)
(incy :integer )
)
diff --git a/src/complex-matrix.lisp b/src/complex-matrix.lisp
index 64382d7..0e02c9b 100644
--- a/src/complex-matrix.lisp
+++ b/src/complex-matrix.lisp
@@ -10,6 +10,9 @@
(deftype complex-matrix-store-type (size)
"The type of the storage structure for a COMPLEX-MATRIX"
`(simple-array double-float (,size)))
+
+(deftype complex-double-float ()
+ '(cl:complex (double-float * *)))
)
;;
@@ -298,7 +301,8 @@
(defun imagpart! (mat)
(cond
+ ((typep mat 'real-matrix) nil)
((typep mat 'complex-matrix) (make-instance 'real-matrix :store (store mat)
:nrows (nrows mat) :ncols (ncols mat)
:row-stride (* 2 (row-stride mat)) :col-stride (* 2 (col-stride mat))
- :head (+ 1 (* 2 (head mat)))))))
\ No newline at end of file
+ :head (+ 1 (* 2 (head mat)))))))
diff --git a/src/copy.lisp b/src/copy.lisp
index 68dff10..cf987d9 100644
--- a/src/copy.lisp
+++ b/src/copy.lisp
@@ -204,7 +204,7 @@ don't know how to coerce a COMPLEX to a REAL"))
(generate-typed-copy!-func complex-double-copy!-typed complex-matrix-store-type complex-matrix blas:zcopy)
(generate-typed-num-copy!-func complex-double-num-copy!-typed
- (complex (double-float * *)) complex-matrix-store-type complex-matrix
+ complex-double-float complex-matrix-store-type complex-matrix
blas:zcopy
(num
(1x1-z-array
diff --git a/src/gemm.lisp b/src/gemm.lisp
index c9e8736..c42b892 100644
--- a/src/gemm.lisp
+++ b/src/gemm.lisp
@@ -202,8 +202,8 @@
(defmethod gemm! ((alpha cl:real) (a real-matrix) (b real-matrix)
(beta cl:real) (c real-matrix)
&optional (job :nn))
- (real-double-gemm!-typed (coerce alpha 'real-matrix-element-type) a b
- (coerce beta 'real-matrix-element-type) c
+ (real-double-gemm!-typed (coerce alpha 'double-float) a b
+ (coerce beta 'double-float) c
job))
;;
@@ -215,6 +215,26 @@
(complex-double-gemm!-typed (complex-coerce alpha) a b
(complex-coerce beta) c job))
+
+(defmethod gemm! ((alpha cl:real) (a real-matrix) (b complex-matrix)
+ (beta number) (c complex-matrix)
+ &optional (job :nn))
+ (scal! beta c)
+ (real-double-gemm!-typed (coerce alpha 'double-float) a (realpart! b)
+ 1d0 (realpart! c) job)
+ (real-double-gemm!-typed (coerce alpha 'double-float) a (imagpart! b)
+ 1d0 (imagpart! c) job))
+
+(defmethod gemm! ((alpha complex) (a real-matrix) (b complex-matrix)
+ (beta number) (c complex-matrix)
+ &optional (job :nn))
+ (scal! beta c)
+ (real-double-gemm!-typed (coerce alpha 'double-float) a (realpart! b)
+ 1d0 (realpart! c) job)
+ (real-double-gemm!-typed (coerce alpha 'double-float) a (imagpart! b)
+ 1d0 (imagpart! c) job))
+
+
;;
(defmethod gemm! ((alpha number) (a standard-matrix) (b standard-matrix)
(beta number) (c complex-matrix)
@@ -306,4 +326,4 @@
(gemm! (complex-coerce alpha) a b
(complex-coerce beta) c
- job)))
+ job)))
\ No newline at end of file
diff --git a/src/gemv.lisp b/src/gemv.lisp
new file mode 100644
index 0000000..9adb6cc
--- /dev/null
+++ b/src/gemv.lisp
@@ -0,0 +1,35 @@
+(in-package :matlisp)
+
+(defun real-double-gemv!-typed (alpha A x beta y job)
+ (declare ((type double-float alpha beta)
+ (type real-matrix A x y)
+ (type symbol job)))
+ (symbol-macrolet
+ ((common-block
+ (mlet* (((st-a hd-a nr-a nc-a rs-a cs-a) (slot-values A '(store head number-of-rows number-of-cols row-stride col-stride))
+ :type ((real-matrix-store-type *) fixnum fixnum fixnum fixnum fixnum))
+ ((st-x hd-x rs-x) (slot-values x '(store head row-stride))
+ :type ((real-matrix-store-type *) fixnum fixnum))
+ ((st-y hd-y rs-y) (slot-values y '(store head row-stride))
+ :type ((real-matrix-store-type *) fixnum fixnum))
+ ((sym lda tf-job) (blas-matrix-compatible-p A) :type (nil fixnum (string 1))))
+ (if (not (string= tf-job "?"))
+ (blas:dgemv tf-job nr-a nc-a alpha st-a lda st-x rs-x beta st-y rs-y :head-a hd-a :head-x hd-x :head-y hd-y)
+ (dotimes (i nr-a)
+ (setf (matrix-ref-2d y i 0) (+ (* alpha (blas:ddot nc-a st-a cs-a st-x rs-x :head-x (+ hd-a (* i rs-a)) :head-y hd-x))
+ (* beta (matrix-ref-2d y i 0))))
+
+ (mlet* ((fortran-job (ecase job (:n "N") (:t "T")) :type ((string 1)))
+ ((st-a hd-a nr-a nc-a rs-a cs-a) (slot-values A '(store head number-of-rows number-of-cols row-stride col-stride))
+ :type ((real-matrix-store-type *) fixnum fixnum fixnum fixnum fixnum))
+ ((st-x hd-x rs-x) (slot-values x '(store head row-stride))
+ :type ((real-matrix-store-type *) fixnum fixnum))
+ ((st-y hd-y rs-y) (slot-values y '(store head row-stride))
+ :type ((real-matrix-store-type *) fixnum fixnum))
+ ((sym lda tf-job) (blas-matrix-compatible-p A fortran-job) :type (nil fixnum (string 1))))
+ (if (not (string= tf-job "?"))
+ (blas:dgemv tf-job nr-a nc-a alpha st-a lda st-x rs-x beta st-y rs-y :head-a hd-a :head-x hd-x :head-y hd-y)
+ (dotimes (i nr-a)
+ (setf (matrix-ref-2d y i 0) (+ (* alpha (blas:ddot nc-a st-a cs-a st-x rs-x :head-x (+ hd-a (* i rs-a)) :head-y hd-x))
+ (* beta (matrix-ref-2d y i 0))))
+
\ No newline at end of file
diff --git a/src/scal.lisp b/src/scal.lisp
index 8548a26..88da573 100644
--- a/src/scal.lisp
+++ b/src/scal.lisp
@@ -123,7 +123,7 @@ how to coerce COMPLEX to REAL"))
;;
(generate-typed-scal!-func complex-double-dscal!-typed double-float complex-matrix-store-type complex-matrix blas:zdscal)
-(generate-typed-scal!-func complex-double-zscal!-typed (complex (double-float * *)) complex-matrix-store-type complex-matrix blas:zscal)
+(generate-typed-scal!-func complex-double-zscal!-typed complex-double-float complex-matrix-store-type complex-matrix blas:zscal)
(defmethod scal! ((alpha cl:real) (x complex-matrix))
(complex-double-dscal!-typed (coerce alpha 'double-float) x))
diff --git a/src/standard-matrix.lisp b/src/standard-matrix.lisp
index 47c4ff9..3c2cfce 100644
--- a/src/standard-matrix.lisp
+++ b/src/standard-matrix.lisp
@@ -277,21 +277,6 @@ matrix and a number"))
(format stream "~%")))
;;
-(defun get-order-stride (matrix &optional (fortran-op "N"))
- (check-type matrix standard-matrix)
- (let ((rs (row-stride matrix))
- (cs (col-stride matrix)))
- (declare (type fixnum rs cs))
- (cond
- ((= cs 1) (values :row-major rs (cond
- ((string= fortran-op "N" ) "T")
- ((string= fortran-op "T" ) "N"))))
- ((= rs 1) (values :col-major cs (cond
- ((string= fortran-op "N" ) "N")
- ((string= fortran-op "T" ) "T"))))
- (t nil))))
-
-;;
(defun make-sub-matrix (matrix i j nrows ncols)
(declare (type standard-matrix matrix)
(type fixnum i j nrows ncols))
@@ -326,8 +311,20 @@ matrix and a number"))
(cs (col-stride matrix) :type fixnum)
(ne (number-of-elements matrix) :type fixnum))
(cond
- ((= nc 1) (values t rs ne))
- ((= nr 1) (values t cs ne))
- ((= rs (* nc cs)) (values t cs ne))
- ((= cs (* nr rs)) (values t rs ne))
- (t (values nil -1 -1)))))
\ No newline at end of file
+ ((or (= nc 1) (= cs (* nr rs))) (values t rs ne))
+ ((or (= nr 1) (= rs (* nc cs))) (values t cs ne))
+ (t (values nil -1 -1)))))
+
+;;
+(defun blas-matrix-compatible-p (matrix &optional (fortran-op "N"))
+ (declare (optimize (safety 0) (speed 3))
+ (type (or real-matrix complex-matrix) mat))
+ (mlet* (((rs cs) (slot-values mat '(row-stride col-stride))
+ :type (fixnum fixnum)))
+ (cond
+ ((= cs 1) (values :row-major rs (cond
+ ((string= fortran-op "N" ) "T")
+ ((string= fortran-op "T" ) "N"))))
+ ((= rs 1) (values :col-major cs fortran-op))
+ ;;Lets not confound lisp's type declaration.
+ (t (values nil -1 "?")))))
\ No newline at end of file
commit f4b94df81eab2264c071a5a6f592e6becf76f770
Merge: 2cbbcb6 8a2722a
Author: Akshay Srinivasan <aks...@gm...>
Date: Tue Mar 13 17:59:35 2012 +0530
Merge branch 'master' into matlisp-cffi
Conflicts:
src/blas.lisp
diff --cc src/blas.lisp
index 964f7b0,18c22d1..2875a5f
--- a/src/blas.lisp
+++ b/src/blas.lisp
@@@ -494,7 -495,24 +495,24 @@@
(incy :integer :input)
)
+ (def-fortran-routine mzdotu :void
+ (result (* :complex-double-float) :output)
+ (n :integer :input)
+ (zx (* :complex-double-float) :input)
+ (incx :integer :input)
+ (zy (* :complex-double-float) :input)
+ (incy :integer :input)
+ )
-(defun zdotu (n zx incx zy incy)
- (let ((result (make-array 2 :element-type 'double-float)))
++(let ((result (make-array 2 :element-type 'double-float)))
++ (defun zdotu (n zx incx zy incy)
+ (with-vector-data-addresses ((addr-result result)
+ (addr-zx zx)
+ (addr-zy zy))
+ (mzdotu addr-result n addr-zx incx addr-zy incy)
+ (complex (aref result 0) (aref result 1)))))
+
+ #-(and)
(def-fortran-routine zdotc :complex-double-float
"
Syntax
@@@ -538,6 -555,523 +556,24 @@@
(incy :integer :input)
)
+ (def-fortran-routine mzdotc :void
+ (result (* :complex-double-float) :output)
+ (n :integer :input)
+ (zx (* :complex-double-float) :input)
+ (incx :integer :input)
+ (zy (* :complex-double-float) :input)
+ (incy :integer :input)
+ )
+
-(defun zdotc (n zx incx zy incy)
- (let ((result (make-array 2 :element-type 'double-float)))
++(let ((result (make-array 2 :element-type 'double-float)))
++ (defun zdotc (n zx incx zy incy)
+ (with-vector-data-addresses ((addr-result result)
+ (addr-zx zx)
+ (addr-zy zy))
+ (mzdotc addr-result n addr-zx incx addr-zy incy)
+ (complex (aref result 0) (aref result 1)))))
+
+
-#|
-(declaim (inline fortran-daxpy))
-(def-alien-routine ("daxpy_" fortran-daxpy) void
- (n int :in-out)
- (da double-float :in-out)
- (dx (* double-float))
- (incx int :in-out)
- (dy (* double-float))
- (incy int :in-out)
-)
-(defun daxpy (
- n
- da
- dx
- incx
- dy
- incy
- )
- (declare
- (type fixnum n)
- (type double-float da)
- (type (simple-array double-float (*)) dx)
- (type fixnum incx)
- (type (simple-array double-float (*)) dy)
- (type fixnum incy)
- )
- (with-vector-data-addresses (
- (addr-dx dx)
- (addr-dy dy)
- )
- (multiple-value-bind ( return-value
- new-n
- new-incx
- new-incy
- )
- (fortran-daxpy
- n
- da
- addr-dx
- incx
- addr-dy
- incy
- )
- (declare (ignore return-value new-n new-incx new-incy))
- (values
- dy
- ))))
-
-(declaim (inline fortran-dcopy))
-(def-alien-routine ("dcopy_" fortran-dcopy) void
- (n int :in-out)
- (dx (* double-float))
- (incx int :in-out)
- (dy (* double-float))
- (incy int :in-out)
-)
-(defun dcopy (
- n
- dx
- incx
- dy
- incy
- )
- (declare
- (type fixnum n)
- (type (simple-array double-float (*)) dx)
- (type fixnum incx)
- (type (simple-array double-float (*)) dy)
- (type fixnum incy)
- )
- (with-vector-data-addresses (
- (addr-dx dx)
- (addr-dy dy)
- )
- (multiple-value-bind ( return-value
- new-n
- new-incx
- new-incy
- )
- (fortran-dcopy
- n
- addr-dx
- incx
- addr-dy
- incy
- )
- (declare (ignore return-value new-n new-incx new-incy))
- (values
- dy
- ))))
-
-
-(declaim (inline fortran-drot))
-(def-alien-routine ("drot_" fortran-drot) void
- (n int :in-out)
- (dx (* double-float))
- (incx int :in-out)
- (dy (* double-float))
- (incy int :in-out)
- (c double-float :in-out)
- (s double-float :in-out)
-)
-(defun drot (
- n
- dx
- incx
- dy
- incy
- c
- s
- )
- (declare
- (type fixnum n)
- (type (simple-array double-float (*)) dx)
- (type fixnum incx)
- (type (simple-array double-float (*)) dy)
- (type fixnum incy)
- (type double-float c)
- (type double-float s)
- )
- (with-vector-data-addresses (
- (addr-dx dx)
- (addr-dy dy)
- )
- (multiple-value-bind ( return-value
- new-n
- new-incx
- new-incy
- new-c
- new-s
- )
-
- (fortran-drot
- n
- addr-dx
- incx
- addr-dy
- incy
- c
- s
- )
- (declare (ignore return-value new-n new-incx new-incy))
- (values
- dx dy new-c new-s ))))
-
-(declaim (inline fortran-drotg))
-(def-alien-routine ("drotg_" fortran-drotg) void
- (da double-float :in-out)
- (db double-float :in-out)
- (c double-float :in-out)
- (s double-float :in-out)
-)
-(defun drotg (
- da
- db
- c
- s
- )
- (declare
- (type double-float da)
- (type double-float db)
- (type double-float c)
- (type double-float s)
- )
- (with-vector-data-addresses (
- )
- (multiple-value-bind ( return-value new-da new-db new-c new-s
- )
- (fortran-drotg
- da
- db
- c
- s
- )
- (declare (ignore return-value))
- (values
- new-da new-db new-c new-s ))))
-
-(declaim (inline fortran-dscal))
-(def-alien-routine ("dscal_" fortran-dscal) void
- (n int :in-out)
- (da double-float :in-out)
- (dx (* double-float))
- (incx int :in-out)
-)
-(defun dscal (
- n
- da
- dx
- incx
- )
- (declare
- (type fixnum n)
- (type double-float da)
- (type (simple-array double-float (*)) dx)
- (type fixnum incx)
- )
- (with-vector-data-addresses (
- (addr-dx dx)
- )
- (multiple-value-bind ( return-value
- new-n
- new-incx
- )
- (fortran-dscal
- n
- da
- addr-dx
- incx
- )
- (declare (ignore return-value new-n new-incx))
- (values
- dx ))))
-
-(declaim (inline fortran-dswap))
-(def-alien-routine ("dswap_" fortran-dswap) void
- (n int :in-out)
- (dx (* double-float))
- (incx int :in-out)
- (dy (* double-float))
- (incy int :in-out)
-)
-(defun dswap (
- n
- dx
- incx
- dy
- incy
- )
- (declare
- (type fixnum n)
- (type (simple-array double-float (*)) dx)
- (type fixnum incx)
- (type (simple-array double-float (*)) dy)
- (type fixnum incy)
- )
- (with-vector-data-addresses (
- (addr-dx dx)
- (addr-dy dy)
- )
- (multiple-value-bind ( return-value
- new-n
- new-incx
- new-incy
- )
- (fortran-dswap
- n
- addr-dx
- incx
- addr-dy
- incy
- )
- (declare (ignore return-value new-n new-incx new-incy))
- (values
- dx ))))
-
-
-(declaim (inline fortran-zaxpy))
-(def-alien-routine ("zaxpy_" fortran-zaxpy) void
- (n int :in-out)
- (za (* double-float))
- (zx (* double-float))
- (incx int :in-out)
- (zy (* double-float))
- (incy int :in-out)
-)
-(defun zaxpy (
- n
- za
- zx
- incx
- zy
- incy
- )
- (declare
- (type fixnum n)
- (type (simple-array double-float (*)) za)
- (type (simple-array double-float (*)) zx)
- (type fixnum incx)
- (type (simple-array double-float (*)) zy)
- (type fixnum incy)
- )
- (with-vector-data-addresses (
- (addr-za za)
- (addr-zx zx)
- (addr-zy zy)
- )
- (multiple-value-bind ( return-value
- new-n
- new-incx
- new-incy
- )
- (fortran-zaxpy
- n
- addr-za
- addr-zx
- incx
- addr-zy
- incy
- )
- (declare (ignore return-value new-n new-incx new-incy))
- (values
- zy ))))
-
-(declaim (inline fortran-zcopy))
-(def-alien-routine ("zcopy_" fortran-zcopy) void
- (n int :in-out)
- (zx (* double-float))
- (incx int :in-out)
- (zy (* double-float))
- (incy int :in-out)
-)
-(defun zcopy (
- n
- zx
- incx
- zy
- incy
- )
- (declare
- (type fixnum n)
- (type (simple-array double-float (*)) zx)
- (type fixnum incx)
- (type (simple-array double-float (*)) zy)
- (type fixnum incy)
- )
- (with-vector-data-addresses (
- (addr-zx zx)
- (addr-zy zy)
- )
- (multiple-value-bind ( return-value
- new-n
- new-incx
- new-incy
- )
- (fortran-zcopy
- n
- addr-zx
- incx
- addr-zy
- incy
- )
- (declare (ignore return-value new-n new-incx new-incy))
- (values
- zy ))))
-
-
-(declaim (inline fortran-zdscal))
-(def-alien-routine ("zdscal_" fortran-zdscal) void
- (n int :in-out)
- (da double-float :in-out)
- (zx (* double-float))
- (incx int :in-out)
-)
-(defun zdscal (
- n
- da
- zx
- incx
- )
- (declare
- (type fixnum n)
- (type double-float da)
- (type (simple-array double-float (*)) zx)
- (type fixnum incx)
- )
- (with-vector-data-addresses (
- (addr-zx zx)
- )
- (multiple-value-bind ( return-value
- new-n
- new-incx
- )
- (fortran-zdscal
- n
- da
- addr-zx
- incx
- )
- (declare (ignore return-value new-n new-incx))
- (values
- zx ))))
-
-(declaim (inline fortran-zrotg))
-(def-alien-routine ("zrotg_" fortran-zrotg) void
- (ca (* double-float))
- (cb (* double-float))
- (c double-float :in-out)
- (s (* double-float))
-)
-(defun zrotg (
- ca
- cb
- c
- s
- )
- (declare
- (type (simple-array double-float (*)) ca)
- (type (simple-array double-float (*)) cb)
- (type double-float c)
- (type (simple-array double-float (*)) s)
- )
- (with-vector-data-addresses (
- (addr-ca ca)
- (addr-cb cb)
- (addr-s s)
- )
- (multiple-value-bind ( return-value new-c
- )
- (fortran-zrotg
- addr-ca
- addr-cb
- c
- addr-s
- )
- (declare (ignore return-value))
- (values
- ca cb new-c s ))))
-
-(declaim (inline fortran-zscal))
-(def-alien-routine ("zscal_" fortran-zscal) void
- (n int :in-out)
- (za (* double-float))
- (zx (* double-float))
- (incx int :in-out)
-)
-(defun zscal (
- n
- za
- zx
- incx
- )
- (declare
- (type fixnum n)
- (type (simple-array double-float (*)) za)
- (type (simple-array double-float (*)) zx)
- (type fixnum incx)
- )
- (with-vector-data-addresses (
- (addr-za za)
- (addr-zx zx)
- )
- (multiple-value-bind ( return-value
- new-n
- new-incx
- )
- (fortran-zscal
- n
- addr-za
- addr-zx
- incx
- )
- (declare (ignore return-value new-n new-incx))
- (values
- zx ))))
-
-(declaim (inline fortran-zswap))
-(def-alien-routine ("zswap_" fortran-zswap) void
- (n int :in-out)
- (zx (* double-float))
- (incx int :in-out)
- (zy (* double-float))
- (incy int :in-out)
-)
-(defun zswap (
- n
- zx
- incx
- zy
- incy
- )
- (declare
- (type fixnum n)
- (type (simple-array double-float (*)) zx)
- (type fixnum incx)
- (type (simple-array double-float (*)) zy)
- (type fixnum incy)
- )
- (with-vector-data-addresses (
- (addr-zx zx)
- (addr-zy zy)
- )
- (multiple-value-bind ( return-value
- new-n
- new-incx
- new-incy
- )
- (fortran-zswap
- n
- addr-zx
- incx
- addr-zy
- incy
- )
- (declare (ignore return-value new-n new-incx new-incy))
- (values
- zx ))))
-
-|#
(def-fortran-routine idamax :integer
"
"
-----------------------------------------------------------------------
Summary of changes:
Makefile.am | 1 +
configure | 130 ++++++----------------
configure.ac | 254 +++++++++++++++++++++++++++---------------
lib-src/compat/Makefile.am | 12 ++
lib-src/compat/compat.f | 34 ++++++
lib/lazy-loader.lisp.in | 61 +++++++---
src/axpy.lisp | 269 ++++++++++++++------------------------------
src/blas.lisp | 56 ++++++++--
src/complex-matrix.lisp | 6 +-
src/copy.lisp | 57 +++++-----
src/gemm.lisp | 26 ++++-
src/gemv.lisp | 32 +++++
src/scal.lisp | 26 +++--
src/standard-matrix.lisp | 37 +++----
tests/blas.lisp | 8 +-
15 files changed, 541 insertions(+), 468 deletions(-)
create mode 100644 lib-src/compat/Makefile.am
create mode 100644 lib-src/compat/compat.f
create mode 100644 src/gemv.lisp
hooks/post-receive
--
matlisp
|
|
From: Raymond T. <toy...@gm...> - 2012-03-14 02:03:23
|
On 3/13/12 12:08 AM, Akshay Srinivasan wrote: >> hooks/post-receive > I tried doing something similar earlier, but my initial attempt threw > up some odd stuff. For some odd reason compiling the F77 wrapper with > -ff2c and linking this with ATLAS works, but throws up a segfault when > linked with the BLAS compiled by Matlisp (with -ff2c) (?). I'm a bit confused on what you were trying. I think that they need to be consistent in the calling convention. When building LAPACK ourselves, we're consistent. For ATLAS, we need to check the calling convention of zdotu and make sure we compile the wrapper with (or without) -ff2c. Of course, if ATLAS is using -ff2c, then we don't need the wrapper at all. Ray |
|
From: Akshay S. <aks...@gm...> - 2012-03-13 07:11:38
|
On 03/13/2012 10:24 AM, Raymond Toy wrote:
> This is an automated email from the git hooks/post-receive script.
> It was generated because a ref change was pushed to the repository
> containing the project "matlisp".
>
> The branch, master has been updated via
> 8a2722a7452cac1de3edd5fc8cc65d97c0d1a0a8 (commit) via
> ecab7fc4f75f615eabed198e0e7a325b2c055fd0 (commit) from
> 14707c10ced9ee08431313323cc1905ac77e2a21 (commit)
>
> Those revisions listed above that are new to this repository have
> not appeared on any other notification email; so we list those
> revisions in full, below.
>
> - Log
> -----------------------------------------------------------------
> commit 8a2722a7452cac1de3edd5fc8cc65d97c0d1a0a8 Author: Raymond
> Toy <toy...@gm...> Date: Mon Mar 12 21:54:26 2012 -0700
>
> Add wrapper for ZDOTC and ZDOTU to handle the differing possible
> Fortran return conventions. Tests pass. Still needs testing for
> ATLAS.
>
> Makefile.am: o Compile and install the new matlisp library.
>
> configure: o Regenerated
>
> configure.ac: o Remove checks if we need -ff2c. o Remove check for
> ATLAS compatible with f2c conventions.
>
> lib-src/compat/Makefile.am: o New file to compile and build the
> matlisp compat library.
>
> lib-src/compat/compat.f: o New file adding wrapper functions for
> ZDOTU and ZDOTC.
>
> lib/lazy-loader.lisp.in: o Define and use our libmatlisp
> compatibility library for both the standard LAPACK build and for
> ATLAS.
>
> src/blas.lisp: o Comment out old definitions for ZDOTU and ZDOTC.
> o Add definitions for our wrapper functions and add new functions
> ZDOTU and ZDOTC to call our wrappers.
>
> diff --git a/Makefile.am b/Makefile.am index 01c9e07..b073e06
> 100644 --- a/Makefile.am +++ b/Makefile.am @@ -15,6 +15,7 @@ all :
> (cd dfftpack; $(MAKE) install) (cd lib-src/toms715; $(MAKE)
> install) (cd lib-src/odepack; $(MAKE) install) + (cd
> lib-src/compat; $(MAKE) install) if !ATLAS (cd LAPACK/BLAS/SRC;
> $(MAKE) install) (cd LAPACK/SRC; $(MAKE) install) diff --git
> a/configure b/configure index 5465828..35ee372 100755 ---
> a/configure +++ b/configure @@ -15205,75 +15205,6 @@ else $as_echo
> "no" >&6; } fi
>
> -# Check to see if gfortran needs -ff2c. If needed, this option
> forces -# gfortran to be compatible with f2c. This is needed by
> matlisp's FFI -# so that it can correctly handle Fortran functions
> that return -# complex numbers. With f2c, complex numbers cause a
> hidden initial -# parameter that is a pointer to a structure where
> the complex result -# is stored. - -if test x"$F77" = x"gfortran";
> then - { $as_echo "$as_me:${as_lineno-$LINENO}: checking if
> gfortran needs -ff2c for functions returning complex numbers" >&5
> -$as_echo_n "checking if gfortran needs -ff2c for functions
> returning complex numbers... " >&6; } - # Just some silly Fortran
> function that returns a complex number. - cat > conftest.f <<EOF
> - complex*16 function f77_name(x) - complex*16 x -
> f77_name = 3*x - return - end -EOF - #echo $F77 $FFLAGS
> -c conftest.f - $F77 $FFLAGS -c conftest.f - # A C main routine
> that calls the fortran routine. This assumes the - # f2c calling
> convention of a hidden result structure. If the real - # part of
> the result matches our expectations, then main returns 0 - # (f2c
> compatible). Any other exit code means incompatible. - cat >
> ctest.c <<EOF -#include <stdio.h> -#include <math.h> - -struct c16
> -{ - double real; - double imag; -}; - - -extern void
> ${f77_name}(struct c16* result, struct c16* x); - -int main() -{ -
> int rc; - - struct c16 x, y; - - x.real = 1; - x.imag = 2; - -
> ${f77_name}(&y, &x); - - if (fabs(y.real - 3) < 1e-5) { - rc =
> 0; - } else { - rc = 1; - } - - return rc; -} -EOF - #echo
> $CC $CFLAGS -o a.out ctest.c conftest.o $LIBS - $CC $CFLAGS -o
> a.out ctest.c conftest.o $LIBS - NEED_FF2C=yes - F2C="-ff2c" - if
> a.out; then - NEED_FF2c=no - F2C="" - fi - rm -f ctest.c
> a.out - { $as_echo "$as_me:${as_lineno-$LINENO}: result:
> $NEED_FF2C" >&5 -$as_echo "$NEED_FF2C" >&6; } -fi
>
> # Allow user to use ATLAS if available. # We assume the standard
> names for the ATLAS libraries. @@ -15305,80 +15236,6 @@ else fi
>
>
> -# Check to see if the ATLAS libraries are compatible with
> matlisp's -# ffi. Basically the same test as above that checks to
> see if -ff2c -# is needed. We call zdotu which is a Fortran
> function returning a -# complex number. Matlisp assumes such
> functions return the result by -# storing the answer at address
> given by a hidden first parameter to -# the function. - -if test
> x"$atlas" = xtrue; then - { $as_echo
> "$as_me:${as_lineno-$LINENO}: checking if ATLAS is compatible with
> matlisp's FFI" >&5 -$as_echo_n "checking if ATLAS is compatible
> with matlisp's FFI... " >&6; } - # From the value of f77_name,
> figure out the actual name for - # Fortran's zdotu. - case
> $f77_name in - F77*) case $f77_name in - F77_NAME)
> ZDOTU="ZDOTU" ;; - F77_NAME_) ZDOTU="ZDOTU_" ;; -
> F77_NAME__) ZDOTU="ZDOTU_" ;; - esac - ;; - f77*) case
> $f77_name in - f77_name) ZDOTU="zdotu" ;; - f77_name_)
> ZDOTU="zdotu_" ;; - f77_name__) ZDOTU="zdotu_" ;; - esac - ;;
> - esac - - cat > conftest.c <<EOF - -#include <stdio.h> -#include
> <math.h> - -extern void ${ZDOTU}(double *, int *, double *, int *,
> double *, int *); - -int main() -{ - int rc; - - int n = 2; -
> int incx = 1; - double x[4]; - double out[2]; - x[0] = 1; -
> x[1] = 0; - x[2] = 2; - x[3] = 0; - out[0] = 0; - out[1] = 0; -
> - ${ZDOTU}(out, &n, x, &incx, x, &incx); - - if (fabs(out[0] - 5)
> < 1e-10) { - rc = 0; - } else { - printf("Actual output = %lg,
> %lg, instead of 5, 0\n", out[0], out[1]); - rc = 1; - } - -
> return rc; -} - -EOF - $CC $CFLAGS -c conftest.c - $F77 $FFLAGS
> -o a.out conftest.o -L${ATLAS_DIR} -latlas -lcblas -lf77blas
> -llapack - if a.out; then - { $as_echo
> "$as_me:${as_lineno-$LINENO}: result: yes" >&5 -$as_echo "yes" >&6;
> } - else - { $as_echo "$as_me:${as_lineno-$LINENO}: result: no"
> >&5 -$as_echo "no" >&6; } - as_fn_error $? "ATLAS libraries are
> not compatible with matlisp." "$LINENO" 5 - fi -fi -
>
>
>
> @@ -15407,7 +15264,7 @@ case $host in *) share_ext=so ;; esac
>
> -ac_config_files="$ac_config_files matlisp.mk Makefile start.lisp
> config.lisp lib/lazy-loader.lisp LAPACK/SRC/Makefile
> LAPACK/BLAS/SRC/Makefile dfftpack/Makefile
> lib-src/toms715/Makefile lib-src/odepack/Makefile
> src/f77-mangling.lisp" +ac_config_files="$ac_config_files
> matlisp.mk Makefile start.lisp config.lisp lib/lazy-loader.lisp
> LAPACK/SRC/Makefile LAPACK/BLAS/SRC/Makefile dfftpack/Makefile
> lib-src/toms715/Makefile lib-src/odepack/Makefile
> lib-src/compat/Makefile src/f77-mangling.lisp"
>
>
> echo FLIBS = $FLIBS @@ -16541,6 +16398,7 @@ do
> "dfftpack/Makefile") CONFIG_FILES="$CONFIG_FILES dfftpack/Makefile"
> ;; "lib-src/toms715/Makefile") CONFIG_FILES="$CONFIG_FILES
> lib-src/toms715/Makefile" ;; "lib-src/odepack/Makefile")
> CONFIG_FILES="$CONFIG_FILES lib-src/odepack/Makefile" ;; +
> "lib-src/compat/Makefile") CONFIG_FILES="$CONFIG_FILES
> lib-src/compat/Makefile" ;; "src/f77-mangling.lisp")
> CONFIG_FILES="$CONFIG_FILES src/f77-mangling.lisp" ;;
>
> *) as_fn_error $? "invalid argument: \`$ac_config_target'"
> "$LINENO" 5;; diff --git a/configure.ac b/configure.ac index
> fad1307..ef543fb 100644 --- a/configure.ac +++ b/configure.ac @@
> -227,74 +227,74 @@ else AC_MSG_RESULT([no]) fi
>
> -# Check to see if gfortran needs -ff2c. If needed, this option
> forces -# gfortran to be compatible with f2c. This is needed by
> matlisp's FFI -# so that it can correctly handle Fortran functions
> that return -# complex numbers. With f2c, complex numbers cause a
> hidden initial -# parameter that is a pointer to a structure where
> the complex result -# is stored. - -if test x"$F77" = x"gfortran";
> then - AC_MSG_CHECKING([if gfortran needs -ff2c for functions
> returning complex numbers]) - # Just some silly Fortran function
> that returns a complex number. - cat > conftest.f <<EOF -
> complex*16 function f77_name(x) - complex*16 x -
> f77_name = 3*x - return - end -EOF - #echo $F77 $FFLAGS
> -c conftest.f - $F77 $FFLAGS -c conftest.f - # A C main routine
> that calls the fortran routine. This assumes the - # f2c calling
> convention of a hidden result structure. If the real - # part of
> the result matches our expectations, then main returns 0 - # (f2c
> compatible). Any other exit code means incompatible. - cat >
> ctest.c <<EOF -#include <stdio.h> -#include <math.h> - -struct c16
> -{ - double real; - double imag; -}; - - -extern void
> ${f77_name}(struct c16* result, struct c16* x); - -int main() -{ -
> int rc; - - struct c16 x, y; - - x.real = 1; - x.imag = 2; - -
> ${f77_name}(&y, &x); - - if (fabs(y.real - 3) < 1e-5) { - rc =
> 0; - } else { - rc = 1; - } - - return rc; -} -EOF - #echo
> $CC $CFLAGS -o a.out ctest.c conftest.o $LIBS - $CC $CFLAGS -o
> a.out ctest.c conftest.o $LIBS - NEED_FF2C=yes - F2C="-ff2c" - if
> a.out; then - NEED_FF2c=no - F2C="" - fi - rm -f ctest.c
> a.out - AC_MSG_RESULT([$NEED_FF2C]) -fi - +dnl # Check to see if
> gfortran needs -ff2c. If needed, this option forces +dnl #
> gfortran to be compatible with f2c. This is needed by matlisp's
> FFI +dnl # so that it can correctly handle Fortran functions that
> return +dnl # complex numbers. With f2c, complex numbers cause a
> hidden initial +dnl # parameter that is a pointer to a structure
> where the complex result +dnl # is stored. + +dnl if test x"$F77"
> = x"gfortran"; then +dnl AC_MSG_CHECKING([if gfortran needs
> -ff2c for functions returning complex numbers]) +dnl # Just some
> silly Fortran function that returns a complex number. +dnl cat >
> conftest.f <<EOF +dnl complex*16 function f77_name(x) +dnl
> complex*16 x +dnl f77_name = 3*x +dnl return +dnl end
> +dnl EOF +dnl #echo $F77 $FFLAGS -c conftest.f +dnl $F77
> $FFLAGS -c conftest.f +dnl # A C main routine that calls the
> fortran routine. This assumes the +dnl # f2c calling convention
> of a hidden result structure. If the real +dnl # part of the
> result matches our expectations, then main returns 0 +dnl # (f2c
> compatible). Any other exit code means incompatible. +dnl cat >
> ctest.c <<EOF +dnl #include <stdio.h> +dnl #include <math.h> +dnl
> +dnl struct c16 +dnl { +dnl double real; +dnl double imag;
> +dnl }; +dnl +dnl +dnl extern void ${f77_name}(struct c16* result,
> struct c16* x); +dnl +dnl int main() +dnl { +dnl int rc; +dnl
> +dnl struct c16 x, y; +dnl +dnl x.real = 1; +dnl x.imag = 2;
> +dnl +dnl ${f77_name}(&y, &x); +dnl +dnl if (fabs(y.real - 3)
> < 1e-5) { +dnl rc = 0; +dnl } else { +dnl rc = 1; +dnl
> } +dnl +dnl return rc; +dnl } +dnl EOF +dnl #echo $CC $CFLAGS
> -o a.out ctest.c conftest.o $LIBS +dnl $CC $CFLAGS -o a.out
> ctest.c conftest.o $LIBS +dnl NEED_FF2C=yes +dnl F2C="-ff2c"
> +dnl if a.out; then +dnl NEED_FF2c=no +dnl F2C="" +dnl
> fi +dnl rm -f ctest.c a.out +dnl AC_MSG_RESULT([$NEED_FF2C]) +dnl
> fi +dnl # Allow user to use ATLAS if available. # We assume the
> standard names for the ATLAS libraries. AC_ARG_WITH([atlas], @@
> -314,77 +314,77 @@ AC_HELP_STRING([--with-atlas=libpath], [Location
> of the ATLAS libraries]), ]) AM_CONDITIONAL([ATLAS], [test x$atlas
> = xtrue])
>
> -# Check to see if the ATLAS libraries are compatible with
> matlisp's -# ffi. Basically the same test as above that checks to
> see if -ff2c -# is needed. We call zdotu which is a Fortran
> function returning a -# complex number. Matlisp assumes such
> functions return the result by -# storing the answer at address
> given by a hidden first parameter to -# the function. - -if test
> x"$atlas" = xtrue; then - AC_MSG_CHECKING([if ATLAS is compatible
> with matlisp's FFI]) - # From the value of f77_name, figure out
> the actual name for - # Fortran's zdotu. - case $f77_name in -
> F77*) case $f77_name in - F77_NAME) ZDOTU="ZDOTU" ;; -
> F77_NAME_) ZDOTU="ZDOTU_" ;; - F77_NAME__) ZDOTU="ZDOTU_" ;; -
> esac - ;; - f77*) case $f77_name in - f77_name)
> ZDOTU="zdotu" ;; - f77_name_) ZDOTU="zdotu_" ;; - f77_name__)
> ZDOTU="zdotu_" ;; - esac - ;; - esac - - cat > conftest.c
> <<EOF -[ -#include <stdio.h> -#include <math.h> - -extern void
> ${ZDOTU}(double *, int *, double *, int *, double *, int *); - -int
> main() -{ - int rc; - - int n = 2; - int incx = 1; - double
> x[4]; - double out[2]; - x[0] = 1; - x[1] = 0; - x[2] = 2; -
> x[3] = 0; - out[0] = 0; - out[1] = 0; - - ${ZDOTU}(out, &n, x,
> &incx, x, &incx); - - if (fabs(out[0] - 5) < 1e-10) { - rc = 0;
> - } else { - printf("Actual output = %lg, %lg, instead of 5,
> 0\n", out[0], out[1]); - rc = 1; - } - - return rc; -} -] -EOF
> - $CC $CFLAGS -c conftest.c - $F77 $FFLAGS -o a.out conftest.o
> -L${ATLAS_DIR} -latlas -lcblas -lf77blas -llapack - if a.out; then
> - AC_MSG_RESULT([yes]) - else - AC_MSG_RESULT([no]) -
> AC_MSG_ERROR([ATLAS libraries are not compatible with matlisp.]) -
> fi -fi - +dnl # Check to see if the ATLAS libraries are compatible
> with matlisp's +dnl # ffi. Basically the same test as above that
> checks to see if -ff2c +dnl # is needed. We call zdotu which is a
> Fortran function returning a +dnl # complex number. Matlisp
> assumes such functions return the result by +dnl # storing the
> answer at address given by a hidden first parameter to +dnl # the
> function. +dnl +dnl if test x"$atlas" = xtrue; then +dnl
> AC_MSG_CHECKING([if ATLAS is compatible with matlisp's FFI]) +dnl
> # From the value of f77_name, figure out the actual name for +dnl
> # Fortran's zdotu. +dnl case $f77_name in +dnl F77*) case
> $f77_name in +dnl F77_NAME) ZDOTU="ZDOTU" ;; +dnl
> F77_NAME_) ZDOTU="ZDOTU_" ;; +dnl F77_NAME__) ZDOTU="ZDOTU_"
> ;; +dnl esac +dnl ;; +dnl f77*) case $f77_name in +dnl
> f77_name) ZDOTU="zdotu" ;; +dnl f77_name_) ZDOTU="zdotu_" ;; +dnl
> f77_name__) ZDOTU="zdotu_" ;; +dnl esac +dnl ;; +dnl esac
> +dnl +dnl cat > conftest.c <<EOF +dnl [ +dnl #include <stdio.h>
> +dnl #include <math.h> +dnl +dnl extern void ${ZDOTU}(double *, int
> *, double *, int *, double *, int *); +dnl +dnl int main() +dnl {
> +dnl int rc; +dnl +dnl int n = 2; +dnl int incx = 1; +dnl
> double x[4]; +dnl double out[2]; +dnl x[0] = 1; +dnl x[1] =
> 0; +dnl x[2] = 2; +dnl x[3] = 0; +dnl out[0] = 0; +dnl out[1]
> = 0; +dnl +dnl ${ZDOTU}(out, &n, x, &incx, x, &incx); +dnl +dnl
> if (fabs(out[0] - 5) < 1e-10) { +dnl rc = 0; +dnl } else {
> +dnl printf("Actual output = %lg, %lg, instead of 5, 0\n", out[0],
> out[1]); +dnl rc = 1; +dnl } +dnl +dnl return rc; +dnl }
> +dnl ] +dnl EOF +dnl $CC $CFLAGS -c conftest.c +dnl $F77
> $FFLAGS -o a.out conftest.o -L${ATLAS_DIR} -latlas -lcblas
> -lf77blas -llapack +dnl if a.out; then +dnl AC_MSG_RESULT([yes])
> +dnl else +dnl AC_MSG_RESULT([no]) +dnl AC_MSG_ERROR([ATLAS
> libraries are not compatible with matlisp.]) +dnl fi +dnl fi +dnl
> dnl The following variables will be substituted into the .in files
> AC_SUBST(LISPEXEC) AC_SUBST(BLAS_OBJS) @@ -404,7 +404,7 @@
> AC_SUBST(FLIBS_OPTS) AC_SUBST(share_ext) AC_SUBST(GIT_VERSION)
> AC_SUBST(HAVE_QL) -AC_SUBST(F2C) +dnl AC_SUBST(F2C)
>
> echo host = $host
>
> @@ -425,6 +425,7 @@ AC_CONFIG_FILES([ dfftpack/Makefile
> lib-src/toms715/Makefile lib-src/odepack/Makefile +
> lib-src/compat/Makefile src/f77-mangling.lisp ])
>
> diff --git a/lib-src/compat/Makefile.am
> b/lib-src/compat/Makefile.am new file mode 100644 index
> 0000000..afa1180 --- /dev/null +++ b/lib-src/compat/Makefile.am @@
> -0,0 +1,12 @@ +lib_LTLIBRARIES = libmatlisp.la + +AM_FFLAGS =
> $(F2C) +if LIB32 +AM_FFLAGS += -m32 +endif +
> +libmatlisp_la_LDFLAGS = ../../LAPACK/BLAS/SRC/libblas.la
> $(FLIBS_OPTS) +libmatlisp_la_LIBADD = $(FLIBS) +
> +libmatlisp_la_SOURCES = \ +compat.f diff --git
> a/lib-src/compat/compat.f b/lib-src/compat/compat.f new file mode
> 100644 index 0000000..a4aabd5 --- /dev/null +++
> b/lib-src/compat/compat.f @@ -0,0 +1,34 @@ +c Compatibility
> routines for Matlisp. +c +c Fortran compilers differ on how a
> function returns a complex +c result. A compiler compatible with
> f2c will add a hidden first +c parameter where the result of the
> function will be stored. Other +c compilers will use the
> standard C structure return convention. +c +c Currently there
> are two BLAS functions the return complex numbers: +c ZDOTC and
> ZDOTU. + +c We don't want to deal with this difference in
> Matlisp, so we +c define the following wrapper functions to
> call the real BLAS +c routines using whatever the default
> convention is and then return +c the result in the first
> argument to routine. + subroutine mzdotc(result, n, x, incx, y,
> incy) + complex*16 result + complex*16 x(*), y(*) +
> integer n, incx, incy + external zdotc + complex*16 zdotc +
> result = zdotc(n, x, incx, y, incy) + return + end + +
> subroutine mzdotu(result, n, x, incx, y, incy) + complex*16
> result + complex*16 x(*), y(*) + integer n, incx, incy +
> external zdotu + complex*16 zdotu + result = zdotu(n, x,
> incx, y, incy) + return + end + diff --git
> a/lib/lazy-loader.lisp.in b/lib/lazy-loader.lisp.in index
> 7df2449..d558183 100644 --- a/lib/lazy-loader.lisp.in +++
> b/lib/lazy-loader.lisp.in @@ -137,6 +137,10 @@ (:darwin
> "libodepack.dylib") (t (:default "@libdir@/libodepack")))
>
> +(cffi:define-foreign-library matlisp + (:darwin
> "libmatlisp.dylib") + (t (:default "@libdir@/libmatlisp"))) + (if
> @ATLAS_P@ (progn ;; Define the ATLAS libraries and their locations.
> @@ -172,9 +176,11 @@ (cffi:use-foreign-library cblas)
> (cffi:use-foreign-library f77blas) (cffi:use-foreign-library
> lapack) + (cffi:use-foreign-library matlisp) ) (t
> (cffi:use-foreign-library blas) + (cffi:use-foreign-library
> matlisp) (cffi:use-foreign-library lapack))))
>
> #+:allegro diff --git a/src/blas.lisp b/src/blas.lisp index
> 827ab13..18c22d1 100644 --- a/src/blas.lisp +++ b/src/blas.lisp @@
> -451,6 +451,7 @@ (incy :integer :input) )
>
> +#-(and) (def-fortran-routine zdotu :complex-double-float " Syntax
> @@ -494,7 +495,24 @@ (incy :integer :input) )
>
> +(def-fortran-routine mzdotu :void + (result (*
> :complex-double-float) :output) + (n :integer :input) + (zx (*
> :complex-double-float) :input) + (incx :integer :input) + (zy (*
> :complex-double-float) :input) + (incy :integer :input) + )
>
> +(defun zdotu (n zx incx zy incy) + (let ((result (make-array 2
> :element-type 'double-float))) + (with-vector-data-addresses
> ((addr-result result) + (addr-zx zx) + (addr-zy zy)) +
> (mzdotu addr-result n addr-zx incx addr-zy incy) + (complex
> (aref result 0) (aref result 1))))) + +#-(and)
> (def-fortran-routine zdotc :complex-double-float " Syntax @@ -537,6
> +555,24 @@ (incy :integer :input) )
>
> +(def-fortran-routine mzdotc :void + (result (*
> :complex-double-float) :output) + (n :integer :input) + (zx (*
> :complex-double-float) :input) + (incx :integer :input) + (zy (*
> :complex-double-float) :input) + (incy :integer :input) + ) +
> +(defun zdotc (n zx incx zy incy) + (let ((result (make-array 2
> :element-type 'double-float))) + (with-vector-data-addresses
> ((addr-result result) + (addr-zx zx) + (addr-zy zy)) +
> (mzdotc addr-result n addr-zx incx addr-zy incy) + (complex
> (aref result 0) (aref result 1))))) + + #| (declaim (inline
> fortran-daxpy)) (def-alien-routine ("daxpy_" fortran-daxpy) void
>
> commit ecab7fc4f75f615eabed198e0e7a325b2c055fd0 Author: Raymond
> Toy <toy...@gm...> Date: Sun Mar 11 10:56:12 2012 -0700
>
> Regenerated.
>
> diff --git a/configure b/configure index 0f6e275..5465828 100755
> --- a/configure +++ b/configure @@ -615,8 +615,6 @@
> ac_subst_vars='am__EXEEXT_FALSE am__EXEEXT_TRUE LTLIBOBJS LIBOBJS
> -ATLAS_FALSE -ATLAS_TRUE F2C HAVE_QL GIT_VERSION @@ -634,6 +632,8
> @@ ATLAS_DIR NO_ATLAS_LAPACK_OBJS BLAS_OBJS LISPEXEC +ATLAS_FALSE
> +ATLAS_TRUE LIB32_FALSE LIB32_TRUE CPP @@ -4664,7 +4664,10 @@
> ac_link='$F77 -o conftest$ac_exeext $FFLAGS $LDFLAGS
> conftest.$ac_ext $LIBS >&5'
> ac_compiler_gnu=$ac_cv_f77_compiler_gnu
>
>
> -# Setup our environment variables based on the value of f77_name
> +# Setup our environment variables based on the value of f77_name.
> For +# the record, F77_EXTRA_UNDERSCORE means that any embedded
> underscores +# in the Fortran name are converted to two
> underscores for the C +# equivalent name. F77_LOWER_CASE=t
> F77_UNDERSCORE=t F77_EXTRA_UNDERSCORE=nil @@ -15155,7 +15158,8 @@
> FLIBS_OPTS=`echo $FLIBS | tr -s ' ' '\n' | $GREP -v gfortranbegin |
> $GREP -v x86
>
> GIT_VERSION=`git describe --dirty`
>
> -# Check that quicklisp exists. We need that currently to get
> cffi. +# Check that quicklisp exists. We need that currently to
> get cffi if +# it's not already available. as_ac_File=`$as_echo
> "ac_cv_file_$HOME/quicklisp/setup.lisp" | $as_tr_sh` { $as_echo
> "$as_me:${as_lineno-$LINENO}: checking for
> $HOME/quicklisp/setup.lisp" >&5 $as_echo_n "checking for
> $HOME/quicklisp/setup.lisp... " >&6; } @@ -15271,37 +15275,6 @@ EOF
> $as_echo "$NEED_FF2C" >&6; } fi
>
> - - - - - - - - - - - - - - - - - - - - -echo host = $host - -#
> Set the extension for shared libraries. This is not very robust.
> -case $host in - *darwin*) share_ext=dylib ;; - *) share_ext=so
> ;; -esac - -ac_config_files="$ac_config_files matlisp.mk Makefile
> start.lisp config.lisp lib/lazy-loader.lisp LAPACK/SRC/Makefile
> LAPACK/BLAS/SRC/Makefile dfftpack/Makefile
> lib-src/toms715/Makefile lib-src/odepack/Makefile
> src/f77-mangling.lisp" - - # Allow user to use ATLAS if available.
> # We assume the standard names for the ATLAS libraries.
>
> @@ -15334,7 +15307,10 @@ fi
>
> # Check to see if the ATLAS libraries are compatible with matlisp's
> # ffi. Basically the same test as above that checks to see if
> -ff2c -# is needed. +# is needed. We call zdotu which is a Fortran
> function returning a +# complex number. Matlisp assumes such
> functions return the result by +# storing the answer at address
> given by a hidden first parameter to +# the function.
>
> if test x"$atlas" = xtrue; then { $as_echo
> "$as_me:${as_lineno-$LINENO}: checking if ATLAS is compatible with
> matlisp's FFI" >&5 @@ -15345,13 +15321,13 @@ $as_echo_n "checking
> if ATLAS is compatible with matlisp's FFI... " >&6; } F77*) case
> $f77_name in F77_NAME) ZDOTU="ZDOTU" ;; F77_NAME_) ZDOTU="ZDOTU_"
> ;; - F77_NAME__) ZDOTU="ZDOTU__" ;; + F77_NAME__)
> ZDOTU="ZDOTU_" ;; esac ;; f77*) case $f77_name in f77_name)
> ZDOTU="zdotu" ;; f77_name_) ZDOTU="zdotu_" ;; - f77_name__)
> ZDOTU="zdotu__" ;; + f77_name__) ZDOTU="zdotu_" ;; esac ;; esac
> @@ -15383,6 +15359,7 @@ int main() if (fabs(out[0] - 5) < 1e-10) {
> rc = 0; } else { + printf("Actual output = %lg, %lg, instead of
> 5, 0\n", out[0], out[1]); rc = 1; }
>
> @@ -15402,6 +15379,37 @@ $as_echo "no" >&6; } fi fi
>
> + + + + + + + + + + + + + + + + + + + + +echo host = $host + +#
> Set the extension for shared libraries. This is not very robust.
> +case $host in + *darwin*) share_ext=dylib ;; + *) share_ext=so
> ;; +esac + +ac_config_files="$ac_config_files matlisp.mk Makefile
> start.lisp config.lisp lib/lazy-loader.lisp LAPACK/SRC/Makefile
> LAPACK/BLAS/SRC/Makefile dfftpack/Makefile
> lib-src/toms715/Makefile lib-src/odepack/Makefile
> src/f77-mangling.lisp" + + echo FLIBS = $FLIBS echo HAVE_QL =
> $HAVE_QL cat >confcache <<\_ACEOF
>
> -----------------------------------------------------------------------
>
>
>
Summary of changes: Makefile.am | 1 + configure
> | 198 +++++-------------------------- configure.ac
> | 281 ++++++++++++++++++++++----------------------
> lib-src/compat/Makefile.am | 12 ++ lib-src/compat/compat.f |
> 34 ++++++ lib/lazy-loader.lisp.in | 6 + src/blas.lisp | 36
> ++++++ 7 files changed, 262 insertions(+), 306 deletions(-) create
> mode 100644 lib-src/compat/Makefile.am create mode 100644
> lib-src/compat/compat.f
>
>
> hooks/post-receive
I tried doing something similar earlier, but my initial attempt threw
up some odd stuff. For some odd reason compiling the F77 wrapper with
-ff2c and linking this with ATLAS works, but throws up a segfault when
linked with the BLAS compiled by Matlisp (with -ff2c) (?).
Akshay
|
|
From: Raymond T. <rt...@us...> - 2012-03-13 04:54:52
|
This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "matlisp".
The branch, master has been updated
via 8a2722a7452cac1de3edd5fc8cc65d97c0d1a0a8 (commit)
via ecab7fc4f75f615eabed198e0e7a325b2c055fd0 (commit)
from 14707c10ced9ee08431313323cc1905ac77e2a21 (commit)
Those revisions listed above that are new to this repository have
not appeared on any other notification email; so we list those
revisions in full, below.
- Log -----------------------------------------------------------------
commit 8a2722a7452cac1de3edd5fc8cc65d97c0d1a0a8
Author: Raymond Toy <toy...@gm...>
Date: Mon Mar 12 21:54:26 2012 -0700
Add wrapper for ZDOTC and ZDOTU to handle the differing possible
Fortran return conventions. Tests pass. Still needs testing for
ATLAS.
Makefile.am:
o Compile and install the new matlisp library.
configure:
o Regenerated
configure.ac:
o Remove checks if we need -ff2c.
o Remove check for ATLAS compatible with f2c conventions.
lib-src/compat/Makefile.am:
o New file to compile and build the matlisp compat library.
lib-src/compat/compat.f:
o New file adding wrapper functions for ZDOTU and ZDOTC.
lib/lazy-loader.lisp.in:
o Define and use our libmatlisp compatibility library for both the
standard LAPACK build and for ATLAS.
src/blas.lisp:
o Comment out old definitions for ZDOTU and ZDOTC.
o Add definitions for our wrapper functions and add new functions
ZDOTU and ZDOTC to call our wrappers.
diff --git a/Makefile.am b/Makefile.am
index 01c9e07..b073e06 100644
--- a/Makefile.am
+++ b/Makefile.am
@@ -15,6 +15,7 @@ all :
(cd dfftpack; $(MAKE) install)
(cd lib-src/toms715; $(MAKE) install)
(cd lib-src/odepack; $(MAKE) install)
+ (cd lib-src/compat; $(MAKE) install)
if !ATLAS
(cd LAPACK/BLAS/SRC; $(MAKE) install)
(cd LAPACK/SRC; $(MAKE) install)
diff --git a/configure b/configure
index 5465828..35ee372 100755
--- a/configure
+++ b/configure
@@ -15205,75 +15205,6 @@ else
$as_echo "no" >&6; }
fi
-# Check to see if gfortran needs -ff2c. If needed, this option forces
-# gfortran to be compatible with f2c. This is needed by matlisp's FFI
-# so that it can correctly handle Fortran functions that return
-# complex numbers. With f2c, complex numbers cause a hidden initial
-# parameter that is a pointer to a structure where the complex result
-# is stored.
-
-if test x"$F77" = x"gfortran"; then
- { $as_echo "$as_me:${as_lineno-$LINENO}: checking if gfortran needs -ff2c for functions returning complex numbers" >&5
-$as_echo_n "checking if gfortran needs -ff2c for functions returning complex numbers... " >&6; }
- # Just some silly Fortran function that returns a complex number.
- cat > conftest.f <<EOF
- complex*16 function f77_name(x)
- complex*16 x
- f77_name = 3*x
- return
- end
-EOF
- #echo $F77 $FFLAGS -c conftest.f
- $F77 $FFLAGS -c conftest.f
- # A C main routine that calls the fortran routine. This assumes the
- # f2c calling convention of a hidden result structure. If the real
- # part of the result matches our expectations, then main returns 0
- # (f2c compatible). Any other exit code means incompatible.
- cat > ctest.c <<EOF
-#include <stdio.h>
-#include <math.h>
-
-struct c16
-{
- double real;
- double imag;
-};
-
-
-extern void ${f77_name}(struct c16* result, struct c16* x);
-
-int main()
-{
- int rc;
-
- struct c16 x, y;
-
- x.real = 1;
- x.imag = 2;
-
- ${f77_name}(&y, &x);
-
- if (fabs(y.real - 3) < 1e-5) {
- rc = 0;
- } else {
- rc = 1;
- }
-
- return rc;
-}
-EOF
- #echo $CC $CFLAGS -o a.out ctest.c conftest.o $LIBS
- $CC $CFLAGS -o a.out ctest.c conftest.o $LIBS
- NEED_FF2C=yes
- F2C="-ff2c"
- if a.out; then
- NEED_FF2c=no
- F2C=""
- fi
- rm -f ctest.c a.out
- { $as_echo "$as_me:${as_lineno-$LINENO}: result: $NEED_FF2C" >&5
-$as_echo "$NEED_FF2C" >&6; }
-fi
# Allow user to use ATLAS if available.
# We assume the standard names for the ATLAS libraries.
@@ -15305,80 +15236,6 @@ else
fi
-# Check to see if the ATLAS libraries are compatible with matlisp's
-# ffi. Basically the same test as above that checks to see if -ff2c
-# is needed. We call zdotu which is a Fortran function returning a
-# complex number. Matlisp assumes such functions return the result by
-# storing the answer at address given by a hidden first parameter to
-# the function.
-
-if test x"$atlas" = xtrue; then
- { $as_echo "$as_me:${as_lineno-$LINENO}: checking if ATLAS is compatible with matlisp's FFI" >&5
-$as_echo_n "checking if ATLAS is compatible with matlisp's FFI... " >&6; }
- # From the value of f77_name, figure out the actual name for
- # Fortran's zdotu.
- case $f77_name in
- F77*) case $f77_name in
- F77_NAME) ZDOTU="ZDOTU" ;;
- F77_NAME_) ZDOTU="ZDOTU_" ;;
- F77_NAME__) ZDOTU="ZDOTU_" ;;
- esac
- ;;
- f77*) case $f77_name in
- f77_name) ZDOTU="zdotu" ;;
- f77_name_) ZDOTU="zdotu_" ;;
- f77_name__) ZDOTU="zdotu_" ;;
- esac
- ;;
- esac
-
- cat > conftest.c <<EOF
-
-#include <stdio.h>
-#include <math.h>
-
-extern void ${ZDOTU}(double *, int *, double *, int *, double *, int *);
-
-int main()
-{
- int rc;
-
- int n = 2;
- int incx = 1;
- double x[4];
- double out[2];
- x[0] = 1;
- x[1] = 0;
- x[2] = 2;
- x[3] = 0;
- out[0] = 0;
- out[1] = 0;
-
- ${ZDOTU}(out, &n, x, &incx, x, &incx);
-
- if (fabs(out[0] - 5) < 1e-10) {
- rc = 0;
- } else {
- printf("Actual output = %lg, %lg, instead of 5, 0\n", out[0], out[1]);
- rc = 1;
- }
-
- return rc;
-}
-
-EOF
- $CC $CFLAGS -c conftest.c
- $F77 $FFLAGS -o a.out conftest.o -L${ATLAS_DIR} -latlas -lcblas -lf77blas -llapack
- if a.out; then
- { $as_echo "$as_me:${as_lineno-$LINENO}: result: yes" >&5
-$as_echo "yes" >&6; }
- else
- { $as_echo "$as_me:${as_lineno-$LINENO}: result: no" >&5
-$as_echo "no" >&6; }
- as_fn_error $? "ATLAS libraries are not compatible with matlisp." "$LINENO" 5
- fi
-fi
-
@@ -15407,7 +15264,7 @@ case $host in
*) share_ext=so ;;
esac
-ac_config_files="$ac_config_files matlisp.mk Makefile start.lisp config.lisp lib/lazy-loader.lisp LAPACK/SRC/Makefile LAPACK/BLAS/SRC/Makefile dfftpack/Makefile lib-src/toms715/Makefile lib-src/odepack/Makefile src/f77-mangling.lisp"
+ac_config_files="$ac_config_files matlisp.mk Makefile start.lisp config.lisp lib/lazy-loader.lisp LAPACK/SRC/Makefile LAPACK/BLAS/SRC/Makefile dfftpack/Makefile lib-src/toms715/Makefile lib-src/odepack/Makefile lib-src/compat/Makefile src/f77-mangling.lisp"
echo FLIBS = $FLIBS
@@ -16541,6 +16398,7 @@ do
"dfftpack/Makefile") CONFIG_FILES="$CONFIG_FILES dfftpack/Makefile" ;;
"lib-src/toms715/Makefile") CONFIG_FILES="$CONFIG_FILES lib-src/toms715/Makefile" ;;
"lib-src/odepack/Makefile") CONFIG_FILES="$CONFIG_FILES lib-src/odepack/Makefile" ;;
+ "lib-src/compat/Makefile") CONFIG_FILES="$CONFIG_FILES lib-src/compat/Makefile" ;;
"src/f77-mangling.lisp") CONFIG_FILES="$CONFIG_FILES src/f77-mangling.lisp" ;;
*) as_fn_error $? "invalid argument: \`$ac_config_target'" "$LINENO" 5;;
diff --git a/configure.ac b/configure.ac
index fad1307..ef543fb 100644
--- a/configure.ac
+++ b/configure.ac
@@ -227,74 +227,74 @@ else
AC_MSG_RESULT([no])
fi
-# Check to see if gfortran needs -ff2c. If needed, this option forces
-# gfortran to be compatible with f2c. This is needed by matlisp's FFI
-# so that it can correctly handle Fortran functions that return
-# complex numbers. With f2c, complex numbers cause a hidden initial
-# parameter that is a pointer to a structure where the complex result
-# is stored.
-
-if test x"$F77" = x"gfortran"; then
- AC_MSG_CHECKING([if gfortran needs -ff2c for functions returning complex numbers])
- # Just some silly Fortran function that returns a complex number.
- cat > conftest.f <<EOF
- complex*16 function f77_name(x)
- complex*16 x
- f77_name = 3*x
- return
- end
-EOF
- #echo $F77 $FFLAGS -c conftest.f
- $F77 $FFLAGS -c conftest.f
- # A C main routine that calls the fortran routine. This assumes the
- # f2c calling convention of a hidden result structure. If the real
- # part of the result matches our expectations, then main returns 0
- # (f2c compatible). Any other exit code means incompatible.
- cat > ctest.c <<EOF
-#include <stdio.h>
-#include <math.h>
-
-struct c16
-{
- double real;
- double imag;
-};
-
-
-extern void ${f77_name}(struct c16* result, struct c16* x);
-
-int main()
-{
- int rc;
-
- struct c16 x, y;
-
- x.real = 1;
- x.imag = 2;
-
- ${f77_name}(&y, &x);
-
- if (fabs(y.real - 3) < 1e-5) {
- rc = 0;
- } else {
- rc = 1;
- }
-
- return rc;
-}
-EOF
- #echo $CC $CFLAGS -o a.out ctest.c conftest.o $LIBS
- $CC $CFLAGS -o a.out ctest.c conftest.o $LIBS
- NEED_FF2C=yes
- F2C="-ff2c"
- if a.out; then
- NEED_FF2c=no
- F2C=""
- fi
- rm -f ctest.c a.out
- AC_MSG_RESULT([$NEED_FF2C])
-fi
-
+dnl # Check to see if gfortran needs -ff2c. If needed, this option forces
+dnl # gfortran to be compatible with f2c. This is needed by matlisp's FFI
+dnl # so that it can correctly handle Fortran functions that return
+dnl # complex numbers. With f2c, complex numbers cause a hidden initial
+dnl # parameter that is a pointer to a structure where the complex result
+dnl # is stored.
+
+dnl if test x"$F77" = x"gfortran"; then
+dnl AC_MSG_CHECKING([if gfortran needs -ff2c for functions returning complex numbers])
+dnl # Just some silly Fortran function that returns a complex number.
+dnl cat > conftest.f <<EOF
+dnl complex*16 function f77_name(x)
+dnl complex*16 x
+dnl f77_name = 3*x
+dnl return
+dnl end
+dnl EOF
+dnl #echo $F77 $FFLAGS -c conftest.f
+dnl $F77 $FFLAGS -c conftest.f
+dnl # A C main routine that calls the fortran routine. This assumes the
+dnl # f2c calling convention of a hidden result structure. If the real
+dnl # part of the result matches our expectations, then main returns 0
+dnl # (f2c compatible). Any other exit code means incompatible.
+dnl cat > ctest.c <<EOF
+dnl #include <stdio.h>
+dnl #include <math.h>
+dnl
+dnl struct c16
+dnl {
+dnl double real;
+dnl double imag;
+dnl };
+dnl
+dnl
+dnl extern void ${f77_name}(struct c16* result, struct c16* x);
+dnl
+dnl int main()
+dnl {
+dnl int rc;
+dnl
+dnl struct c16 x, y;
+dnl
+dnl x.real = 1;
+dnl x.imag = 2;
+dnl
+dnl ${f77_name}(&y, &x);
+dnl
+dnl if (fabs(y.real - 3) < 1e-5) {
+dnl rc = 0;
+dnl } else {
+dnl rc = 1;
+dnl }
+dnl
+dnl return rc;
+dnl }
+dnl EOF
+dnl #echo $CC $CFLAGS -o a.out ctest.c conftest.o $LIBS
+dnl $CC $CFLAGS -o a.out ctest.c conftest.o $LIBS
+dnl NEED_FF2C=yes
+dnl F2C="-ff2c"
+dnl if a.out; then
+dnl NEED_FF2c=no
+dnl F2C=""
+dnl fi
+dnl rm -f ctest.c a.out
+dnl AC_MSG_RESULT([$NEED_FF2C])
+dnl fi
+dnl
# Allow user to use ATLAS if available.
# We assume the standard names for the ATLAS libraries.
AC_ARG_WITH([atlas],
@@ -314,77 +314,77 @@ AC_HELP_STRING([--with-atlas=libpath], [Location of the ATLAS libraries]),
])
AM_CONDITIONAL([ATLAS], [test x$atlas = xtrue])
-# Check to see if the ATLAS libraries are compatible with matlisp's
-# ffi. Basically the same test as above that checks to see if -ff2c
-# is needed. We call zdotu which is a Fortran function returning a
-# complex number. Matlisp assumes such functions return the result by
-# storing the answer at address given by a hidden first parameter to
-# the function.
-
-if test x"$atlas" = xtrue; then
- AC_MSG_CHECKING([if ATLAS is compatible with matlisp's FFI])
- # From the value of f77_name, figure out the actual name for
- # Fortran's zdotu.
- case $f77_name in
- F77*) case $f77_name in
- F77_NAME) ZDOTU="ZDOTU" ;;
- F77_NAME_) ZDOTU="ZDOTU_" ;;
- F77_NAME__) ZDOTU="ZDOTU_" ;;
- esac
- ;;
- f77*) case $f77_name in
- f77_name) ZDOTU="zdotu" ;;
- f77_name_) ZDOTU="zdotu_" ;;
- f77_name__) ZDOTU="zdotu_" ;;
- esac
- ;;
- esac
-
- cat > conftest.c <<EOF
-[
-#include <stdio.h>
-#include <math.h>
-
-extern void ${ZDOTU}(double *, int *, double *, int *, double *, int *);
-
-int main()
-{
- int rc;
-
- int n = 2;
- int incx = 1;
- double x[4];
- double out[2];
- x[0] = 1;
- x[1] = 0;
- x[2] = 2;
- x[3] = 0;
- out[0] = 0;
- out[1] = 0;
-
- ${ZDOTU}(out, &n, x, &incx, x, &incx);
-
- if (fabs(out[0] - 5) < 1e-10) {
- rc = 0;
- } else {
- printf("Actual output = %lg, %lg, instead of 5, 0\n", out[0], out[1]);
- rc = 1;
- }
-
- return rc;
-}
-]
-EOF
- $CC $CFLAGS -c conftest.c
- $F77 $FFLAGS -o a.out conftest.o -L${ATLAS_DIR} -latlas -lcblas -lf77blas -llapack
- if a.out; then
- AC_MSG_RESULT([yes])
- else
- AC_MSG_RESULT([no])
- AC_MSG_ERROR([ATLAS libraries are not compatible with matlisp.])
- fi
-fi
-
+dnl # Check to see if the ATLAS libraries are compatible with matlisp's
+dnl # ffi. Basically the same test as above that checks to see if -ff2c
+dnl # is needed. We call zdotu which is a Fortran function returning a
+dnl # complex number. Matlisp assumes such functions return the result by
+dnl # storing the answer at address given by a hidden first parameter to
+dnl # the function.
+dnl
+dnl if test x"$atlas" = xtrue; then
+dnl AC_MSG_CHECKING([if ATLAS is compatible with matlisp's FFI])
+dnl # From the value of f77_name, figure out the actual name for
+dnl # Fortran's zdotu.
+dnl case $f77_name in
+dnl F77*) case $f77_name in
+dnl F77_NAME) ZDOTU="ZDOTU" ;;
+dnl F77_NAME_) ZDOTU="ZDOTU_" ;;
+dnl F77_NAME__) ZDOTU="ZDOTU_" ;;
+dnl esac
+dnl ;;
+dnl f77*) case $f77_name in
+dnl f77_name) ZDOTU="zdotu" ;;
+dnl f77_name_) ZDOTU="zdotu_" ;;
+dnl f77_name__) ZDOTU="zdotu_" ;;
+dnl esac
+dnl ;;
+dnl esac
+dnl
+dnl cat > conftest.c <<EOF
+dnl [
+dnl #include <stdio.h>
+dnl #include <math.h>
+dnl
+dnl extern void ${ZDOTU}(double *, int *, double *, int *, double *, int *);
+dnl
+dnl int main()
+dnl {
+dnl int rc;
+dnl
+dnl int n = 2;
+dnl int incx = 1;
+dnl double x[4];
+dnl double out[2];
+dnl x[0] = 1;
+dnl x[1] = 0;
+dnl x[2] = 2;
+dnl x[3] = 0;
+dnl out[0] = 0;
+dnl out[1] = 0;
+dnl
+dnl ${ZDOTU}(out, &n, x, &incx, x, &incx);
+dnl
+dnl if (fabs(out[0] - 5) < 1e-10) {
+dnl rc = 0;
+dnl } else {
+dnl printf("Actual output = %lg, %lg, instead of 5, 0\n", out[0], out[1]);
+dnl rc = 1;
+dnl }
+dnl
+dnl return rc;
+dnl }
+dnl ]
+dnl EOF
+dnl $CC $CFLAGS -c conftest.c
+dnl $F77 $FFLAGS -o a.out conftest.o -L${ATLAS_DIR} -latlas -lcblas -lf77blas -llapack
+dnl if a.out; then
+dnl AC_MSG_RESULT([yes])
+dnl else
+dnl AC_MSG_RESULT([no])
+dnl AC_MSG_ERROR([ATLAS libraries are not compatible with matlisp.])
+dnl fi
+dnl fi
+dnl
dnl The following variables will be substituted into the .in files
AC_SUBST(LISPEXEC)
AC_SUBST(BLAS_OBJS)
@@ -404,7 +404,7 @@ AC_SUBST(FLIBS_OPTS)
AC_SUBST(share_ext)
AC_SUBST(GIT_VERSION)
AC_SUBST(HAVE_QL)
-AC_SUBST(F2C)
+dnl AC_SUBST(F2C)
echo host = $host
@@ -425,6 +425,7 @@ AC_CONFIG_FILES([
dfftpack/Makefile
lib-src/toms715/Makefile
lib-src/odepack/Makefile
+ lib-src/compat/Makefile
src/f77-mangling.lisp
])
diff --git a/lib-src/compat/Makefile.am b/lib-src/compat/Makefile.am
new file mode 100644
index 0000000..afa1180
--- /dev/null
+++ b/lib-src/compat/Makefile.am
@@ -0,0 +1,12 @@
+lib_LTLIBRARIES = libmatlisp.la
+
+AM_FFLAGS = $(F2C)
+if LIB32
+AM_FFLAGS += -m32
+endif
+
+libmatlisp_la_LDFLAGS = ../../LAPACK/BLAS/SRC/libblas.la $(FLIBS_OPTS)
+libmatlisp_la_LIBADD = $(FLIBS)
+
+libmatlisp_la_SOURCES = \
+compat.f
diff --git a/lib-src/compat/compat.f b/lib-src/compat/compat.f
new file mode 100644
index 0000000..a4aabd5
--- /dev/null
+++ b/lib-src/compat/compat.f
@@ -0,0 +1,34 @@
+c Compatibility routines for Matlisp.
+c
+c Fortran compilers differ on how a function returns a complex
+c result. A compiler compatible with f2c will add a hidden first
+c parameter where the result of the function will be stored. Other
+c compilers will use the standard C structure return convention.
+c
+c Currently there are two BLAS functions the return complex numbers:
+c ZDOTC and ZDOTU.
+
+c We don't want to deal with this difference in Matlisp, so we
+c define the following wrapper functions to call the real BLAS
+c routines using whatever the default convention is and then return
+c the result in the first argument to routine.
+ subroutine mzdotc(result, n, x, incx, y, incy)
+ complex*16 result
+ complex*16 x(*), y(*)
+ integer n, incx, incy
+ external zdotc
+ complex*16 zdotc
+ result = zdotc(n, x, incx, y, incy)
+ return
+ end
+
+ subroutine mzdotu(result, n, x, incx, y, incy)
+ complex*16 result
+ complex*16 x(*), y(*)
+ integer n, incx, incy
+ external zdotu
+ complex*16 zdotu
+ result = zdotu(n, x, incx, y, incy)
+ return
+ end
+
diff --git a/lib/lazy-loader.lisp.in b/lib/lazy-loader.lisp.in
index 7df2449..d558183 100644
--- a/lib/lazy-loader.lisp.in
+++ b/lib/lazy-loader.lisp.in
@@ -137,6 +137,10 @@
(:darwin "libodepack.dylib")
(t (:default "@libdir@/libodepack")))
+(cffi:define-foreign-library matlisp
+ (:darwin "libmatlisp.dylib")
+ (t (:default "@libdir@/libmatlisp")))
+
(if @ATLAS_P@
(progn
;; Define the ATLAS libraries and their locations.
@@ -172,9 +176,11 @@
(cffi:use-foreign-library cblas)
(cffi:use-foreign-library f77blas)
(cffi:use-foreign-library lapack)
+ (cffi:use-foreign-library matlisp)
)
(t
(cffi:use-foreign-library blas)
+ (cffi:use-foreign-library matlisp)
(cffi:use-foreign-library lapack))))
#+:allegro
diff --git a/src/blas.lisp b/src/blas.lisp
index 827ab13..18c22d1 100644
--- a/src/blas.lisp
+++ b/src/blas.lisp
@@ -451,6 +451,7 @@
(incy :integer :input)
)
+#-(and)
(def-fortran-routine zdotu :complex-double-float
"
Syntax
@@ -494,7 +495,24 @@
(incy :integer :input)
)
+(def-fortran-routine mzdotu :void
+ (result (* :complex-double-float) :output)
+ (n :integer :input)
+ (zx (* :complex-double-float) :input)
+ (incx :integer :input)
+ (zy (* :complex-double-float) :input)
+ (incy :integer :input)
+ )
+(defun zdotu (n zx incx zy incy)
+ (let ((result (make-array 2 :element-type 'double-float)))
+ (with-vector-data-addresses ((addr-result result)
+ (addr-zx zx)
+ (addr-zy zy))
+ (mzdotu addr-result n addr-zx incx addr-zy incy)
+ (complex (aref result 0) (aref result 1)))))
+
+#-(and)
(def-fortran-routine zdotc :complex-double-float
"
Syntax
@@ -537,6 +555,24 @@
(incy :integer :input)
)
+(def-fortran-routine mzdotc :void
+ (result (* :complex-double-float) :output)
+ (n :integer :input)
+ (zx (* :complex-double-float) :input)
+ (incx :integer :input)
+ (zy (* :complex-double-float) :input)
+ (incy :integer :input)
+ )
+
+(defun zdotc (n zx incx zy incy)
+ (let ((result (make-array 2 :element-type 'double-float)))
+ (with-vector-data-addresses ((addr-result result)
+ (addr-zx zx)
+ (addr-zy zy))
+ (mzdotc addr-result n addr-zx incx addr-zy incy)
+ (complex (aref result 0) (aref result 1)))))
+
+
#|
(declaim (inline fortran-daxpy))
(def-alien-routine ("daxpy_" fortran-daxpy) void
commit ecab7fc4f75f615eabed198e0e7a325b2c055fd0
Author: Raymond Toy <toy...@gm...>
Date: Sun Mar 11 10:56:12 2012 -0700
Regenerated.
diff --git a/configure b/configure
index 0f6e275..5465828 100755
--- a/configure
+++ b/configure
@@ -615,8 +615,6 @@ ac_subst_vars='am__EXEEXT_FALSE
am__EXEEXT_TRUE
LTLIBOBJS
LIBOBJS
-ATLAS_FALSE
-ATLAS_TRUE
F2C
HAVE_QL
GIT_VERSION
@@ -634,6 +632,8 @@ ATLAS_DIR
NO_ATLAS_LAPACK_OBJS
BLAS_OBJS
LISPEXEC
+ATLAS_FALSE
+ATLAS_TRUE
LIB32_FALSE
LIB32_TRUE
CPP
@@ -4664,7 +4664,10 @@ ac_link='$F77 -o conftest$ac_exeext $FFLAGS $LDFLAGS conftest.$ac_ext $LIBS >&5'
ac_compiler_gnu=$ac_cv_f77_compiler_gnu
-# Setup our environment variables based on the value of f77_name
+# Setup our environment variables based on the value of f77_name. For
+# the record, F77_EXTRA_UNDERSCORE means that any embedded underscores
+# in the Fortran name are converted to two underscores for the C
+# equivalent name.
F77_LOWER_CASE=t
F77_UNDERSCORE=t
F77_EXTRA_UNDERSCORE=nil
@@ -15155,7 +15158,8 @@ FLIBS_OPTS=`echo $FLIBS | tr -s ' ' '\n' | $GREP -v gfortranbegin | $GREP -v x86
GIT_VERSION=`git describe --dirty`
-# Check that quicklisp exists. We need that currently to get cffi.
+# Check that quicklisp exists. We need that currently to get cffi if
+# it's not already available.
as_ac_File=`$as_echo "ac_cv_file_$HOME/quicklisp/setup.lisp" | $as_tr_sh`
{ $as_echo "$as_me:${as_lineno-$LINENO}: checking for $HOME/quicklisp/setup.lisp" >&5
$as_echo_n "checking for $HOME/quicklisp/setup.lisp... " >&6; }
@@ -15271,37 +15275,6 @@ EOF
$as_echo "$NEED_FF2C" >&6; }
fi
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-echo host = $host
-
-# Set the extension for shared libraries. This is not very robust.
-case $host in
- *darwin*) share_ext=dylib ;;
- *) share_ext=so ;;
-esac
-
-ac_config_files="$ac_config_files matlisp.mk Makefile start.lisp config.lisp lib/lazy-loader.lisp LAPACK/SRC/Makefile LAPACK/BLAS/SRC/Makefile dfftpack/Makefile lib-src/toms715/Makefile lib-src/odepack/Makefile src/f77-mangling.lisp"
-
-
# Allow user to use ATLAS if available.
# We assume the standard names for the ATLAS libraries.
@@ -15334,7 +15307,10 @@ fi
# Check to see if the ATLAS libraries are compatible with matlisp's
# ffi. Basically the same test as above that checks to see if -ff2c
-# is needed.
+# is needed. We call zdotu which is a Fortran function returning a
+# complex number. Matlisp assumes such functions return the result by
+# storing the answer at address given by a hidden first parameter to
+# the function.
if test x"$atlas" = xtrue; then
{ $as_echo "$as_me:${as_lineno-$LINENO}: checking if ATLAS is compatible with matlisp's FFI" >&5
@@ -15345,13 +15321,13 @@ $as_echo_n "checking if ATLAS is compatible with matlisp's FFI... " >&6; }
F77*) case $f77_name in
F77_NAME) ZDOTU="ZDOTU" ;;
F77_NAME_) ZDOTU="ZDOTU_" ;;
- F77_NAME__) ZDOTU="ZDOTU__" ;;
+ F77_NAME__) ZDOTU="ZDOTU_" ;;
esac
;;
f77*) case $f77_name in
f77_name) ZDOTU="zdotu" ;;
f77_name_) ZDOTU="zdotu_" ;;
- f77_name__) ZDOTU="zdotu__" ;;
+ f77_name__) ZDOTU="zdotu_" ;;
esac
;;
esac
@@ -15383,6 +15359,7 @@ int main()
if (fabs(out[0] - 5) < 1e-10) {
rc = 0;
} else {
+ printf("Actual output = %lg, %lg, instead of 5, 0\n", out[0], out[1]);
rc = 1;
}
@@ -15402,6 +15379,37 @@ $as_echo "no" >&6; }
fi
fi
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+echo host = $host
+
+# Set the extension for shared libraries. This is not very robust.
+case $host in
+ *darwin*) share_ext=dylib ;;
+ *) share_ext=so ;;
+esac
+
+ac_config_files="$ac_config_files matlisp.mk Makefile start.lisp config.lisp lib/lazy-loader.lisp LAPACK/SRC/Makefile LAPACK/BLAS/SRC/Makefile dfftpack/Makefile lib-src/toms715/Makefile lib-src/odepack/Makefile src/f77-mangling.lisp"
+
+
echo FLIBS = $FLIBS
echo HAVE_QL = $HAVE_QL
cat >confcache <<\_ACEOF
-----------------------------------------------------------------------
Summary of changes:
Makefile.am | 1 +
configure | 198 +++++--------------------------
configure.ac | 281 ++++++++++++++++++++++----------------------
lib-src/compat/Makefile.am | 12 ++
lib-src/compat/compat.f | 34 ++++++
lib/lazy-loader.lisp.in | 6 +
src/blas.lisp | 36 ++++++
7 files changed, 262 insertions(+), 306 deletions(-)
create mode 100644 lib-src/compat/Makefile.am
create mode 100644 lib-src/compat/compat.f
hooks/post-receive
--
matlisp
|
|
From: Akshay S. <ak...@us...> - 2012-03-12 17:28:26
|
This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "matlisp".
The branch, matlisp-cffi has been updated
via 2cbbcb64c997e7ed95c2e9a344347ca9606fce66 (commit)
via b0f1cb2dd42338c9189c83cbcbcb177eaf1c7845 (commit)
via 98cdecd68f57b4a561ac8f68a0ede4a0374a6a95 (commit)
via 097c9251aaa15f702cc8b9701832dbdc1d9bcb13 (commit)
via b0415f8f7e3a4af5682f22aed1885e20f3484188 (commit)
from 96f2abfd9395d6f25520c3b828c4c69a0f35d8a3 (commit)
Those revisions listed above that are new to this repository have
not appeared on any other notification email; so we list those
revisions in full, below.
- Log -----------------------------------------------------------------
commit 2cbbcb64c997e7ed95c2e9a344347ca9606fce66
Author: Akshay Srinivasan <aks...@gm...>
Date: Mon Mar 12 22:48:55 2012 +0530
-> Support for copy(!) and scal(!) are reasonably complete.
-> Restricted the store type of {real/complex}-matrix to one-dimension.
-> Usual churning in utilities.lisp, made mlet* even more incomprehensible (but nicer!).
diff --git a/packages.lisp b/packages.lisp
index 383af89..566dbd6 100644
--- a/packages.lisp
+++ b/packages.lisp
@@ -298,6 +298,7 @@
"*PRINT-MATRIX*"
"AXPY!"
"AXPY"
+ "BLAS-COPYABLE-P"
"COL-VECTOR-P"
"COMPLEX-COERCE"
"COMPLEX-MATRIX"
@@ -420,6 +421,7 @@
"SWAP!"
"TR"
"TRANSPOSE"
+ "TRANSPOSE!"
"VEC"
"UNLOAD-BLAS-&-LAPACK-LIBRARIES"
"ZEROS"
diff --git a/src/axpy.lisp b/src/axpy.lisp
index c3a6baf..57104cf 100644
--- a/src/axpy.lisp
+++ b/src/axpy.lisp
@@ -150,8 +150,8 @@
(store-y (store y))
(store-result (store result)))
(declare (type fixnum n m nxm)
- (type (real-matrix-store-type (*)) store-y)
- (type (complex-matrix-store-type (*)) store-x store-result))
+ (type (real-matrix-store-type *) store-y)
+ (type (complex-matrix-store-type *) store-x store-result))
(zcopy nxm store-x 1 store-result 1) ;; same as (COPY! x result)
(zdscal nxm alpha store-result 1) ;; same as (SCAL! alpha result)
@@ -192,8 +192,8 @@
(c-alpha (complex-coerce alpha)))
(declare (type complex-double-float c-alpha)
(type fixnum n m nxm)
- (type (real-matrix-store-type (*)) store-x)
- (type (complex-matrix-store-type (*)) store-y store-result))
+ (type (real-matrix-store-type *) store-x)
+ (type (complex-matrix-store-type *) store-y store-result))
(dcopy nxm store-x 1 store-result 2)
(zscal nxm c-alpha store-result 1)
@@ -209,7 +209,7 @@
(c-alpha (complex-coerce alpha)))
(declare (type complex-double-float c-alpha)
(type fixnum nxm)
- (type (complex-matrix-store-type (*)) store-result))
+ (type (complex-matrix-store-type *) store-result))
(zscal nxm c-alpha store-result 1)
(daxpy nxm 1.0d0 (store y) 1 store-result 2)
@@ -283,8 +283,8 @@ don't know how to coerce COMPLEX to REAL"))
(c-alpha (complex-coerce alpha)))
(declare (type complex-double-float c-alpha)
(type fixnum nxm)
- (type (real-matrix-store-type (*)) store-x)
- (type (complex-matrix-store-type (*)) store-y))
+ (type (real-matrix-store-type *) store-x)
+ (type (complex-matrix-store-type *) store-y))
(daxpy nxm (realpart c-alpha) store-x 1 store-y 2)
(with-vector-data-addresses ((addr-y store-y)
diff --git a/src/blas.lisp b/src/blas.lisp
index 999b389..964f7b0 100644
--- a/src/blas.lisp
+++ b/src/blas.lisp
@@ -1,4 +1,4 @@
-t;;; -*- Mode: lisp; Syntax: ansi-common-lisp; Package: :blas; Base: 10 -*-
+;;; -*- Mode: lisp; Syntax: ansi-common-lisp; Package: :blas; Base: 10 -*-
;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;
@@ -204,7 +204,7 @@ t;;; -*- Mode: lisp; Syntax: ansi-common-lisp; Package: :blas; Base: 10 -*-
"
(n :integer :input)
(da :double-float :input)
- (dx (* :double-float :inc head-dx) :output)
+ (dx (* :double-float :inc head-x) :output)
(incx :integer :input)
)
@@ -363,7 +363,7 @@ t;;; -*- Mode: lisp; Syntax: ansi-common-lisp; Package: :blas; Base: 10 -*-
"
(n :integer :input)
(da :double-float :input)
- (zx (* :complex-double-float :inc head-zx) :output)
+ (zx (* :complex-double-float :inc head-x) :output)
(incx :integer :input)
)
@@ -405,7 +405,7 @@ t;;; -*- Mode: lisp; Syntax: ansi-common-lisp; Package: :blas; Base: 10 -*-
"
(n :integer :input)
(za :complex-double-float)
- (zx (* :complex-double-float) :output)
+ (zx (* :complex-double-float :inc head-x) :output)
(incx :integer :input)
)
@@ -530,7 +530,8 @@ t;;; -*- Mode: lisp; Syntax: ansi-common-lisp; Package: :blas; Base: 10 -*-
considered in the operations are:
Y(0),Y(2*INCY), ... , Y(2*(N-1)*INCY)
-" (n :integer :input)
+"
+ (n :integer :input)
(zx (* :complex-double-float) :input)
(incx :integer :input)
(zy (* :complex-double-float) :input)
@@ -547,6 +548,11 @@ t;;; -*- Mode: lisp; Syntax: ansi-common-lisp; Package: :blas; Base: 10 -*-
(def-fortran-routine dasum :double-float
"
+ Purpose
+ =======
+
+ Takes the sum of the absolute values.
+
"
(n :integer :input)
(dx (* :double-float) :input)
diff --git a/src/compat.lisp b/src/compat.lisp
index 4b0f329..b30383e 100644
--- a/src/compat.lisp
+++ b/src/compat.lisp
@@ -93,11 +93,11 @@
(deftype float-matrix-array-type (size)
"Defines the same type as (REAL-MATRIX-STORE-TYPE (*))"
- `(simple-array double-float ,size))
+ `(simple-array double-float (,size)))
(deftype complex-matrix-array-type (size)
"Defines the same type as (COMPLEX-MATRIX-STORE-TYPE (*))"
- `(simple-array double-float ,size))
+ `(simple-array double-float (,size)))
(deftype float-matrix ()
"Defines the same type as REAL-MATRIX"
diff --git a/src/complex-matrix.lisp b/src/complex-matrix.lisp
index 48e7d0f..64382d7 100644
--- a/src/complex-matrix.lisp
+++ b/src/complex-matrix.lisp
@@ -9,7 +9,7 @@
(deftype complex-matrix-store-type (size)
"The type of the storage structure for a COMPLEX-MATRIX"
- `(simple-array double-float ,size))
+ `(simple-array double-float (,size)))
)
;;
@@ -37,7 +37,7 @@
(defclass complex-matrix (standard-matrix)
((store
:initform nil
- :type (simple-array complex-matrix-element-type (*))))
+ :type (complex-matrix-store-type *)))
(:documentation "A class of matrices with complex elements."))
;;
@@ -55,13 +55,13 @@
;;
(defmethod matrix-ref-1d ((matrix complex-matrix) (idx fixnum))
(let ((store (store matrix)))
- (declare (type (complex-matrix-store-type (*)) store))
+ (declare (type (complex-matrix-store-type *) store))
(complex (aref store (* 2 idx)) (aref store (+ 1 (* 2 idx))))))
(defmethod (setf matrix-ref-1d) ((value number) (matrix complex-matrix) (idx fixnum))
(let ((store (store matrix))
(coerced-value (complex-coerce value)))
- (declare (type (complex-matrix-store-type (*)) store))
+ (declare (type (complex-matrix-store-type *) store))
(setf (aref store (* 2 idx)) (realpart coerced-value)
(aref store (+ 1 (* 2 idx))) (imagpart coerced-value))))
@@ -121,7 +121,7 @@
(size (* n m))
(store (allocate-complex-store size)))
(declare (type fixnum n m size)
- (type (complex-matrix-store-type (*)) store))
+ (type (complex-matrix-store-type *) store))
(multiple-value-bind (row-stride col-stride)
(ecase order
(:row-major (values m 1))
@@ -151,7 +151,7 @@
(size (* n m))
(store (allocate-complex-store size)))
(declare (type fixnum n m size)
- (type (complex-matrix-store-type (*)) store))
+ (type (complex-matrix-store-type *) store))
(multiple-value-bind (row-stride col-stride)
(ecase order
(:row-major (values m 1))
@@ -182,7 +182,7 @@
(let* ((n (length seq))
(store (allocate-complex-store n)))
(declare (type fixnum n)
- (type (complex-matrix-store-type (*)) store))
+ (type (complex-matrix-store-type *) store))
(dotimes (k n)
(declare (type fixnum k))
(let* ((val (complex-coerce (elt seq k)))
diff --git a/src/copy.lisp b/src/copy.lisp
index 586aaf0..68dff10 100644
--- a/src/copy.lisp
+++ b/src/copy.lisp
@@ -79,29 +79,6 @@
(in-package "MATLISP")
;;
-(defun blas-copyable-p (matrix)
- (declare (optimize (safety 0) (speed 3))
- (type (or real-matrix complex-matrix) matrix))
- (mlet* ((nr (nrows matrix) :type fixnum)
- (nc (ncols matrix) :type fixnum)
- (rs (row-stride matrix) :type fixnum)
- (cs (col-stride matrix) :type fixnum)
- (ne (number-of-elements matrix) :type fixnum))
- (cond
- ((= nc 1) (values t rs ne))
- ((= nr 1) (values t cs ne))
- ((= rs (* nc cs)) (values t cs ne))
- ((= cs (* nr rs)) (values t rs ne))
- (t (values nil -1 -1)))))
-
-;;
-(defmacro with-transpose! (matlst &rest body)
- `(progn
- ,@(mapcar #'(lambda (mat) `(transpose! ,mat)) matlst)
- ,@body
- ,@(mapcar #'(lambda (mat) `(transpose! ,mat)) matlst)))
-
-;;
(defmacro generate-typed-copy!-func (func store-type matrix-type blas-func)
`(defun ,func (mat-a mat-b)
(declare (type ,matrix-type mat-a mat-b)
@@ -111,7 +88,7 @@
((hd-a st-a sz) (slot-values mat-a '(head store number-of-elements)) :type (fixnum (,store-type *) fixnum))
((hd-b st-b) (slot-values mat-b '(head store)) :type (fixnum (,store-type *))))
(if (and cp-a cp-b)
- (,blas-func sz (store mat-a) inc-a (store mat-b) inc-b :head-x hd-a :head-y hd-b)
+ (,blas-func sz st-a inc-a st-b inc-b :head-x hd-a :head-y hd-b)
(symbol-macrolet
((common-code
(mlet* (((nr-a nc-a rs-a cs-a) (slot-values mat-a '(number-of-rows number-of-cols row-stride col-stride))
@@ -128,8 +105,33 @@
mat-b)))
;;
-(defvar *1x1-real-array* (make-array 1 :element-type 'double-float))
-(defvar *1x1-complex-array* (make-array 2 :element-type 'double-float))
+(defmacro generate-typed-num-copy!-func (func element-type store-type matrix-type blas-func
+ array-decl)
+ (let ((num-arg (car array-decl)))
+ (destructuring-bind (var var-maker-form var-setf-form &key type) (cadr array-decl)
+ `(mlet* (((,var) (,@var-maker-form) ,@(if type
+ `(:type (,type)))))
+ (defun ,func (,num-arg mat-x)
+ (declare (type ,element-type ,num-arg)
+ (type ,matrix-type mat-x)
+ (optimize (safety 0) (speed 3)))
+ (,@var-setf-form)
+ (mlet* (((cp-x inc-x sz-x) (blas-copyable-p mat-x) :type (boolean fixnum nil))
+ ((hd-x st-x sz) (slot-values mat-x '(head store number-of-elements)) :type (fixnum (,store-type *) fixnum)))
+ (if cp-x
+ (,blas-func sz ,var 0 st-x inc-x :head-y hd-x)
+ (symbol-macrolet
+ ((common-code
+ (mlet* (((nr-x nc-x rs-x cs-x) (slot-values mat-x '(number-of-rows number-of-cols row-stride col-stride))
+ :type (fixnum fixnum fixnum fixnum)))
+ (loop for i from 0 below nr-x
+ do (,blas-func nc-x ,var 0 st-x cs-x :head-y (+ hd-x (* i rs-x)))))))
+ ;;Choose the smaller of the loops
+ (if (> (nrows mat-x) (ncols mat-x))
+ (with-transpose! (mat-x)
+ common-code)
+ common-code)))
+ mat-x))))))
;;
(defgeneric copy! (matrix new-matrix)
@@ -157,102 +159,72 @@
"))
(defmethod copy! :before ((x standard-matrix) (y standard-matrix))
- (let ((nxm-x (number-of-elements x))
- (nxm-y (number-of-elements y)))
- (declare (type fixnum nxm-x nxm-y))
- (if (not (= nxm-x nxm-y))
- (warn "arguments X,Y to COPY! are of different sizes"))))
+ (mlet* (((nr-x nc-x) (slot-values x '(number-of-rows number-of-cols)) :type (fixnum fixnum))
+ ((nr-y nc-y) (slot-values y '(number-of-rows number-of-cols)) :type (fixnum fixnum)))
+ (unless (and (= nr-x nr-y) (= nc-x nc-y))
+ (error "Arguments X,Y to COPY! are of different dimensions."))))
;;
-(generate-typed-copy!-func real-double-copy!-typed real-matrix-store-type real-matrix blas:dcopy)
-
-(defmethod copy! ((x real-matrix) (y real-matrix))
- (let* ((nxm-x (number-of-elements x))
- (nxm-y (number-of-elements y))
- (nxm (min nxm-x nxm-y)))
- (declare (type fixnum nxm-x nxm-y nxm))
- (dcopy nxm (store x) 1 (store y) 1)
+(defmethod copy! ((x standard-matrix) (y standard-matrix))
+ (mlet* (((nr-x nc-x) (slot-values x '(number-of-rows number-of-cols))
+ :type (fixnum fixnum)))
+ (dotimes (i nr-x)
+ (dotimes (j nc-x)
+ (declare (type fixnum i j))
+ (setf (matrix-ref-2d y i j) (matrix-ref-2d x i j))))
y))
-(defmethod copy! ((x real-matrix) (y complex-matrix))
- (let* ((nxm-x (number-of-elements x))
- (nxm-y (number-of-elements y))
- (nxm (min nxm-x nxm-y)))
- (declare (type fixnum nxm-x nxm-y nxm))
+;;
+(generate-typed-copy!-func real-double-copy!-typed real-matrix-store-type real-matrix blas:dcopy)
- ;; Set the imaginary parts of Y to zero.
- (zdscal nxm 0.0d0 (store y) 1)
+(generate-typed-num-copy!-func real-double-num-copy!-typed
+ double-float real-matrix-store-type real-matrix
+ blas:dcopy
+ (num
+ (1x1-array
+ (allocate-real-store 1)
+ (setf (aref 1x1-array 0) num)
+ :type (real-matrix-store-type 1))))
- ;; Copy the elements of X to the real parts of Y.
- (dcopy nxm (store x) 1 (store y) 2)
- y))
-
(defmethod copy! ((x complex-matrix) (y real-matrix))
- (error "cannot copy a COMPLEX-MATRIX into a REAL-MATRIX,
+ (error "Cannot copy a COMPLEX-MATRIX into a REAL-MATRIX,
don't know how to coerce a COMPLEX to a REAL"))
+(defmethod copy! ((x complex) (y real-matrix))
+ (error "Cannot copy ~a to ~a, don't know how to coerce COMPLEX to REAL"
+ x y))
-;;
-(generate-typed-copy!-func complex-double-copy!-typed complex-matrix-store-type complex-matrix blas:zcopy)
-
-(defmethod copy! ((x complex-matrix) (y complex-matrix))
- (let* ((nxm-x (number-of-elements x))
- (nxm-y (number-of-elements y))
- (nxm (min nxm-x nxm-y)))
- (declare (type fixnum nxm-x nxm-y nxm))
- (dcopy (* 2 nxm) (store x) 1 (store y) 1)
- y))
-
-(defmethod copy! ((x standard-matrix) (y standard-matrix))
- (let* ((nxm-x (number-of-elements x))
- (nxm-y (number-of-elements y))
- (nxm (min nxm-x nxm-y)))
- (declare (type fixnum nxm-x nxm-y nxm))
-
- (dotimes (i nxm)
- (declare (type fixnum i))
- (setf (matrix-ref y i) (matrix-ref x i)))
-
- y))
-
-(defmethod copy! ((x #+(or cmu sbcl) double-float #-(or cmu sbcl) float) (y real-matrix))
- (let ((nxm (number-of-elements y)))
- (setf (aref *1x1-real-array* 0) x)
- (dcopy nxm *1x1-real-array* 0 (store y) 1)
- y))
+(defmethod copy! ((x real-matrix) (y real-matrix))
+ (real-double-copy!-typed x y))
(defmethod copy! ((x cl:real) (y real-matrix))
- (let ((nxm (number-of-elements y)))
- (setf (aref *1x1-real-array* 0) (coerce x 'real-matrix-element-type))
- (dcopy nxm *1x1-real-array* 0 (store y) 1)
- y))
+ (real-double-num-copy!-typed (coerce x 'double-float) y))
-(defmethod copy! ((x complex) (y real-matrix))
- (error "cannot copy ~a to ~a, don't know how to coerce COMPLEX to REAL"
- x
- y))
+;;
+(generate-typed-copy!-func complex-double-copy!-typed complex-matrix-store-type complex-matrix blas:zcopy)
-(defmethod copy! ((x #+:cmu kernel::complex-double-float
- #+:sbcl sb-kernel::complex-double-float
- #-(or cmu sbcl) complex) (y complex-matrix))
- (let ((nxm (number-of-elements y)))
+(generate-typed-num-copy!-func complex-double-num-copy!-typed
+ (complex (double-float * *)) complex-matrix-store-type complex-matrix
+ blas:zcopy
+ (num
+ (1x1-z-array
+ (allocate-complex-store 1)
+ (setf (aref 1x1-z-array 0) (realpart num)
+ (aref 1x1-z-array 1) (imagpart num))
+ :type (complex-matrix-store-type 2))))
- #-(or cmu sbcl) (setq x (complex-coerce x))
+(defmethod copy! ((x complex-matrix) (y complex-matrix))
+ (complex-double-copy!-typed x y))
- (setf (aref *1x1-complex-array* 0) (realpart x))
- (setf (aref *1x1-complex-array* 1) (imagpart x))
- (zcopy nxm *1x1-complex-array* 0 (store y) 1)
- y))
+(defmethod copy! ((x real-matrix) (y complex-matrix))
+ (real-double-copy!-typed x (realpart! y))
+ (scal! 0d0 (imagpart! y))
+ y)
(defmethod copy! ((x number) (y complex-matrix))
- (let ((nxm (number-of-elements y)))
- (setq x (complex-coerce x))
- (setf (aref *1x1-complex-array* 0) (realpart x))
- (setf (aref *1x1-complex-array* 1) (imagpart x))
- (zcopy nxm *1x1-complex-array* 0 (store y) 1)
- y))
+ (complex-double-num-copy!-typed (complex-coerce x) y))
-;;
+;;;;
(defgeneric copy (matrix)
(:documentation
"
@@ -264,32 +236,23 @@ don't know how to coerce a COMPLEX to a REAL"))
=======
Return a copy of the matrix X"))
-(defmethod copy ((matrix standard-matrix))
- (make-instance 'standard-matrix :nrows (nrows matrix) :ncols (ncols matrix) :store (copy-seq (store matrix))))
-
(defmethod copy ((matrix real-matrix))
- (let* ((size (number-of-elements matrix))
- (n (nrows matrix))
- (m (ncols matrix))
+ (let* ((n (nrows matrix))
+ (m (nrows matrix))
(result (make-real-matrix-dim n m)))
- (declare (type fixnum size n m))
- (blas:dcopy size (store matrix) 1 (store result) 1)
- result))
+ (declare (type fixnum n m))
+ (copy! matrix result)))
(defmethod copy ((matrix complex-matrix))
- (let* ((size (number-of-elements matrix))
- (n (nrows matrix))
+ (let* ((n (nrows matrix))
(m (ncols matrix))
(result (make-complex-matrix-dim n m)))
- (declare (type fixnum size n m))
- (blas:zcopy size (store matrix) 1 (store result) 1)
- result))
+ (declare (type fixnum n m))
+ (copy! matrix result)))
(defmethod copy ((matrix number))
matrix)
-
-
;;
(defgeneric convert-to-lisp-array (matrix)
(:documentation
diff --git a/src/dot.lisp b/src/dot.lisp
index 4b643c2..badff89 100644
--- a/src/dot.lisp
+++ b/src/dot.lisp
@@ -133,8 +133,8 @@
(store-x (store x))
(store-y (store y)))
(declare (type fixnum nxm)
- (type (real-matrix-store-type (*)) store-x)
- (type (complex-matrix-store-type (*)) store-y))
+ (type (real-matrix-store-type *) store-x)
+ (type (complex-matrix-store-type *) store-y))
(let ((realpart (ddot nxm store-x 1 store-y 2))
(imagpart (with-vector-data-addresses ((addr-x store-x)
@@ -186,8 +186,8 @@
(store-x (store x))
(store-y (store y)))
(declare (type fixnum nxm)
- (type (real-matrix-store-type (*)) store-y)
- (type (complex-matrix-store-type (*)) store-x))
+ (type (real-matrix-store-type *) store-y)
+ (type (complex-matrix-store-type *) store-x))
(let ((realpart (ddot nxm store-x 2 store-y 1))
(imagpart (with-vector-data-addresses ((addr-x store-x)
diff --git a/src/map.lisp b/src/map.lisp
index 8206e55..5b2afa4 100644
--- a/src/map.lisp
+++ b/src/map.lisp
@@ -124,7 +124,7 @@
(let ((nxm (number-of-elements mat))
(store (store mat)))
(declare (type fixnum nxm)
- (type (real-matrix-store-type (*)) store))
+ (type (real-matrix-store-type *) store))
(dotimes (k nxm mat)
(declare (type fixnum k))
(setf (aref store k) (funcall func (aref store k))))))
@@ -143,7 +143,7 @@
(let ((nxm (number-of-elements mat))
(store (store mat)))
(declare (type fixnum nxm)
- (type (real-matrix-store-type (*)) store))
+ (type (real-matrix-store-type *) store))
(dotimes (k nxm mat)
(declare (type fixnum k))
(setf (aref store k) (funcall func (aref store k))))))
diff --git a/src/norm.lisp b/src/norm.lisp
index b03b512..7087dbe 100644
--- a/src/norm.lisp
+++ b/src/norm.lisp
@@ -128,7 +128,7 @@
(nxm (number-of-elements a))
(store (store a)))
(declare (type fixnum n m nxm)
- (type (real-matrix-store-type (*)) store))
+ (type (real-matrix-store-type *) store))
(if (row-or-col-vector-p a)
(case p
@@ -199,7 +199,7 @@
(nxm (number-of-elements a))
(store (store a)))
(declare (type fixnum n m nxm)
- (type (complex-matrix-store-type (*)) store))
+ (type (complex-matrix-store-type *) store))
(if (row-or-col-vector-p a)
(case p
diff --git a/src/real-matrix.lisp b/src/real-matrix.lisp
index f6a0e0e..87b6684 100644
--- a/src/real-matrix.lisp
+++ b/src/real-matrix.lisp
@@ -9,14 +9,13 @@
(deftype real-matrix-store-type (size)
"The type of the storage structure for a REAL-MATRIX"
- `(simple-array double-float ,size))
+ `(simple-array double-float (,size)))
)
-
;;
(defclass real-matrix (standard-matrix)
((store
:initform nil
- :type (simple-array real-matrix-element-type (*))))
+ :type (real-matrix-store-type *)))
(:documentation "A class of matrices with real elements."))
;;
@@ -34,13 +33,13 @@
;;
(defmethod matrix-ref-1d ((matrix real-matrix) (idx fixnum))
(let ((store (store matrix)))
- (declare (type (real-matrix-store-type (*)) store))
+ (declare (type (real-matrix-store-type *) store))
(aref store idx)))
(defmethod (setf matrix-ref-1d) ((value cl:real) (matrix real-matrix) (idx fixnum))
(let ((store (store matrix)))
- (declare (type (real-matrix-store-type (*)) store))
+ (declare (type (real-matrix-store-type *) store))
(setf (aref store idx) value)))
;;
@@ -103,7 +102,7 @@ don't know how to coerce COMPLEX to REAL"))
(size (* n m))
(store (allocate-real-store size)))
(declare (type fixnum n m size)
- (type (real-matrix-store-type (*)) store))
+ (type (real-matrix-store-type *) store))
(multiple-value-bind (row-stride col-stride)
(ecase order
(:row-major (values m 1))
@@ -126,7 +125,7 @@ don't know how to coerce COMPLEX to REAL"))
(size (* n m))
(store (allocate-real-store size)))
(declare (type fixnum n m size)
- (type (real-matrix-store-type (*)) store))
+ (type (real-matrix-store-type *) store))
(multiple-value-bind (row-stride col-stride)
(ecase order
(:row-major (values m 1))
diff --git a/src/realimag.lisp b/src/realimag.lisp
index 0b7ae11..2b35319 100644
--- a/src/realimag.lisp
+++ b/src/realimag.lisp
@@ -112,8 +112,8 @@
(store (store mat))
(new-store (allocate-real-store nxm)))
(declare (type fixnum n m nxm)
- (type (complex-matrix-store-type (*)) store)
- (type (real-matrix-store-type (*)) new-store))
+ (type (complex-matrix-store-type *) store)
+ (type (real-matrix-store-type *) new-store))
(dcopy nxm store 2 new-store 1)
@@ -140,8 +140,8 @@ its element types are unknown"))
(store (store mat))
(new-store (allocate-real-store nxm)))
(declare (type fixnum n m nxm)
- (type (complex-matrix-store-type (*)) store)
- (type (real-matrix-store-type (*)) new-store))
+ (type (complex-matrix-store-type *) store)
+ (type (real-matrix-store-type *) new-store))
(with-vector-data-addresses ((addr-store store)
(addr-new-store new-store))
diff --git a/src/ref.lisp b/src/ref.lisp
index ce21f79..c5619ef 100644
--- a/src/ref.lisp
+++ b/src/ref.lisp
@@ -152,7 +152,7 @@
(store (allocate-real-store k)))
(declare (optimize (speed 3) (safety 0))
(type fixnum n m k)
- (type (real-matrix-store-type (*)) idx-store mat-store store))
+ (type (real-matrix-store-type *) idx-store mat-store store))
(if (and (> n 1)
(> m 1))
@@ -183,7 +183,7 @@
(store (allocate-real-store k)))
(declare (optimize (speed 3) (safety 0))
(type fixnum n m k)
- (type (real-matrix-store-type (*)) mat-store store))
+ (type (real-matrix-store-type *) mat-store store))
(if (and (> n 1)
(> m 1))
@@ -216,7 +216,7 @@
(mat-store (store mat))
(store (store slice)))
(declare (type fixnum l n m)
- (type (real-matrix-store-type (*))
+ (type (real-matrix-store-type *)
row-idx-store
col-idx-store
mat-store
@@ -247,7 +247,7 @@
(store (store slice)))
(declare (type fixnum l n m)
- (type (real-matrix-store-type (*))
+ (type (real-matrix-store-type *)
mat-store
store))
@@ -269,7 +269,7 @@
(store (store mat)))
(declare (optimize (speed 3) (safety 0))
(type fixnum n m)
- (type (real-matrix-store-type (*)) store))
+ (type (real-matrix-store-type *) store))
(dotimes (i n)
(declare (type fixnum i))
@@ -303,7 +303,7 @@
(m (ncols matrix))
(store (store matrix)))
(declare (type fixnum n m)
- (type (real-matrix-store-type (*)) store))
+ (type (real-matrix-store-type *) store))
(labels ((consistent-i (i)
(and (integerp i)
@@ -472,7 +472,7 @@
(declare (optimize (speed 3) (safety 0))
(type fixnum n m k)
- (type (real-matrix-store-type (*)) idx-store new-store mat-store))
+ (type (real-matrix-store-type *) idx-store new-store mat-store))
(if (and (> n 1)
(> m 1))
@@ -496,7 +496,7 @@
(mat-store (store mat)))
(declare (optimize (speed 3) (safety 0))
(type fixnum n m)
- (type (real-matrix-store-type (*)) new-store mat-store))
+ (type (real-matrix-store-type *) new-store mat-store))
(if (and (> n 1)
(> m 1))
@@ -519,7 +519,7 @@
(new-store (store new))
(mat-store (store mat)))
(declare (type fixnum l n m)
- (type (real-matrix-store-type (*))
+ (type (real-matrix-store-type *)
row-idx-store
col-idx-store
new-store
@@ -545,7 +545,7 @@
(mat-store (store mat)))
(declare (type fixnum l n)
- (type (real-matrix-store-type (*))
+ (type (real-matrix-store-type *)
new-store
mat-store))
@@ -572,7 +572,7 @@
(declare (optimize (speed 3) (safety 0))
(type fixnum n m k)
- (type (real-matrix-store-type (*)) idx-store mat-store))
+ (type (real-matrix-store-type *) idx-store mat-store))
(if (and (> n 1)
(> m 1))
@@ -594,7 +594,7 @@
(mat-store (store mat)))
(declare (optimize (speed 3) (safety 0))
(type fixnum n m)
- (type (real-matrix-store-type (*)) mat-store))
+ (type (real-matrix-store-type *) mat-store))
(if (and (> n 1)
(> m 1))
@@ -618,7 +618,7 @@
(col-idx-store (store col-idx))
(mat-store (store mat)))
(declare (type fixnum l n m)
- (type (real-matrix-store-type (*))
+ (type (real-matrix-store-type *)
row-idx-store
col-idx-store
mat-store))
@@ -642,7 +642,7 @@
(mat-store (store mat)))
(declare (type fixnum l)
- (type (real-matrix-store-type (*))
+ (type (real-matrix-store-type *)
mat-store))
@@ -676,7 +676,7 @@
(new-store (store new)))
(declare (type fixnum n m new-n new-m)
- (type (real-matrix-store-type (*)) store new-store))
+ (type (real-matrix-store-type *) store new-store))
(let ((p (if (integerp i)
@@ -775,7 +775,7 @@
(new-store (store new)))
(declare (type fixnum n m new-n new-m)
- (type (real-matrix-store-type (*)) store new-store))
+ (type (real-matrix-store-type *) store new-store))
(let ((p (if (integerp i)
@@ -838,7 +838,7 @@
(new-store (store new)))
(declare (type fixnum n)
- (type (real-matrix-store-type (*)) store new-store))
+ (type (real-matrix-store-type *) store new-store))
(setf (aref store (fortran-matrix-indexing i 0 n))
(aref new-store 0))))
@@ -888,7 +888,7 @@
(new-store (store new)))
(declare (type fixnum n m new-n new-m)
- (type (real-matrix-store-type (*)) store new-store))
+ (type (real-matrix-store-type *) store new-store))
(let ((p (if (integerp i)
@@ -1025,7 +1025,7 @@
(store (store matrix)))
(declare (type fixnum n m)
- (type (real-matrix-store-type (*)) store))
+ (type (real-matrix-store-type *) store))
(labels ((consistent-i (i)
(and (integerp i)
@@ -1100,7 +1100,7 @@
(store (store matrix)))
(declare (type fixnum n m)
- (type (real-matrix-store-type (*)) store))
+ (type (real-matrix-store-type *) store))
(labels ((consistent-i (i)
(and (integerp i)
@@ -1134,7 +1134,7 @@
(store (store matrix)))
(declare (type fixnum n m)
- (type (real-matrix-store-type (*)) store))
+ (type (real-matrix-store-type *) store))
(labels ((consistent-i (i)
(and (integerp i)
@@ -1235,8 +1235,8 @@
(store (allocate-complex-store k)))
(declare (optimize (speed 3) (safety 0))
(type fixnum n m k)
- (type (real-matrix-store-type (*)) idx-store)
- (type (complex-matrix-store-type (*)) mat-store store))
+ (type (real-matrix-store-type *) idx-store)
+ (type (complex-matrix-store-type *) mat-store store))
(if (and (> n 1)
(> m 1))
@@ -1272,7 +1272,7 @@
(store (allocate-complex-store k)))
(declare (optimize (speed 3) (safety 0))
(type fixnum n m k)
- (type (complex-matrix-store-type (*)) mat-store store))
+ (type (complex-matrix-store-type *) mat-store store))
(if (and (> n 1)
(> m 1))
@@ -1311,10 +1311,10 @@
(mat-store (store mat))
(store (store slice)))
(declare (type fixnum l n m)
- (type (real-matrix-store-type (*))
+ (type (real-matrix-store-type *)
row-idx-store
col-idx-store)
- (type (complex-matrix-store-type (*))
+ (type (complex-matrix-store-type *)
mat-store
store))
@@ -1347,7 +1347,7 @@
(store (store slice)))
(declare (type fixnum l n m)
- (type (complex-matrix-store-type (*))
+ (type (complex-matrix-store-type *)
mat-store
store))
@@ -1375,7 +1375,7 @@
(m (ncols matrix))
(store (store matrix)))
(declare (type fixnum n m)
- (type (complex-matrix-store-type (*)) store))
+ (type (complex-matrix-store-type *) store))
(labels ((consistent-i (i)
(and (integerp i)
@@ -1452,7 +1452,7 @@
(m (ncols matrix))
(store (store matrix)))
(declare (type fixnum n m)
- (type (complex-matrix-store-type (*)) store))
+ (type (complex-matrix-store-type *) store))
(labels ((consistent-i (i)
(and (integerp i)
@@ -1498,7 +1498,7 @@
(m (ncols matrix))
(store (store matrix)))
(declare (type fixnum n m)
- (type (complex-matrix-store-type (*)) store))
+ (type (complex-matrix-store-type *) store))
(labels ((consistent-i (i)
(and (integerp i)
@@ -1595,8 +1595,8 @@
(declare (optimize (speed 3) (safety 0))
(type fixnum n m k)
- (type (real-matrix-store-type (*)) idx-store)
- (type (complex-matrix-store-type (*)) new-store mat-store))
+ (type (real-matrix-store-type *) idx-store)
+ (type (complex-matrix-store-type *) new-store mat-store))
(if (and (> n 1)
(> m 1))
@@ -1624,7 +1624,7 @@
(mat-store (store mat)))
(declare (optimize (speed 3) (safety 0))
(type fixnum n m) ;; k)
- (type (complex-matrix-store-type (*)) new-store mat-store))
+ (type (complex-matrix-store-type *) new-store mat-store))
(if (and (> n 1)
(> m 1))
@@ -1651,10 +1651,10 @@
(new-store (store new))
(mat-store (store mat)))
(declare (type fixnum l n m)
- (type (real-matrix-store-type (*))
+ (type (real-matrix-store-type *)
row-idx-store
col-idx-store)
- (type (complex-matrix-store-type (*))
+ (type (complex-matrix-store-type *)
new-store
mat-store))
@@ -1683,7 +1683,7 @@
(mat-store (store mat)))
(declare (type fixnum l n) ;;m)
- (type (complex-matrix-store-type (*))
+ (type (complex-matrix-store-type *)
new-store
mat-store))
@@ -1715,8 +1715,8 @@
(declare (optimize (speed 3) (safety 0))
(type fixnum n m k)
- (type (real-matrix-store-type (*)) idx-store new-store)
- (type (complex-matrix-store-type (*)) mat-store))
+ (type (real-matrix-store-type *) idx-store new-store)
+ (type (complex-matrix-store-type *) mat-store))
(if (and (> n 1)
(> m 1))
@@ -1744,8 +1744,8 @@
(mat-store (store mat)))
(declare (optimize (speed 3) (safety 0))
(type fixnum n m) ;; k)
- (type (real-matrix-store-type (*)) new-store)
- (type (complex-matrix-store-type (*)) mat-store))
+ (type (real-matrix-store-type *) new-store)
+ (type (complex-matrix-store-type *) mat-store))
(if (and (> n 1)
(> m 1))
@@ -1772,11 +1772,11 @@
(new-store (store new))
(mat-store (store mat)))
(declare (type fixnum l n m)
- (type (real-matrix-store-type (*))
+ (type (real-matrix-store-type *)
new-store
row-idx-store
col-idx-store)
- (type (complex-matrix-store-type (*))
+ (type (complex-matrix-store-type *)
mat-store))
(dotimes (i n)
@@ -1804,9 +1804,9 @@
(mat-store (store mat)))
(declare (type fixnum l n) ;; m)
- (type (real-matrix-store-type (*))
+ (type (real-matrix-store-type *)
new-store)
- (type (complex-matrix-store-type (*))
+ (type (complex-matrix-store-type *)
mat-store))
(let ((i -1))
@@ -1836,8 +1836,8 @@
(declare (optimize (speed 3) (safety 0))
(type fixnum n m k)
- (type (real-matrix-store-type (*)) idx-store)
- (type (complex-matrix-store-type (*)) mat-store))
+ (type (real-matrix-store-type *) idx-store)
+ (type (complex-matrix-store-type *) mat-store))
(if (and (> n 1)
(> m 1))
@@ -1864,7 +1864,7 @@
(mat-store (store mat)))
(declare (optimize (speed 3) (safety 0))
(type fixnum n m) ;; k)
- (type (complex-matrix-store-type (*)) mat-store))
+ (type (complex-matrix-store-type *) mat-store))
(if (and (> n 1)
(> m 1))
@@ -1892,10 +1892,10 @@
(col-idx-store (store col-idx))
(mat-store (store mat)))
(declare (type fixnum l n m)
- (type (real-matrix-store-type (*))
+ (type (real-matrix-store-type *)
row-idx-store
col-idx-store)
- (type (complex-matrix-store-type (*))
+ (type (complex-matrix-store-type *)
mat-store))
(dotimes (i n)
@@ -1922,7 +1922,7 @@
(mat-store (store mat)))
(declare (type fixnum l) ;; n m)
- (type (complex-matrix-store-type (*))
+ (type (complex-matrix-store-type *)
mat-store))
;; Hmm, another redundant i,j (see above)
@@ -1954,7 +1954,7 @@
(new-store (store new)))
(declare (type fixnum n m new-n new-m)
- (type (complex-matrix-store-type (*)) store new-store))
+ (type (complex-matrix-store-type *) store new-store))
(let ((p (if (integerp i)
@@ -2059,7 +2059,7 @@
(new-store (store new)))
(declare (type fixnum n m new-n new-m)
- (type (complex-matrix-store-type (*)) store new-store))
+ (type (complex-matrix-store-type *) store new-store))
(let ((p (if (integerp i)
@@ -2126,7 +2126,7 @@
(new-store (store new)))
(declare (type fixnum n m new-n new-m)
- (type (complex-matrix-store-type (*)) store new-store))
+ (type (complex-matrix-store-type *) store new-store))
(let ((p (if (integerp i)
@@ -2239,8 +2239,8 @@
(new-store (store new)))
(declare (type fixnum n m new-n new-m)
- (type (real-matrix-store-type (*)) new-store)
- (type (complex-matrix-store-type (*)) store))
+ (type (real-matrix-store-type *) new-store)
+ (type (complex-matrix-store-type *) store))
(let ((p (if (integerp i)
@@ -2344,8 +2344,8 @@
(new-store (store new)))
(declare (type fixnum n m new-n new-m)
- (type (real-matrix-store-type (*)) new-store)
- (type (complex-matrix-store-type (*)) store))
+ (type (real-matrix-store-type *) new-store)
+ (type (complex-matrix-store-type *) store))
(let ((p (if (integerp i)
@@ -2395,8 +2395,8 @@
(new-store (store new)))
(declare (type fixnum n m new-n new-m)
- (type (real-matrix-store-type (*)) new-store)
- (type (complex-matrix-store-type (*)) store))
+ (type (real-matrix-store-type *) new-store)
+ (type (complex-matrix-store-type *) store))
(let ((p (if (integerp i)
@@ -2494,7 +2494,7 @@
(declare (type fixnum n m)
- (type (complex-matrix-store-type (*)) store))
+ (type (complex-matrix-store-type *) store))
#+(or :ccl :allegro) (setq new (complex-coerce new))
@@ -2583,7 +2583,7 @@
(declare (type fixnum n m)
- (type (complex-matrix-store-type (*)) store))
+ (type (complex-matrix-store-type *) store))
#+(or :ccl :allegro) (setq new (complex-coerce new))
@@ -2626,7 +2626,7 @@
(declare (type fixnum n m)
- (type (complex-matrix-store-type (*)) store))
+ (type (complex-matrix-store-type *) store))
#+(or :ccl :allegro) (setq new (complex-coerce new))
@@ -2811,7 +2811,7 @@
(m (ncols item))
(store (store item)))
(declare (type fixnum n m)
- (type (real-matrix-store-type (*)) store))
+ (type (real-matrix-store-type *) store))
(dotimes (i n)
(declare (type fixnum i))
(dotimes (j m)
diff --git a/src/reshape.lisp b/src/reshape.lisp
index 7321600..9681990 100644
--- a/src/reshape.lisp
+++ b/src/reshape.lisp
@@ -115,7 +115,7 @@
(new-size (* new-nrows new-ncols))
(new-store (allocate-real-store new-size)))
(declare (fixnum old-size new-size)
- (type (real-matrix-store-type (*)) new-store))
+ (type (real-matrix-store-type *) new-store))
(dcopy (min old-size new-size) (store mat) 1 new-store 1)
@@ -132,7 +132,7 @@
;; it's bigger than the old size. Allocate it and copy the
;; elements over.
(let ((new-store (allocate-real-store new-size)))
- (declare (type (real-matrix-store-type (*)) new-store))
+ (declare (type (real-matrix-store-type *) new-store))
(dcopy new-size (store mat) 1 new-store 1)
(setf (slot-value mat 'store) new-store)))
@@ -148,7 +148,7 @@
(new-size (* new-nrows new-ncols))
(new-store (allocate-complex-store new-size)))
(declare (fixnum old-size new-size)
- (type (complex-matrix-store-type (*)) new-store))
+ (type (complex-matrix-store-type *) new-store))
(zcopy (min old-size new-size) (store mat) 1 new-store 1)
(make-instance 'complex-matrix :nrows new-nrows :ncols new-ncols :store new-store)))
@@ -164,7 +164,7 @@
;; it's bigger than the old size. Allocate it and copy the
;; elements over.
(let ((new-store (allocate-complex-store new-size)))
- (declare (type (complex-matrix-store-type (*)) new-store))
+ (declare (type (complex-matrix-store-type *) new-store))
(zcopy new-size (store mat) 1 new-store 1)
(setf (slot-value mat 'store) new-store)))
diff --git a/src/scal.lisp b/src/scal.lisp
index 6d2026e..8548a26 100644
--- a/src/scal.lisp
+++ b/src/scal.lisp
@@ -69,34 +69,34 @@
(in-package "MATLISP")
-#+nil (use-package "BLAS")
-#+nil (use-package "LAPACK")
-#+nil (use-package "FORTRAN-FFI-ACCESSORS")
-
-#+nil (export '(scal!
- scal))
-
-(defgeneric scal (alpha x)
- (:documentation
-"
- Sytnax
- ======
- (SCAL alpha x)
-
- Purpose
- =======
- Computes and returns a new matrix equal to
-
- alpha * X
-
- where alpha is a scalar and X is a matrix.
-
-"))
-
+(defmacro generate-typed-scal!-func (func element-type store-type matrix-type blas-func)
+ `(defun ,func (alpha mat-x)
+ (declare (type ,matrix-type mat-x)
+ (type ,element-type alpha)
+ (optimize (safety 0) (speed 3)))
+ (mlet* (((cp-x inc-x sz-x) (blas-copyable-p mat-x)
+ :type (boolean fixnum fixnum))
+ ((hd-x st-x) (slot-values mat-x '(head store))
+ :type (fixnum (,store-type *))))
+ (if cp-x
+ (,blas-func sz-x alpha st-x inc-x :head-x hd-x)
+ (symbol-macrolet
+ ((common-code
+ (mlet* (((nr-x nc-x rs-x cs-x) (slot-values mat-x '(number-of-rows number-of-cols row-stride col-stride))
+ :type (fixnum fixnum fixnum fixnum)))
+ (loop for i from 0 below nr-x
+ do (,blas-func nc-x alpha st-x cs-x :head-x (+ hd-x (* i rs-x)))))))
+ (if (> (nrows mat-x) (ncols mat-x))
+ (with-transpose! (mat-x)
+ common-code)
+ common-code)))
+ mat-x)))
+
+;;
(defgeneric scal! (alpha x)
(:documentation
"
- Sytnax
+ Syntax
======
(SCAL! alpha x)
@@ -106,117 +106,66 @@
stored in X.
"))
-(defmethod scal ((alpha number) (x number))
- (* alpha x))
-
-(defmethod scal ((alpha #+(or cmu sbcl) double-float #-(or cmu sbcl) float) (x real-matrix))
- (let ((nxm (number-of-elements x))
- (result (copy x)))
- (declare (type fixnum nxm))
-
- (dscal nxm alpha (store result) 1)
- result))
-
-(defmethod scal ((alpha cl:real) (x real-matrix))
- (scal (coerce alpha 'real-matrix-element-type) x))
-
-(defmethod scal ((alpha #+:cmu kernel::complex-double-float
- #+:sbcl sb-kernel::complex-double-float
- #-(or cmu sbcl) complex) (x real-matrix))
- (let* ((nxm (number-of-elements x))
- (n (nrows x))
- (m (ncols x))
- (result (make-complex-matrix-dim n m)))
- (declare (type fixnum n m nxm))
-
- #-(or cmu sbcl) (setq alpha (complex-coerce alpha))
-
- (copy! x result)
- (setf (aref *1x1-complex-array* 0) (realpart alpha))
- (setf (aref *1x1-complex-array* 1) (imagpart alpha))
- (zscal nxm *1x1-complex-array* (store result) 1)
-
- result))
-
-#+(or :cmu :sbcl)
-(defmethod scal ((alpha complex) (x real-matrix))
- (scal (complex-coerce alpha) x))
-
-(defmethod scal ((alpha #+(or cmu sbcl) double-float #-(or cmu sbcl) float) (x complex-matrix))
- (let ((nxm (number-of-elements x))
- (result (copy x)))
- (declare (type fixnum nxm))
- (zdscal nxm alpha (store result) 1)
-
- result))
-
-(defmethod scal ((alpha cl:real) (x complex-matrix))
- (scal (coerce alpha 'real-matrix-element-type) x))
-
-(defmethod scal ((alpha #+:cmu kernel::complex-double-float
- #+:sbcl sb-kernel::complex-double-float
- #-(or cmu sbcl) complex) (x complex-matrix))
- (let ((nxm (number-of-elements x))
- (result (copy x)))
- (declare (type fixnum nxm))
-
- #-(or cmu sbcl) (setq alpha (complex-coerce alpha))
-
- (setf (aref *1x1-complex-array* 0) (realpart alpha))
- (setf (aref *1x1-complex-array* 1) (imagpart alpha))
- (zscal nxm *1x1-complex-array* (store result) 1)
-
- result))
-
-#+(or :cmu :sbcl)
-(defmethod scal ((alpha complex) (x complex-matrix))
- (scal (complex-coerce alpha) x))
-
+;;
+(generate-typed-scal!-func real-double-dscal!-typed double-float real-matrix-store-type real-matrix blas:dscal)
(defmethod scal! ((alpha number) (x number))
- (error "cannot SCAL! two scalars, arg X must
+ (error "Cannot SCAL! two scalars, arg X must
be a matrix to SCAL!"))
-(defmethod scal! ((alpha #+(or cmu sbcl) double-float #-(or cmu sbcl) float) (x real-matrix))
- (let ((nxm (number-of-elements x)))
- (declare (type fixnum nxm))
-
- (dscal nxm alpha (store x) 1)
- x))
+(defmethod scal! ((alpha complex) (x real-matrix))
+ (error "Cannot SCAL! a REAL-MATRIX by a COMPLEX, don't know
+how to coerce COMPLEX to REAL"))
(defmethod scal! ((alpha cl:real) (x real-matrix))
- (scal! (coerce alpha 'real-matrix-element-type) x))
+ (real-double-dscal!-typed (coerce alpha 'double-float) x))
-(defmethod scal! ((alpha complex) (x real-matrix))
- (error "cannot SCAL! a REAL-MATRIX by a COMPLEX, don't know
-how to coerce COMPLEX to REAL"))
+;;
+(generate-typed-scal!-func complex-double-dscal!-typed double-float complex-matrix-store-type complex-matrix blas:zdscal)
-(defmethod scal! ((alpha #+(or cmu sbcl) double-float #-(or cmu sbcl) float) (x complex-matrix))
- (let ((nxm (number-of-elements x)))
- (declare (type fixnum nxm))
- (zdscal nxm alpha (store x) 1)
-
- x))
+(generate-typed-scal!-func complex-double-zscal!-typed (complex (double-float * *)) complex-matrix-store-type complex-matrix blas:zscal)
(defmethod scal! ((alpha cl:real) (x complex-matrix))
- (scal! (coerce alpha 'real-matrix-element-type) x))
+ (complex-double-dscal!-typed (coerce alpha 'double-float) x))
+
+(defmethod scal! ((alpha complex) (x complex-matrix))
+ (complex-double-zscal!-typed (complex-coerce alpha) x))
+
+;;;;
+(defgeneric scal (alpha x)
+ (:documentation
+"
+ Syntax
+ ======
+ (SCAL alpha x)
-(defmethod scal! ((alpha #+:cmu kernel::complex-double-float
- #+:sbcl sb-kernel::complex-double-float
- #-(or cmu sbcl) complex) (x complex-matrix))
- (let ((nxm (number-of-elements x)))
- (declare (type fixnum nxm))
+ Purpose
+ =======
+ Computes and returns a new matrix equal to
- #-(or cmu sbcl) (setq alpha (complex-coerce alpha))
+ alpha * X
- (setf (aref *1x1-complex-array* 0) (realpart alpha))
- (setf (aref *1x1-complex-array* 1) (imagpart alpha))
- (zscal nxm *1x1-complex-array* (store x) 1)
+ where alpha is a scalar and X is a matrix.
- x))
+"))
-#+(or :cmu :sbcl)
-(defmethod scal! ((alpha complex) (x complex-matrix))
- (scal! (complex-coerce alpha) x))
+(defmethod scal ((alpha number) (x number))
+ (* alpha x))
+
+;;
+(defmethod scal ((alpha cl:real) (x real-matrix))
+ (let ((result (copy x)))
+ (scal! alpha result)))
+(defmethod scal ((alpha complex) (x real-matrix))
+ (let* ((n (nrows x))
+ (m (ncols x))
+ (result (make-complex-matrix-dim n m)))
+ (declare (type fixnum n m))
+ (copy! x result)
+ (scal! alpha result)))
+;;
+(defmethod scal ((alpha number) (x complex-matrix))
+ (let ((result (copy x)))
+ (scal! alpha result)))
\ No newline at end of file
diff --git a/src/special.lisp b/src/special.lisp
index 93ff726..20138bc 100644
--- a/src/special.lisp
+++ b/src/special.lisp
@@ -153,7 +153,7 @@
(unity #.(coerce 1 'real-matrix-element-type)))
(declare (fixnum size)
- (type (real-matrix-store-type (*)) store))
+ (type (real-matrix-store-type *) store))
(dotimes (k size)
(declare (fixnum k))
(setf (aref store k) (random unity state)))
diff --git a/src/standard-matrix.lisp b/src/standard-matrix.lisp
index 142639e..47c4ff9 100644
--- a/src/standard-matrix.lisp
+++ b/src/standard-matrix.lisp
@@ -307,4 +307,27 @@ matrix and a number"))
:nrows nrows :ncols ncols
:store (store matrix)
:head (store-indexing i j hd rs cs)
- :row-stride rs :col-stride cs)))
\ No newline at end of file
+ :row-stride rs :col-stride cs)))
+
+;;
+(defmacro with-transpose! (matlst &rest body)
+ `(progn
+ ,@(mapcar #'(lambda (mat) `(transpose! ,mat)) matlst)
+ ,@body
+ ,@(mapcar #'(lambda (mat) `(transpose! ,mat)) matlst)))
+
+;;
+(defun blas-copyable-p (matrix)
+ (declare (optimize (safety 0) (speed 3))
+ (type (or real-matrix complex-matrix) matrix))
+ (mlet* ((nr (nrows matrix) :type fixnum)
+ (nc (ncols matrix) :type fixnum)
+ (rs (row-stride matrix) :type fixnum)
+ (cs (col-stride matrix) :type fixnum)
+ (ne (number-of-elements matrix) :type fixnum))
+ (cond
+ ((= nc 1) (values t rs ne))
+ ((= nr 1) (values t cs ne))
+ ((= rs (* nc cs)) (values t cs ne))
+ ((= cs (* nr rs)) (values t rs ne))
+ (t (values nil -1 -1)))))
\ No newline at end of file
diff --git a/src/sum.lisp b/src/sum.lisp
index ad05d8f..49d0907 100644
--- a/src/sum.lisp
+++ b/src/sum.lisp
@@ -87,7 +87,7 @@
(nxm (number-of-elements a)))
(declare (type fixnum nxm)
(type real-matrix-element-type result)
- (type (real-matrix-store-type (*)) store))
+ (type (real-matrix-store-type *) store))
(dotimes (i nxm)
(declare (type fixnum i))
(incf result (aref store i)))
@@ -99,7 +99,7 @@
(store-a (store a))
(store-result (store result)))
(declare (type fixnum n m)
- (type (real-matrix-store-type (*)) store-a store-result))
+ (type (real-matrix-store-type *) store-a store-result))
(dotimes (i n)
(declare (type fixnum i))
(setf (aref store-result i)
@@ -118,7 +118,7 @@
(nxm (number-of-elements a)))
(declare (type fixnum nxm)
(type complex-matrix-element-type imagpart realpart)
- (type (complex-matrix-store-type (*)) store))
+ (type (complex-matrix-store-type *) store))
(dotimes (i nxm)
(declare (type fixnum i))
(let ((k (* 2 i)))
@@ -140,7 +140,7 @@
(store-a (store a))
(store-result (store result)))
(declare (type fixnum n m)
- (type (complex-matrix-store-type (*)) store-a store-result))
+ (type (complex-matrix-store-type *) store-a store-result))
(dotimes (i n)
(declare (type fixnum i))
(let ((realpart 0.0d0)
diff --git a/src/trans.lisp b/src/trans.lisp
index 48b6496..192d548 100644
--- a/src/trans.lisp
+++ b/src/trans.lisp
@@ -141,7 +141,7 @@
(new-store (allocate-real-store nxm)))
(declare (type fixnum n m nxm)
- (type (real-matrix-store-type (*)) store new-store))
+ (type (real-matrix-store-type *) store new-store))
(dotimes (i n)
(declare (type fixnum i))
@@ -160,7 +160,7 @@
(new-store (allocate-complex-store nxm)))
(declare (type fixnum n m nxm)
- (type (complex-matrix-store-type (*)) store new-store))
+ (type (complex-matrix-store-type *) store new-store))
(dotimes (i n)
(declare (type fixnum i))
@@ -201,7 +201,7 @@ a STANDARD-MATRIX, element types are not known"))
(new-store (allocate-complex-store nxm)))
(declare (type fixnum n m nxm)
- (type (complex-matrix-store-type (*)) store new-store))
+ (type (complex-matrix-store-type *) store new-store))
(dotimes (i n)
(declare (type fixnum i))
diff --git a/src/utilities.lisp b/src/utilities.lisp
index 2dcbf64..76caf8c 100644
--- a/src/utilities.lisp
+++ b/src/utilities.lisp
@@ -1,6 +1,63 @@
(in-package :utilities)
;;
+(defmacro mlet* (decls &rest body)
+" mlet* ({ {(var*) | var} values-form &keyform declare type}*) form*
+
+ o var is just one symbol -> expands into let
+ o var is a list -> expands into multiple-value-bind
+
+ This macro also handles type declarations.
+
+ Example:
+ (mlet* ((x 2 :type fixnum :declare ((optimize (safety 0) (speed 3))))
+ ((a b) (floor 3) :type (nil fixnum)))
+ (+ x b))
+
+ expands into:
+
+ (let ((x 2))
+ (declare (optimize (safety 0) (speed 3))
+ (type fixnum x))
+ (multiple-value-bind (a b)
+ (floor 3)
+ (declare (ignore a)
+ (type fixnum b))
+ (+ x b)))
+"
+ (labels ((mlet-decl (vars type decls)
+ (when (or type decls)
+ `((declare ,@decls
+ ,@(when type
+ (mapcar #'(lambda (tv) (if (null (first tv))
+ `(ignore ,(second tv))
+ `(type ,(first tv) ,(second tv))))
+ (map 'list #'list type vars)))))))
+ ;;
+ (mlet-transform (elst nest-code)
+ (destructuring-bind (vars form &key declare type) elst
+ `(,(append (cond
+ ;;If there is only one element use let
+ ;;instead of multiple-value-bind
+ ((or (symbolp vars) (null (cdr vars)))
+ `(let ((,(car (ensure-list vars)) ,form))))
+ ;;
+ (t
+ `(multiple-value-bind (,@vars) ,form)))
+ (mlet-decl (ensure-list vars) (ensure-list type) declare)
+ nest-code))))
+ ;;
+ (mlet-walk (elst body)
+ (if (null elst)
+ `(,@body)
+ (mlet-transform (car elst) (mlet-walk (cdr elst) body)))))
+ ;;
+ (if decls
+ (car (mlet-walk decls body))
+ `(progn
+ ,@body))))
+
+;;
(defmacro with-gensyms (symlist &body body)
`(let ,(mapcar #'(lambda (sym)
`(,sym (gensym ,(symbol-name sym))))
@@ -55,64 +112,6 @@
lst
`(,lst)))
-(defmacro mlet* (decls &rest body)
-"
-mlet* ({ {(var*) | var} values-form &keyform declare type}*) form*
-
-o var is just one symbol -> expands into let
-o var is a list -> expands into multiple-value-bind
-
-This macro also handles type declarations.
-
-Example:
-(mlet* ((x 2 :type fixnum :declare ((optimize (safety 0) (speed 3))))
- ((a b) (floor 3) :type (nil fixnum)))
- (+ x b))
-
-expands into:
-
-(let ((x 2))
- (declare (optimize (safety 0) (speed 3))
- (type fixnum x))
- (multiple-value-bind (a b)
- (floor 3)
- (declare (ignore a)
- (type fixnum b))
- (+ x b)))
-"
- (labels ((mlet-decl (vars type decls)
- (when (or type decls)
- `((declare ,@decls
- ,@(when type
- (mapcar #'(lambda (tv) (if (null (first tv))
- `(ignore ,(second tv))
- `(type ,(first tv) ,(second tv))))
- (map 'list #'list type vars)))))))
- ;;
- (mlet-transform (elst nest-code)
- (destructuring-bind (vars form &key declare type) elst
- `(,(append (cond
- ;;If there is only one element use let
- ;;instead of multiple-value-bind
- ((or (symbolp vars) (null (cdr vars)))
- `(let ((,(car (ensure-list vars)) ,form))))
- ;;
- (t
- `(multiple-value-bind (,@vars) ,form)))
- (mlet-decl (ensure-list vars) (ensure-list type) declare)
- nest-code))))
- ;;
- (mlet-walk (elst body)
- (if (null elst)
- `(,@body)
- (mlet-transform (car elst) (mlet-walk (cdr elst) body)))))
- ;;
- (if decls
- (car (mlet-walk decls body))
- `(progn
- ,@body))))
-
-;;
(defmacro if-ret (form &rest else-body)
"if-ret (form &rest else-body)
Evaluate form, and if the form is not nil, then return it,
@@ -155,4 +154,6 @@ else run else-body"
,@(app-equal (cdr lst))))))
`(let ((,key-eval ,keyform))
(cond
- ,@(app-equal body))))))
\ No newline at end of file
+ ,@(app-equal body))))))
+
+;;
commit b0f1cb2dd42338c9189c83cbcbcb177eaf1c7845
Author: Akshay Srinivasan <aks...@gm...>
Date: Mon Mar 12 10:07:12 2012 +0530
Added slot-values to the list of symbols exported by utilities
diff --git a/packages.lisp b/packages.lisp
index f95b7ac..383af89 100644
--- a/packages.lisp
+++ b/packages.lisp
@@ -162,6 +162,7 @@
#:get-arg
#:nconsc
#:with-gensyms
+ #:slot-values
#:mlet*))
(defpackage :fortran-ffi-accessors
commit 98cdecd68f57b4a561ac8f68a0ede4a0374a6a95
Author: Akshay Srinivasan <aks...@gm...>
Date: Sun Mar 11 23:58:49 2012 +0530
-> Wrote macros to generate a typed copy! function.
Haven't yet adapted the generic methods yet.
diff --git a/src/blas.lisp b/src/blas.lisp
index f8adeed..999b389 100644
--- a/src/blas.lisp
+++ b/src/blas.lisp
@@ -151,9 +151,9 @@ t;;; -*- Mode: lisp; Syntax: ansi-common-lisp; Package: :blas; Base: 10 -*-
Y(0),Y(INCY), ... , Y((N-1)*INCY)
"
(n :integer :input)
- (dx (* :double-float :inc head-dx))
+ (dx (* :double-float :inc head-x))
(incx :integer :input)
- (dy (* :double-float :inc head-dy) :output)
+ (dy (* :double-float :inc head-y) :output)
(incy :integer :input)
)
@@ -328,9 +328,9 @@ t;;; -*- Mode: lisp; Syntax: ansi-common-lisp; Package: :blas; Base: 10 -*-
Y(0),Y(2*INCY), ... , Y(2*(N-1)*INCY)
"
(n :integer :input)
- (zx (* :complex-double-float :inc head-zx))
+ (zx (* :complex-double-float :inc head-x))
(incx :integer :input)
- (zy (* :complex-double-float :inc head-zy) :output)
+ (zy (* :complex-double-float :inc head-y) :output)
(incy :integer :input)
)
diff --git a/src/complex-matrix.lisp b/src/complex-matrix.lisp
index 28785bc..48e7d0f 100644
--- a/src/complex-matrix.lisp
+++ b/src/complex-matrix.lisp
@@ -284,4 +284,21 @@
"cannot make a ~A x ~A matrix" n m)
(make-complex-matrix-dim n m)))
(t
- (error "require 1 or 2 arguments to make a matrix")))))
\ No newline at end of file
+ (error "require 1 or 2 arguments to make a matrix")))))
+
+;;
+
+(defun realpart! (mat)
+ (cond
+ ((typep mat 'real-matrix) mat)
+ ((typep mat 'complex-matrix) (make-instance 'real-matrix :store (store mat)
+ :nrows (nrows mat) :ncols (ncols mat)
+ :row-stride (* 2 (row-stride mat)) :col-stride (* 2 (col-stride mat))
+ :head (* 2 (head mat))))))
+
+(defun imagpart! (mat)
+ (cond
+ ((typep mat 'complex-matrix) (make-instance 'real-matrix :store (store mat)
+ :nrows (nrows mat) :ncols (ncols mat)
+ :row-stride (* 2 (row-stride mat)) :col-stride (* 2 (col-stride mat))
+ :head (+ 1 (* 2 (head mat)))))))
\ No newline at end of file
diff --git a/src/copy.lisp b/src/copy.lisp
index 264ce46..586aaf0 100644
--- a/src/copy.lisp
+++ b/src/copy.lisp
@@ -78,6 +78,7 @@
(in-package "MATLISP")
+;;
(defun blas-copyable-p (matrix)
(declare (optimize (safety 0) (speed 3))
(type (or real-matrix complex-matrix) matrix))
@@ -91,54 +92,45 @@
((= nr 1) (values t cs ne))
((= rs (* nc cs)) (values t cs ne))
((= cs (* nr rs)) (values t rs ne))
- (t nil))))
-
+ (t (values nil -1 -1)))))
+;;
+(defmacro with-transpose! (matlst &rest body)
+ `(progn
+ ,@(mapcar #'(lambda (mat) `(transpose! ,mat)) matlst)
+ ,@body
+ ,@(mapcar #'(lambda (mat) `(transpose! ,mat)) matlst)))
-(defmacro generate-typed-copy!-func (func)
- `(defun ,func (matrix-a matrix-b)
- (declare (optimize (safety 0) (speed 3))
- (type (or ,matrix-type matrix-a matrix-b)))))
+;;
+(defmacro generate-typed-copy!-func (func store-type matrix-type blas-func)
+ `(defun ,func (mat-a mat-b)
+ (declare (type ,matrix-type mat-a mat-b)
+ (optimize (safety 0) (speed 3)))
+ (mlet* (((cp-a inc-a sz-a) (blas-copyable-p mat-a) :type (boolean fixnum nil))
+ ((cp-b inc-b sz-b) (blas-copyable-p mat-b) :type (boolean fixnum nil))
+ ((hd-a st-a sz) (slot-values mat-a '(head store number-of-elements)) :type (fixnum (,store-type *) fixnum))
+ ((hd-b st-b) (slot-values mat-b '(head store)) :type (fixnum (,store-type *))))
+ (if (and cp-a cp-b)
+ (,blas-func sz (store mat-a) inc-a (store mat-b) inc-b :head-x hd-a :head-y hd-b)
+ (symbol-macrolet
+ ((common-code
+ (mlet* (((nr-a nc-a rs-a cs-a) (slot-values mat-a '(number-of-rows number-of-cols row-stride col-stride))
+ :type (fixnum fixnum fixnum fixnum))
+ ((rs-b cs-b) (slot-values mat-b '(row-stride col-stride))
+ :type (fixnum fixnum)))
+ (loop for i from 0 below nr-a
+ do (,blas-func nc-a st-a cs-a st-b cs-b :head-x (+ hd-a (* i rs-a)) :head-y (+ hd-b (* i rs-b)))))))
+ ;;Choose the smaller of the loops
+ (if (> (nrows mat-a) (ncols mat-a))
+ (with-transpose! (mat-a mat-b)
+ common-code)
+ common-code)))
+ mat-b)))
;;
(defvar *1x1-real-array* (make-array 1 :element-type 'double-float))
(defvar *1x1-complex-array* (make-array 2 :element-type 'double-float))
-(defgeneric copy (matrix)
- (:documentation
- "
- Syntax
- ======
- (COPY x)
-
- Purpose
- =======
- Return a copy of the matrix X"))
-
-(defmethod copy ((matrix standard-matrix))
- (make-instance 'standard-matrix :nrows (nrows matrix) :ncols (ncols matrix) :store (copy-seq (store matrix))))
-
-(defmethod copy ((matrix real-matrix))
- (let* ((size (number-of-elements matrix))
- (n (nrows matrix))
- (m (ncols matrix))
- (result (make-real-matrix-dim n m)))
- (declare (type fixnum size n m))
- (blas:dcopy size (store matrix) 1 (store result) 1)
- result))
-
-(defmethod copy ((matrix complex-matrix))
- (let* ((size (number-of-elements matrix))
- (n (nrows matrix))
- (m (ncols matrix))
- (result (make-complex-matrix-dim n m)))
- (declare (type fixnum size n m))
- (blas:zcopy size (store matrix) 1 (store result) 1)
- result))
-
-(defmethod copy ((matrix number))
- matrix)
-
;;
(defgeneric copy! (matrix new-matrix)
(:documentation
@@ -164,7 +156,6 @@
REAL-MATRIX but the converse is possible.
"))
-
(defmethod copy! :before ((x standard-matrix) (y standard-matrix))
(let ((nxm-x (number-of-elements x))
(nxm-y (number-of-elements y)))
@@ -172,6 +163,9 @@
(if (not (= nxm-x nxm-y))
(warn "arguments X,Y to COPY! are of different sizes"))))
+;;
+(generate-typed-copy!-func real-double-copy!-typed real-matrix-store-type real-matrix blas:dcopy)
+
(defmethod copy! ((x real-matrix) (y real-matrix))
(let* ((nxm-x (number-of-elements x))
(nxm-y (number-of-elements y))
@@ -197,6 +191,10 @@
(error "cannot copy a COMPLEX-MATRIX into a REAL-MATRIX,
don't know how to coerce a COMPLEX to a REAL"))
+
+;;
+(genera...
[truncated message content] |
|
From: Raymond T. <rt...@us...> - 2012-03-11 17:49:28
|
This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "matlisp".
The branch, master has been updated
via 14707c10ced9ee08431313323cc1905ac77e2a21 (commit)
via 7f45afca17c2501c0b2a777809e1998683e5a17f (commit)
via 8708b5e4f149f2ba16b879bfa87aa3085097cb34 (commit)
from 4d23b62fbced07e732614e290a08e90a731942ab (commit)
Those revisions listed above that are new to this repository have
not appeared on any other notification email; so we list those
revisions in full, below.
- Log -----------------------------------------------------------------
commit 14707c10ced9ee08431313323cc1905ac77e2a21
Author: Raymond Toy <toy...@gm...>
Date: Sun Mar 11 10:49:04 2012 -0700
Just rearranging the code to group related things together.
diff --git a/configure.ac b/configure.ac
index 41a9f22..fad1307 100644
--- a/configure.ac
+++ b/configure.ac
@@ -59,7 +59,10 @@ AC_F77_LIBRARY_LDFLAGS
# This is only for autoconf 2.49d or later!
AC_F77_FUNC(f77_name)
-# Setup our environment variables based on the value of f77_name
+# Setup our environment variables based on the value of f77_name. For
+# the record, F77_EXTRA_UNDERSCORE means that any embedded underscores
+# in the Fortran name are converted to two underscores for the C
+# equivalent name.
F77_LOWER_CASE=t
F77_UNDERSCORE=t
F77_EXTRA_UNDERSCORE=nil
@@ -208,7 +211,8 @@ FLIBS_OPTS=`echo $FLIBS | tr -s ' ' '\n' | $GREP -v gfortranbegin | $GREP -v x86
GIT_VERSION=`git describe --dirty`
-# Check that quicklisp exists. We need that currently to get cffi.
+# Check that quicklisp exists. We need that currently to get cffi if
+# it's not already available.
AC_CHECK_FILE($HOME/quicklisp/setup.lisp, HAVE_QL="(and)", HAVE_QL="(or)")
# Are we using a 32-bit lisp? If so, define LIB32 so that we compile
@@ -291,49 +295,6 @@ EOF
AC_MSG_RESULT([$NEED_FF2C])
fi
-dnl The following variables will be substituted into the .in files
-AC_SUBST(LISPEXEC)
-AC_SUBST(BLAS_OBJS)
-AC_SUBST(NO_ATLAS_LAPACK_OBJS)
-AC_SUBST(ATLAS_DIR)
-AC_SUBST(ATLAS_LIBS)
-AC_SUBST(ATLAS_P)
-AC_SUBST(F77_LOWER_CASE)
-AC_SUBST(F77_UNDERSCORE)
-AC_SUBST(F77_EXTRA_UNDERSCORE)
-AC_SUBST(LISPEXEC)
-AC_SUBST(LISPEVAL)
-AC_SUBST(srcdir)
-AC_SUBST(libdir)
-AC_SUBST(FLIBS_LIBS)
-AC_SUBST(FLIBS_OPTS)
-AC_SUBST(share_ext)
-AC_SUBST(GIT_VERSION)
-AC_SUBST(HAVE_QL)
-AC_SUBST(F2C)
-
-echo host = $host
-
-# Set the extension for shared libraries. This is not very robust.
-case $host in
- *darwin*) share_ext=dylib ;;
- *) share_ext=so ;;
-esac
-
-AC_CONFIG_FILES([
- matlisp.mk
- Makefile
- start.lisp
- config.lisp
- lib/lazy-loader.lisp
- LAPACK/SRC/Makefile
- LAPACK/BLAS/SRC/Makefile
- dfftpack/Makefile
- lib-src/toms715/Makefile
- lib-src/odepack/Makefile
- src/f77-mangling.lisp
-])
-
# Allow user to use ATLAS if available.
# We assume the standard names for the ATLAS libraries.
AC_ARG_WITH([atlas],
@@ -355,7 +316,10 @@ AM_CONDITIONAL([ATLAS], [test x$atlas = xtrue])
# Check to see if the ATLAS libraries are compatible with matlisp's
# ffi. Basically the same test as above that checks to see if -ff2c
-# is needed.
+# is needed. We call zdotu which is a Fortran function returning a
+# complex number. Matlisp assumes such functions return the result by
+# storing the answer at address given by a hidden first parameter to
+# the function.
if test x"$atlas" = xtrue; then
AC_MSG_CHECKING([if ATLAS is compatible with matlisp's FFI])
@@ -421,6 +385,49 @@ EOF
fi
fi
+dnl The following variables will be substituted into the .in files
+AC_SUBST(LISPEXEC)
+AC_SUBST(BLAS_OBJS)
+AC_SUBST(NO_ATLAS_LAPACK_OBJS)
+AC_SUBST(ATLAS_DIR)
+AC_SUBST(ATLAS_LIBS)
+AC_SUBST(ATLAS_P)
+AC_SUBST(F77_LOWER_CASE)
+AC_SUBST(F77_UNDERSCORE)
+AC_SUBST(F77_EXTRA_UNDERSCORE)
+AC_SUBST(LISPEXEC)
+AC_SUBST(LISPEVAL)
+AC_SUBST(srcdir)
+AC_SUBST(libdir)
+AC_SUBST(FLIBS_LIBS)
+AC_SUBST(FLIBS_OPTS)
+AC_SUBST(share_ext)
+AC_SUBST(GIT_VERSION)
+AC_SUBST(HAVE_QL)
+AC_SUBST(F2C)
+
+echo host = $host
+
+# Set the extension for shared libraries. This is not very robust.
+case $host in
+ *darwin*) share_ext=dylib ;;
+ *) share_ext=so ;;
+esac
+
+AC_CONFIG_FILES([
+ matlisp.mk
+ Makefile
+ start.lisp
+ config.lisp
+ lib/lazy-loader.lisp
+ LAPACK/SRC/Makefile
+ LAPACK/BLAS/SRC/Makefile
+ dfftpack/Makefile
+ lib-src/toms715/Makefile
+ lib-src/odepack/Makefile
+ src/f77-mangling.lisp
+])
+
echo FLIBS = $FLIBS
echo HAVE_QL = $HAVE_QL
AC_OUTPUT
commit 7f45afca17c2501c0b2a777809e1998683e5a17f
Author: Raymond Toy <toy...@gm...>
Date: Sun Mar 11 10:48:29 2012 -0700
Don't use full paths for the library on Darwin (OSX).
diff --git a/lib/lazy-loader.lisp.in b/lib/lazy-loader.lisp.in
index d160594..7df2449 100644
--- a/lib/lazy-loader.lisp.in
+++ b/lib/lazy-loader.lisp.in
@@ -115,14 +115,26 @@
;; Tell cffi where our libraries are.
(push "@libdir@/" cffi:*foreign-library-directories*)
-;; Define our libraries
+;; Define our libraries.
+
+;; For some reason, on Darwin (OSX), we can't load the libraries if we
+;; specify the full path. Loading the library fails because there are
+;; undefined symbols (__gfortran_stop_numeric) referenced from the
+;; library. However, everything works if the full path is not given.
+;; This could be a bug in how automake generates libraries on Darwin.
+;; Not really sure.
+;;
+;; This isn't a problem on linux or sparc.
(cffi:define-foreign-library dfftpack
+ (:darwin "libdfftpack.dylib")
(t (:default "@libdir@/libdfftpack")))
(cffi:define-foreign-library toms715
+ (:darwin "libtoms715.dylib")
(t (:default "@libdir@/libtoms715")))
(cffi:define-foreign-library odepack
+ (:darwin "libodepack.dylib")
(t (:default "@libdir@/libodepack")))
(if @ATLAS_P@
@@ -143,8 +155,10 @@
(progn
;; Use our blas and lapack libraries
(cffi:define-foreign-library blas
+ (:darwin "libblas.dylib")
(t (:default "@libdir@/libblas")))
(cffi:define-foreign-library lapack
+ (:darwin "liblapack.dylib")
(t (:default "@libdir@/liblapack")))))
commit 8708b5e4f149f2ba16b879bfa87aa3085097cb34
Author: Raymond Toy <toy...@gm...>
Date: Sun Mar 11 10:42:33 2012 -0700
Can't pass arrays of complex doubles to blas routines in general. Use
the equivalent arrays of doubles.
diff --git a/tests/blas.lisp b/tests/blas.lisp
index 5995f24..b584f72 100644
--- a/tests/blas.lisp
+++ b/tests/blas.lisp
@@ -14,22 +14,22 @@
(list max-error actual expected))))
(rt:deftest blas.zdotu.1.1
- (let ((x (make-array 3 :element-type '(complex double-float) :initial-contents '(#c(1d0 0) #c(2d0 0) #c(3d0 0)))))
+ (let ((x (matlisp::store [#c(1d0 0) #c(2d0 0) #c(3d0 0)])))
(blas:zdotu 2 x 1 x 1))
#c(5d0 0d0))
(rt:deftest blas.zdotu.1.2
- (let ((x (make-array 3 :element-type '(complex double-float) :initial-contents '(#c(1d0 1d0) #c(2d0 2d0) #c(3d0 3d0)))))
+ (let ((x (matlisp::store [#c(1d0 1d0) #c(2d0 2d0) #c(3d0 3d0)])))
(blas:zdotu 2 x 1 x 1))
#c(0d0 10d0))
(rt:deftest blas.zdotc.1.1
- (let ((x (make-array 3 :element-type '(complex double-float) :initial-contents '(#c(1d0 0) #c(2d0 0) #c(3d0 0)))))
+ (let ((x (matlisp::store [#c(1d0 0) #c(2d0 0) #c(3d0 0)])))
(blas:zdotc 2 x 1 x 1))
#c(5d0 0d0))
(rt:deftest blas.zdotc.1.2
- (let ((x (make-array 3 :element-type '(complex double-float) :initial-contents '(#c(1d0 1d0) #c(2d0 2d0) #c(3d0 3d0)))))
+ (let ((x (matlisp::store [#c(1d0 1d0) #c(2d0 2d0) #c(3d0 3d0)])))
(blas:zdotc 2 x 1 x 1))
#c(10d0 0d0))
-----------------------------------------------------------------------
Summary of changes:
configure.ac | 99 +++++++++++++++++++++++++----------------------
lib/lazy-loader.lisp.in | 16 +++++++-
tests/blas.lisp | 8 ++--
3 files changed, 72 insertions(+), 51 deletions(-)
hooks/post-receive
--
matlisp
|
|
From: Raymond T. <rt...@us...> - 2012-03-10 18:51:36
|
This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "matlisp".
The branch, master has been updated
via 4d23b62fbced07e732614e290a08e90a731942ab (commit)
via 579b17087f685d198095ef8616a6374ed59d6848 (commit)
from ec16280b56855cff043ec1fc6bbddc11cf2cebca (commit)
Those revisions listed above that are new to this repository have
not appeared on any other notification email; so we list those
revisions in full, below.
- Log -----------------------------------------------------------------
commit 4d23b62fbced07e732614e290a08e90a731942ab
Merge: 579b170 ec16280
Author: Raymond Toy <toy...@gm...>
Date: Sat Mar 10 10:51:20 2012 -0800
Merge branch 'master' of ssh://matlisp.git.sourceforge.net/gitroot/matlisp/matlisp
commit 579b17087f685d198095ef8616a6374ed59d6848
Author: Raymond Toy <toy...@gm...>
Date: Sat Mar 10 10:48:02 2012 -0800
Fix f77 name mangling and have the ATLAS test program print out the
actual results on failure.
f77_name__ doesn't mean append two underscores. It maeans append one,
and internal underscores are doubled.
diff --git a/configure.ac b/configure.ac
index a5e5e1e..41a9f22 100644
--- a/configure.ac
+++ b/configure.ac
@@ -365,13 +365,13 @@ if test x"$atlas" = xtrue; then
F77*) case $f77_name in
F77_NAME) ZDOTU="ZDOTU" ;;
F77_NAME_) ZDOTU="ZDOTU_" ;;
- F77_NAME__) ZDOTU="ZDOTU__" ;;
+ F77_NAME__) ZDOTU="ZDOTU_" ;;
esac
;;
f77*) case $f77_name in
f77_name) ZDOTU="zdotu" ;;
f77_name_) ZDOTU="zdotu_" ;;
- f77_name__) ZDOTU="zdotu__" ;;
+ f77_name__) ZDOTU="zdotu_" ;;
esac
;;
esac
@@ -403,6 +403,7 @@ int main()
if (fabs(out[0] - 5) < 1e-10) {
rc = 0;
} else {
+ printf("Actual output = %lg, %lg, instead of 5, 0\n", out[0], out[1]);
rc = 1;
}
-----------------------------------------------------------------------
Summary of changes:
configure.ac | 5 +++--
1 files changed, 3 insertions(+), 2 deletions(-)
hooks/post-receive
--
matlisp
|
|
From: Akshay S. <ak...@us...> - 2012-03-10 17:36:21
|
This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "matlisp".
The branch, master has been updated
via ec16280b56855cff043ec1fc6bbddc11cf2cebca (commit)
from 6695b66525322e6817de899b467fb3505ec22775 (commit)
Those revisions listed above that are new to this repository have
not appeared on any other notification email; so we list those
revisions in full, below.
- Log -----------------------------------------------------------------
commit ec16280b56855cff043ec1fc6bbddc11cf2cebca
Author: Akshay Srinivasan <aks...@gm...>
Date: Sat Mar 10 23:01:08 2012 +0530
Make sure lazy-loader.lisp, loads matlisp's version of libraries when asked to.
diff --git a/lib/lazy-loader.lisp.in b/lib/lazy-loader.lisp.in
index f67b5d0..d160594 100644
--- a/lib/lazy-loader.lisp.in
+++ b/lib/lazy-loader.lisp.in
@@ -117,33 +117,36 @@
;; Define our libraries
(cffi:define-foreign-library dfftpack
- (t (:default "libdfftpack")))
+ (t (:default "@libdir@/libdfftpack")))
(cffi:define-foreign-library toms715
- (t (:default "libtoms715")))
+ (t (:default "@libdir@/libtoms715")))
(cffi:define-foreign-library odepack
- (t (:default "libodepack")))
+ (t (:default "@libdir@/libodepack")))
+
+(if @ATLAS_P@
+ (progn
+ ;; Define the ATLAS libraries and their locations.
+ (push "@ATLAS_DIR@" cffi:*foreign-library-directories*)
+ (cffi:define-foreign-library atlas
+ (t (:default "libatlas")))
+
+ (cffi:define-foreign-library cblas
+ (t (:default "libcblas")))
+
+ (cffi:define-foreign-library f77blas
+ (t (:default "libf77blas")))
+
+ (cffi:define-foreign-library lapack
+ (t (:default "liblapack"))))
+ (progn
+ ;; Use our blas and lapack libraries
+ (cffi:define-foreign-library blas
+ (t (:default "@libdir@/libblas")))
+ (cffi:define-foreign-library lapack
+ (t (:default "@libdir@/liblapack")))))
-(when @ATLAS_P@
- ;; Define the ATLAS libraries and their locations.
- (push "@ATLAS_DIR@" cffi:*foreign-library-directories*)
- (cffi:define-foreign-library atlas
- (t (:default "libatlas")))
-
- (cffi:define-foreign-library cblas
- (t (:default "libcblas")))
-
- (cffi:define-foreign-library f77blas
- (t (:default "libf77blas"))))
-
-(unless @ATLAS_P@
- ;; Use our blas and lapack libraries
- (cffi:define-foreign-library blas
- (t (:default "@libdir@/libblas"))))
-
-(cffi:define-foreign-library lapack
- (t (:default "liblapack")))
(defun load-blas-&-lapack-libraries ()
;; Load the additional matlisp libraries
-----------------------------------------------------------------------
Summary of changes:
lib/lazy-loader.lisp.in | 39 +++++++++++++++++++++------------------
1 files changed, 21 insertions(+), 18 deletions(-)
hooks/post-receive
--
matlisp
|
|
From: Raymond T. <rt...@us...> - 2012-03-10 17:24:58
|
This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "matlisp".
The branch, master has been updated
via 6695b66525322e6817de899b467fb3505ec22775 (commit)
via b4da181a45007d71d24371851089844b257e2fa8 (commit)
from b4755e861709360bf208f0b9b479c388003274de (commit)
Those revisions listed above that are new to this repository have
not appeared on any other notification email; so we list those
revisions in full, below.
- Log -----------------------------------------------------------------
commit 6695b66525322e6817de899b467fb3505ec22775
Author: Raymond Toy <toy...@gm...>
Date: Sat Mar 10 09:19:04 2012 -0800
Regenerated.
diff --git a/configure b/configure
index 0ea9407..0f6e275 100755
--- a/configure
+++ b/configure
@@ -15332,6 +15332,76 @@ else
fi
+# Check to see if the ATLAS libraries are compatible with matlisp's
+# ffi. Basically the same test as above that checks to see if -ff2c
+# is needed.
+
+if test x"$atlas" = xtrue; then
+ { $as_echo "$as_me:${as_lineno-$LINENO}: checking if ATLAS is compatible with matlisp's FFI" >&5
+$as_echo_n "checking if ATLAS is compatible with matlisp's FFI... " >&6; }
+ # From the value of f77_name, figure out the actual name for
+ # Fortran's zdotu.
+ case $f77_name in
+ F77*) case $f77_name in
+ F77_NAME) ZDOTU="ZDOTU" ;;
+ F77_NAME_) ZDOTU="ZDOTU_" ;;
+ F77_NAME__) ZDOTU="ZDOTU__" ;;
+ esac
+ ;;
+ f77*) case $f77_name in
+ f77_name) ZDOTU="zdotu" ;;
+ f77_name_) ZDOTU="zdotu_" ;;
+ f77_name__) ZDOTU="zdotu__" ;;
+ esac
+ ;;
+ esac
+
+ cat > conftest.c <<EOF
+
+#include <stdio.h>
+#include <math.h>
+
+extern void ${ZDOTU}(double *, int *, double *, int *, double *, int *);
+
+int main()
+{
+ int rc;
+
+ int n = 2;
+ int incx = 1;
+ double x[4];
+ double out[2];
+ x[0] = 1;
+ x[1] = 0;
+ x[2] = 2;
+ x[3] = 0;
+ out[0] = 0;
+ out[1] = 0;
+
+ ${ZDOTU}(out, &n, x, &incx, x, &incx);
+
+ if (fabs(out[0] - 5) < 1e-10) {
+ rc = 0;
+ } else {
+ rc = 1;
+ }
+
+ return rc;
+}
+
+EOF
+ $CC $CFLAGS -c conftest.c
+ $F77 $FFLAGS -o a.out conftest.o -L${ATLAS_DIR} -latlas -lcblas -lf77blas -llapack
+ if a.out; then
+ { $as_echo "$as_me:${as_lineno-$LINENO}: result: yes" >&5
+$as_echo "yes" >&6; }
+ else
+ { $as_echo "$as_me:${as_lineno-$LINENO}: result: no" >&5
+$as_echo "no" >&6; }
+ as_fn_error $? "ATLAS libraries are not compatible with matlisp." "$LINENO" 5
+ fi
+fi
+
echo FLIBS = $FLIBS
echo HAVE_QL = $HAVE_QL
cat >confcache <<\_ACEOF
commit b4da181a45007d71d24371851089844b257e2fa8
Author: Raymond Toy <toy...@gm...>
Date: Sat Mar 10 09:18:27 2012 -0800
Add check to see if ATLAS libraries are compatible with matlisp.
diff --git a/configure.ac b/configure.ac
index d501cb2..a5e5e1e 100644
--- a/configure.ac
+++ b/configure.ac
@@ -353,6 +353,73 @@ AC_HELP_STRING([--with-atlas=libpath], [Location of the ATLAS libraries]),
])
AM_CONDITIONAL([ATLAS], [test x$atlas = xtrue])
+# Check to see if the ATLAS libraries are compatible with matlisp's
+# ffi. Basically the same test as above that checks to see if -ff2c
+# is needed.
+
+if test x"$atlas" = xtrue; then
+ AC_MSG_CHECKING([if ATLAS is compatible with matlisp's FFI])
+ # From the value of f77_name, figure out the actual name for
+ # Fortran's zdotu.
+ case $f77_name in
+ F77*) case $f77_name in
+ F77_NAME) ZDOTU="ZDOTU" ;;
+ F77_NAME_) ZDOTU="ZDOTU_" ;;
+ F77_NAME__) ZDOTU="ZDOTU__" ;;
+ esac
+ ;;
+ f77*) case $f77_name in
+ f77_name) ZDOTU="zdotu" ;;
+ f77_name_) ZDOTU="zdotu_" ;;
+ f77_name__) ZDOTU="zdotu__" ;;
+ esac
+ ;;
+ esac
+
+ cat > conftest.c <<EOF
+[
+#include <stdio.h>
+#include <math.h>
+
+extern void ${ZDOTU}(double *, int *, double *, int *, double *, int *);
+
+int main()
+{
+ int rc;
+
+ int n = 2;
+ int incx = 1;
+ double x[4];
+ double out[2];
+ x[0] = 1;
+ x[1] = 0;
+ x[2] = 2;
+ x[3] = 0;
+ out[0] = 0;
+ out[1] = 0;
+
+ ${ZDOTU}(out, &n, x, &incx, x, &incx);
+
+ if (fabs(out[0] - 5) < 1e-10) {
+ rc = 0;
+ } else {
+ rc = 1;
+ }
+
+ return rc;
+}
+]
+EOF
+ $CC $CFLAGS -c conftest.c
+ $F77 $FFLAGS -o a.out conftest.o -L${ATLAS_DIR} -latlas -lcblas -lf77blas -llapack
+ if a.out; then
+ AC_MSG_RESULT([yes])
+ else
+ AC_MSG_RESULT([no])
+ AC_MSG_ERROR([ATLAS libraries are not compatible with matlisp.])
+ fi
+fi
+
echo FLIBS = $FLIBS
echo HAVE_QL = $HAVE_QL
AC_OUTPUT
-----------------------------------------------------------------------
Summary of changes:
configure | 70 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
configure.ac | 67 +++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 files changed, 137 insertions(+), 0 deletions(-)
hooks/post-receive
--
matlisp
|
|
From: Akshay S. <ak...@us...> - 2012-03-10 15:55:32
|
This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "matlisp".
The branch, matlisp-cffi has been updated
via 96f2abfd9395d6f25520c3b828c4c69a0f35d8a3 (commit)
via b4755e861709360bf208f0b9b479c388003274de (commit)
from e6c42164c51c0fbef18cb02df5c5d8f37ce16ef5 (commit)
Those revisions listed above that are new to this repository have
not appeared on any other notification email; so we list those
revisions in full, below.
- Log -----------------------------------------------------------------
commit 96f2abfd9395d6f25520c3b828c4c69a0f35d8a3
Merge: e6c4216 b4755e8
Author: Akshay Srinivasan <aks...@gm...>
Date: Sat Mar 10 21:22:14 2012 +0530
Merge branch 'master' into matlisp-cffi
-----------------------------------------------------------------------
Summary of changes:
lib/lazy-loader.lisp.in | 2 +-
1 files changed, 1 insertions(+), 1 deletions(-)
hooks/post-receive
--
matlisp
|
|
From: Akshay S. <ak...@us...> - 2012-03-10 15:53:38
|
This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "matlisp".
The branch, master has been updated
via b4755e861709360bf208f0b9b479c388003274de (commit)
from 035339c307401f593732418ecdba8fa6785b678d (commit)
Those revisions listed above that are new to this repository have
not appeared on any other notification email; so we list those
revisions in full, below.
- Log -----------------------------------------------------------------
commit b4755e861709360bf208f0b9b479c388003274de
Author: Akshay Srinivasan <aks...@gm...>
Date: Sat Mar 10 21:13:48 2012 +0530
Made sure that CFFI loads our version of BLAS, when it is asked to.
diff --git a/lib/lazy-loader.lisp.in b/lib/lazy-loader.lisp.in
index 11f3358..f67b5d0 100644
--- a/lib/lazy-loader.lisp.in
+++ b/lib/lazy-loader.lisp.in
@@ -140,7 +140,7 @@
(unless @ATLAS_P@
;; Use our blas and lapack libraries
(cffi:define-foreign-library blas
- (t (:default "libblas"))))
+ (t (:default "@libdir@/libblas"))))
(cffi:define-foreign-library lapack
(t (:default "liblapack")))
-----------------------------------------------------------------------
Summary of changes:
lib/lazy-loader.lisp.in | 2 +-
1 files changed, 1 insertions(+), 1 deletions(-)
hooks/post-receive
--
matlisp
|
|
From: Akshay S. <ak...@us...> - 2012-03-09 16:23:44
|
This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "matlisp".
The branch, matlisp-cffi has been updated
via e6c42164c51c0fbef18cb02df5c5d8f37ce16ef5 (commit)
via 035339c307401f593732418ecdba8fa6785b678d (commit)
via 068408a227df4427520decc0f10fe0a253a5cdbb (commit)
via 1e8efae3363f6d50e91831cf39da6b64c106cea3 (commit)
via 91cfba559ec57243acb0d09ae0a9840fb7e4082b (commit)
via c889997acda2f13b2d61c442a61f5e3b95004016 (commit)
via c9ac9850a4916f2871fa49efd6beb1d5d38b357b (commit)
via 089a044899e6f499821fbdefa69157e0c329c0cf (commit)
from 228fb9ad5bc2e26dccbfd4f5ad2d2228bab8c39e (commit)
Those revisions listed above that are new to this repository have
not appeared on any other notification email; so we list those
revisions in full, below.
- Log -----------------------------------------------------------------
commit e6c42164c51c0fbef18cb02df5c5d8f37ce16ef5
Merge: 228fb9a 035339c
Author: Akshay Srinivasan <aks...@gm...>
Date: Fri Mar 9 21:50:22 2012 +0530
Merge branch 'master' into matlisp-cffi
Conflicts:
matlisp.asd
src/gemm.lisp
diff --cc matlisp.asd
index 6288219,dc270d0..4632cf4
--- a/matlisp.asd
+++ b/matlisp.asd
@@@ -304,4 -297,20 +309,20 @@@
:components
((:module "src"
:components
- ((:file "dlsode")))))
+ ((:file "dlsode")))))
+
+ (defmethod perform ((op asdf:test-op) (c (eql (asdf:find-system :matlisp))))
+ (oos 'asdf:test-op 'matlisp-tests))
+
+ (asdf:defsystem matlisp-tests
+ :depends-on (matlisp)
+ :in-order-to ((compile-op (load-op :rt))
+ (test-op (load-op :rt :matlisp)))
+ :components
+ ((:module "tests"
+ :components
+ ((:file "blas")))))
+
+ (defmethod perform ((op test-op) (c (eql (asdf:find-system :matlisp-tests))))
+ (or (funcall (intern "DO-TESTS" (find-package "RT")))
- (error "TEST-OP failed for MATLISP-TESTS")))
++ (error "TEST-OP failed for MATLISP-TESTS")))
diff --cc src/ffi-cffi-interpreter-specific.lisp
index 06b7fdb,db2db4c..75afa9a
--- a/src/ffi-cffi-interpreter-specific.lisp
+++ b/src/ffi-cffi-interpreter-specific.lisp
@@@ -93,17 -99,9 +99,17 @@@ Return
;; the alien object but before the alien function is called. Let's
;; be safe rather than sorry.
`(with-fortran-float-modes
- (without-gcing
- (let (,@(mapcar #'(lambda (pair)
- `(,(first pair)
- (vector-data-address ,(second pair))))
- vlist))
- ,@body))))
+ (without-gcing
+ (let (,@(mapcar #'(lambda (lst)
+ (destructuring-bind (addr-var var &key inc-type inc) lst
+ `(,addr-var ,@(if inc
+ `((cffi:inc-pointer (vector-data-address ,var)
+ ,@(case inc-type
+ (:double-float `((* ,inc 8)))
+ (:single-float `((* ,inc 4)))
+ (:complex-double-float `((* ,inc 16)))
+ (:complex-single-float `((* ,inc 8)))
+ (t `(,inc)))))
+ `((vector-data-address ,var))))))
+ vlist))
- ,@body))))
++ ,@body))))
-----------------------------------------------------------------------
Summary of changes:
configure | 1 +
configure.ac | 1 +
lib/lazy-loader.lisp.in | 2 +-
matlisp.asd | 23 ++++++-
src/axpy.lisp | 10 ++--
src/copy.lisp | 6 +-
src/diag.lisp | 6 +-
src/f77-mangling.lisp.in | 1 +
src/ffi-cffi-interpreter-specific.lisp | 10 ++-
src/gemm.lisp | 2 +-
src/mplus.lisp | 34 +++++-----
src/ref.lisp | 125 +++++++++++++++-----------------
src/scal.lisp | 20 +++---
src/specfun.lisp | 4 +-
tests/blas.lisp | 39 ++++++++++
15 files changed, 173 insertions(+), 111 deletions(-)
create mode 100644 tests/blas.lisp
hooks/post-receive
--
matlisp
|
|
From: Raymond T. <rt...@us...> - 2012-03-09 06:24:26
|
This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "matlisp".
The branch, master has been updated
via 035339c307401f593732418ecdba8fa6785b678d (commit)
from 068408a227df4427520decc0f10fe0a253a5cdbb (commit)
Those revisions listed above that are new to this repository have
not appeared on any other notification email; so we list those
revisions in full, below.
- Log -----------------------------------------------------------------
commit 035339c307401f593732418ecdba8fa6785b678d
Author: Raymond Toy <toy...@gm...>
Date: Thu Mar 8 22:24:08 2012 -0800
Add support for test-op. Requires RT.
Usage: (asdf:oos 'asdf:test-op :matlisp)
matlisp.asd:
o Add defsystem for running the tests and support infrastructure to
compile and run the tests.
tests/blas.lisp:
o Some simple tests for blas functions.
diff --git a/matlisp.asd b/matlisp.asd
index d9b4389..dc270d0 100644
--- a/matlisp.asd
+++ b/matlisp.asd
@@ -28,6 +28,11 @@
(in-package #:common-lisp-user)
+(defpackage #:matlisp-system
+ (:use #:cl :asdf))
+
+(in-package #:matlisp-system)
+
(asdf:defsystem matlisp-packages
:pathname #.(translate-logical-pathname "matlisp:srcdir;")
:components
@@ -294,3 +299,18 @@
:components
((:file "dlsode")))))
+(defmethod perform ((op asdf:test-op) (c (eql (asdf:find-system :matlisp))))
+ (oos 'asdf:test-op 'matlisp-tests))
+
+(asdf:defsystem matlisp-tests
+ :depends-on (matlisp)
+ :in-order-to ((compile-op (load-op :rt))
+ (test-op (load-op :rt :matlisp)))
+ :components
+ ((:module "tests"
+ :components
+ ((:file "blas")))))
+
+(defmethod perform ((op test-op) (c (eql (asdf:find-system :matlisp-tests))))
+ (or (funcall (intern "DO-TESTS" (find-package "RT")))
+ (error "TEST-OP failed for MATLISP-TESTS")))
diff --git a/tests/blas.lisp b/tests/blas.lisp
new file mode 100644
index 0000000..5995f24
--- /dev/null
+++ b/tests/blas.lisp
@@ -0,0 +1,39 @@
+(in-package #:matlisp-user)
+
+(asdf:oos 'asdf:load-op :rt)
+
+(defmethod max-matrix-diff ((actual standard-matrix) (expected standard-matrix) &key (allowed-error 0d0))
+ (let ((max-error (reduce #'max
+ (map 'list
+ #'(lambda (x y)
+ (abs (- x y)))
+ (matlisp::store actual)
+ (matlisp::store expected)))))
+ (if (<= max-error allowed-error)
+ t
+ (list max-error actual expected))))
+
+(rt:deftest blas.zdotu.1.1
+ (let ((x (make-array 3 :element-type '(complex double-float) :initial-contents '(#c(1d0 0) #c(2d0 0) #c(3d0 0)))))
+ (blas:zdotu 2 x 1 x 1))
+ #c(5d0 0d0))
+
+(rt:deftest blas.zdotu.1.2
+ (let ((x (make-array 3 :element-type '(complex double-float) :initial-contents '(#c(1d0 1d0) #c(2d0 2d0) #c(3d0 3d0)))))
+ (blas:zdotu 2 x 1 x 1))
+ #c(0d0 10d0))
+
+(rt:deftest blas.zdotc.1.1
+ (let ((x (make-array 3 :element-type '(complex double-float) :initial-contents '(#c(1d0 0) #c(2d0 0) #c(3d0 0)))))
+ (blas:zdotc 2 x 1 x 1))
+ #c(5d0 0d0))
+
+(rt:deftest blas.zdotc.1.2
+ (let ((x (make-array 3 :element-type '(complex double-float) :initial-contents '(#c(1d0 1d0) #c(2d0 2d0) #c(3d0 3d0)))))
+ (blas:zdotc 2 x 1 x 1))
+ #c(10d0 0d0))
+
+(rt:deftest blas.axpy.1
+ (let ((x [1 2 3]))
+ (max-matrix-diff (axpy 2 x x) [3 6 9]))
+ t)
-----------------------------------------------------------------------
Summary of changes:
matlisp.asd | 20 ++++++++++++++++++++
tests/blas.lisp | 39 +++++++++++++++++++++++++++++++++++++++
2 files changed, 59 insertions(+), 0 deletions(-)
create mode 100644 tests/blas.lisp
hooks/post-receive
--
matlisp
|
|
From: Raymond T. <rt...@us...> - 2012-03-09 04:37:50
|
This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "matlisp".
The branch, matlisp-cffi has been updated
via 228fb9ad5bc2e26dccbfd4f5ad2d2228bab8c39e (commit)
from 39833b146ce09d18e1909ded8c81e65c250ae0c3 (commit)
Those revisions listed above that are new to this repository have
not appeared on any other notification email; so we list those
revisions in full, below.
- Log -----------------------------------------------------------------
commit 228fb9ad5bc2e26dccbfd4f5ad2d2228bab8c39e
Author: Raymond Toy <toy...@gm...>
Date: Thu Mar 8 20:37:36 2012 -0800
Remove test text.
diff --git a/README b/README
index 35fa251..bcb5210 100644
--- a/README
+++ b/README
@@ -117,7 +117,3 @@ m/
divide two matrices. (m/ a b) means the same as inv(B)*A.
(m/ a) is the same as inv(A).
---------------------
-Testing for commit emails.
-Test 2
-Test 3
-----------------------------------------------------------------------
Summary of changes:
README | 4 ----
1 files changed, 0 insertions(+), 4 deletions(-)
hooks/post-receive
--
matlisp
|
|
From: Raymond T. <toy...@gm...> - 2012-03-09 03:53:11
|
direct testing. |
|
From: Raymond T. <rt...@us...> - 2012-03-09 03:11:31
|
This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "matlisp".
The branch, matlisp-cffi has been updated
via 681f4ecf69ef92a8b032591c2fd17084cc037688 (commit)
from f588a09ffa90a7263dc14cd5269a8d131babec7f (commit)
Those revisions listed above that are new to this repository have
not appeared on any other notification email; so we list those
revisions in full, below.
- Log -----------------------------------------------------------------
commit 681f4ecf69ef92a8b032591c2fd17084cc037688
Author: Raymond Toy <toy...@gm...>
Date: Thu Mar 8 19:11:15 2012 -0800
Test 1 for commit messages.
diff --git a/README b/README
index 3c397d7..f7dacb3 100644
--- a/README
+++ b/README
@@ -115,4 +115,7 @@ m*
m/
divide two matrices. (m/ a b) means the same as inv(B)*A.
- (m/ a) is the same as inv(A).
\ No newline at end of file
+ (m/ a) is the same as inv(A).
+
+--------------------
+Testing for commit emails.
-----------------------------------------------------------------------
Summary of changes:
README | 5 ++++-
1 files changed, 4 insertions(+), 1 deletions(-)
hooks/post-receive
--
matlisp
|
