[lisplab-cvs] r26 - src/fft src/linalg src/matlisp src/matrix system
Jørn Inge Vestgården
jivestgarden at common-lisp.net
Sun May 17 19:02:05 UTC 2009
Author: jivestgarden
Date: Sun May 17 15:02:03 2009
New Revision: 26
Log:
refactoring almost complete. Needs tidying
Modified:
TODO
src/fft/level3-fft-blas.lisp
src/linalg/level3-linalg-blas-real.lisp
src/linalg/level3-linalg-generic.lisp
src/matlisp/geev.lisp
src/matlisp/inv.lisp
src/matlisp/mul.lisp
src/matrix/level1-classes.lisp
src/matrix/level1-constructors.lisp
src/matrix/level1-generic.lisp
src/matrix/level1-interface.lisp
src/matrix/level1-matrix.lisp
src/matrix/level1-util.lisp
src/matrix/level2-generic.lisp
src/matrix/level2-matrix-dge.lisp
src/matrix/level2-matrix-zge.lisp
system/lisplab.asd
Modified: TODO
==============================================================================
--- TODO (original)
+++ TODO Sun May 17 15:02:03 2009
@@ -1,4 +1,5 @@
TODO
+o Some way to just extend the exact input matrix type
o Documentation.
o Check if there is some non-threadsafe code in the fortran interface,
i.e. array pre-alocation for the workspaces.
Modified: src/fft/level3-fft-blas.lisp
==============================================================================
--- src/fft/level3-fft-blas.lisp (original)
+++ src/fft/level3-fft-blas.lisp Sun May 17 15:02:03 2009
@@ -19,41 +19,59 @@
;;; TODO should use the normal ref-blas-complex-store
+;;; TODO fix the methods so that they use the actual input matrix type, not just
+;;; the eql spezializer.
+
(in-package :lisplab)
-(defmethod fft1 ((x blas-real))
- (fft1! (convert x 'blas-complex)))
+;;;; Real matrices
+
+(defmethod fft1 ((x matrix-lisp-dge))
+ (fft1! (convert x 'matrix-zge)))
+
+(defmethod ifft1 ((x matrix-lisp-dge))
+ (ifft1! (convert x 'matrix-zge)))
+
+(defmethod ifft2 ((x matrix-lisp-dge))
+ (ifft2! (convert x 'matrix-zge)))
+
+(defmethod fft2 ((x matrix-lisp-dge))
+ (fft2! (convert x 'matrix-zge)))
+
+;;; Complex matrices
-(defmethod ifft1 ((x blas-real))
- (ifft1! (convert x 'blas-complex)))
+(defmethod fft1 ((x matrix-lisp-zge))
+ (fft1! (copy x)))
-(defmethod ifft2 ((x blas-real))
- (ifft2! (convert x 'blas-complex)))
+(defmethod ifft1 ((x matrix-lisp-zge))
+ (ifft1! (copy x)))
-(defmethod fft2 ((x blas-real))
- (fft2! (convert x 'blas-complex)))
+(defmethod ifft2 ((x matrix-lisp-zge))
+ (ifft2! (copy x)))
+(defmethod fft2 ((x matrix-lisp-zge))
+ (fft2! (copy x)))
-(defmethod fft1! ((x blas-complex))
+(defmethod fft1! ((x matrix-lisp-zge))
(dotimes (i (cols x))
- (fft-radix-2-blas-complex-store! :f (store x) (rows x) (* (rows x) i) 1))
+ (fft-radix-2-blas-complex-store! :f (matrix-store x) (rows x) (* (rows x) i) 1))
x)
-(defmethod ifft1! ((x blas-complex))
+(defmethod ifft1! ((x matrix-lisp-zge))
(dotimes (i (cols x))
- (fft-radix-2-blas-complex-store! :r (store x) (rows x) (* (rows x) i) 1))
+ (fft-radix-2-blas-complex-store! :r (matrix-store x) (rows x) (* (rows x) i) 1))
x)
-(defmethod fft2! ((x blas-complex))
+(defmethod fft2! ((x matrix-lisp-zge))
(fft1! x)
(dotimes (i (rows x))
- (fft-radix-2-blas-complex-store! :f (store x) (cols x) i (rows x)))
+ (fft-radix-2-blas-complex-store! :f (matrix-store x) (cols x) i (rows x)))
x)
-(defmethod ifft2! ((x blas-complex))
+(defmethod ifft2! ((x matrix-lisp-zge))
(ifft1! x)
(dotimes (i (rows x))
- (fft-radix-2-blas-complex-store! :r (store x) (cols x) i (rows x)))
+ (fft-radix-2-blas-complex-store! :r (matrix-store x) (cols x) i (rows x)))
x)
(declaim (ftype (function
Modified: src/linalg/level3-linalg-blas-real.lisp
==============================================================================
--- src/linalg/level3-linalg-blas-real.lisp (original)
+++ src/linalg/level3-linalg-blas-real.lisp Sun May 17 15:02:03 2009
@@ -19,14 +19,6 @@
(in-package :lisplab)
-;; TODO move these optimized functions to library
-
-(defmacro *df (a b) `(truly-the double-float (* ,a ,b)))
-(defmacro /df (a b) `(truly-the double-float (/ ,a ,b)))
-(defmacro +df (a b) `(truly-the double-float (+ ,a ,b)))
-(defmacro -df (a b) `(truly-the double-float (- ,a ,b)))
-
-
#+todo (defmethod mtr (matrix)
(let ((ans 0))
(dotimes (i (rows matrix))
@@ -34,26 +26,26 @@
ans))
#+todo (defmethod mtp (a)
- (let ((b (create a 0 (list (cols a) (rows a)))))
+ (let ((b (mcreate a 0 (list (cols a) (rows a)))))
(dotimes (i (rows b))
(dotimes (j (cols b))
(setf (mref b i j) (mref a j i))))
b))
-(defmethod mconj ((a blas-real))
+(defmethod mconj ((a matrix-lisp-dge))
(copy a))
-(defmethod mct ((a blas-real))
+(defmethod mct ((a matrix-lisp-dge))
(mtp a))
-(defmethod m* ((a blas-real) (b blas-real))
+(defmethod m* ((a matrix-lisp-dge) (b matrix-lisp-dge))
(let* ((N (rows a))
(M (cols b))
(S (rows b))
- (c (create a 0 (list N M)))
- (a2 (store a))
- (b2 (store b))
- (c2 (store c)))
+ (c (mcreate a 0 (list N M)))
+ (a2 (matrix-store a))
+ (b2 (matrix-store b))
+ (c2 (matrix-store c)))
(declare (type-blas-store a2 b2 c2)
(type-blas-idx N M S))
(macrolet ((refa (i j) `(ref-blas-real-store A2 ,i ,j N))
@@ -68,7 +60,7 @@
(setf (refc i j) cij))))
c)))
-(defmethod LU-factor! ((A blas-real) p)
+(defmethod LU-factor! ((A matrix-lisp-dge) p)
;; Translation from GSL.
;; Destructive LU factorization. The outout is PA=LU,
;; stored in one matrix, where the diagonal elements belong
@@ -76,7 +68,7 @@
;; Assumes the permutation vector to be initilized
(let ((N (rows A))
(sign 1)
- (A2 (store A)))
+ (A2 (matrix-store A)))
(declare (type-blas-idx N)
(fixnum sign)
(type-blas-store a2)
@@ -111,9 +103,9 @@
(defun L-solve!-blas-real (L x col)
;; Solve Lx=b
- (declare (blas-real L x))
- (let ((L2 (store L))
- (x2 (store x))
+ (declare (matrix-lisp-dge L x))
+ (let ((L2 (matrix-store L))
+ (x2 (matrix-store x))
(N (rows x)))
(declare (type-blas-store L2 x2)
(type-blas-idx N col))
@@ -128,9 +120,9 @@
x)
(defun U-solve!-blas-real (U x col)
- (declare (blas-real U x))
- (let* ((U2 (store U))
- (x2 (store x))
+ (declare (matrix-lisp-dge U x))
+ (let* ((U2 (matrix-store U))
+ (x2 (matrix-store x))
(N (rows x))
(N-1 (1- N))
(N-2 (1- N-1)))
@@ -152,7 +144,7 @@
(U-solve!-blas-real LU x col)
x)
-(defmethod lin-solve ((A blas-real) (b blas-real))
+(defmethod lin-solve ((A matrix-lisp-dge) (b matrix-lisp-dge))
(destructuring-bind (LU pvec sign) (LU-factor A)
(let ((b2 (copy b)))
(dotimes (i (rows A))
@@ -170,10 +162,10 @@
(LU-solve!-blas-real LU A (vref p i)))))
A)
-(defmethod minv! ((A blas-real))
+(defmethod minv! ((A matrix-lisp-dge))
(minv!-blas-real A))
-(defmethod minv ((A blas-real))
+(defmethod minv ((A matrix-lisp-dge))
(minv! (copy A)))
Modified: src/linalg/level3-linalg-generic.lisp
==============================================================================
--- src/linalg/level3-linalg-generic.lisp (original)
+++ src/linalg/level3-linalg-generic.lisp Sun May 17 15:02:03 2009
@@ -17,6 +17,9 @@
;;; with this program; if not, write to the Free Software Foundation, Inc.,
;;; 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+;;; TODO clean up. Move read and write out
+
(in-package :lisplab)
(export '(pgmwrite))
@@ -28,14 +31,14 @@
ans))
(defmethod mtp (a)
- (let ((b (create a 0 (list (cols a) (rows a)))))
+ (let ((b (mcreate a 0 (list (cols a) (rows a)))))
(dotimes (i (rows b))
(dotimes (j (cols b))
(setf (mref b i j) (mref a j i))))
b))
(defmethod mconj (a)
- (let ((b (create a #C(0 0) (list (rows a) (cols a)) )))
+ (let ((b (mcreate a #C(0 0) (list (rows a) (cols a)) )))
(dotimes (i (size b))
(setf (vref b i) (conjugate (vref a i))))
b))
@@ -44,7 +47,7 @@
(mconj (mtp a)))
(defmethod m* (a b)
- (let ((c (create a 0 (list (rows a) (cols b)))))
+ (let ((c (mcreate a 0 (list (rows a) (cols b)))))
(dotimes (i (rows c))
(dotimes (j (cols c))
(dotimes (k (cols a))
@@ -197,7 +200,7 @@
(make-permutation-vector (rows A)))
(let ((L A)
(U (copy A))
- (Pmat (create A 0)))
+ (Pmat (mcreate A 0)))
(w/mat L (x i j) (cond ((> i j) x) ((= i j) 1) (t 0)))
(w/mat U (x i j) (cond ((<= i j) x) (t 0)))
(dotimes (i (rows P))
@@ -265,82 +268,5 @@
-
-;;; TRASH
-
-#+nil (defun apply-permutation (a p)
- this is inverse
- (let ((b (create a 0)))
- (dotimes (row (rows a))
- (let ((s (vref p row)))
- (dotimes (col (cols a))
- (setf (mref b row col) (mref a s col)))))
- b))
-
-#+nil (defun apply-inverse-permutation (a p)
- this is forward
- (let ((b (create a 0)))
- (dotimes (row (rows a))
- (let ((s (vref p row)))
- (dotimes (col (cols a))
- (setf (mref b s col) (mref a row col)))))
- b))
-
-
-#+nil (defmethod LU-factor! (A p)
- ;; Versjon fra boka!
- (let* ((N (rows A))
- (N-1 (1- N))
- (det 1))
- (dotimes (i N)
- (setf (vref p i) i))
- (dotimes (i N-1)
- (let ((i-pivot i))
- (loop for j from (1+ i) below N do
- (when (> (abs (mref A j i))
- (abs (mref A i-pivot i)))
- (setf i-pivot j)))
- (unless (= i-pivot i)
- (let ((tmp (vref p i)))
- (setf (vref p i) (vref p i-pivot)
- (vref p i-pivot) tmp
- det (- det)))
- (dotimes (j N)
- (let ((tmp (mref A i j)))
- (setf (mref A i j) (mref A i-pivot j)
- (mref A i-pivot j) tmp)))))
- ;; Now reduce all elementsbelow the i'th row
- (unless (zerop (mref A i i))
- (loop for r from (1+ i) below N do
- (print (list 'foerr 'i i 'r r A))
- (setf (mref A r i) (./ (mref A r i) (mref A i i)))
- (loop for c from (1+ i) below N do
- (setf (mref A r c) (.- (mref A r c)
- (.* (mref A r i) (mref A i c))))
- (print (list 'mellom 'r r 'c c (mref A r c))))
- (print (list 'etter 'i i 'r r A)))))
- (list A p det)))
-
-#+nil (defun tmp-LU-mul (A)
- ;; Test code. TODO move and make to an automated test
- (destructuring-bind (LU p det)
- (LU-factor A)
- (let ((L (create LU 0))
- (U (create LU 0)))
- (dotimes (i (rows A))
- (setf (mref L i i) 1)
- (loop for j from 0 below i do
- (setf (mref L i j)
- (mref LU i j))))
- (dotimes (i (rows A))
- #+nil (setf (mref U i i) 1)
- (loop for j from i below (cols A) do
- (setf (mref U i j)
- (mref LU i j))))
- (list A
- (apply-inverse-permutation (m* L U) p)
- LU
- p L U))))
-
Modified: src/matlisp/geev.lisp
==============================================================================
--- src/matlisp/geev.lisp (original)
+++ src/matlisp/geev.lisp Sun May 17 15:02:03 2009
@@ -32,12 +32,12 @@
(in-package :lisplab)
-(defmethod eigenvectors ((a blas-real))
+(defmethod eigenvectors ((a matrix-blas-dge))
(destructuring-bind (evals vl-mat vr-mat)
(dgeev (copy a) nil (create a 0))
(list evals vr-mat)))
-(defmethod eigenvalues ((a blas-real))
+(defmethod eigenvalues ((a matrix-blas-dge))
(destructuring-bind (evals vl-mat vr-mat)
(dgeev (copy a) nil nil)
evals))
@@ -47,7 +47,7 @@
(defmethod rearrange-eigenvector-matrix (v p)
p)
-(defmethod rearrange-eigenvector-matrix ((evals blas-complex) (p blas-real ))
+(defmethod rearrange-eigenvector-matrix ((evals matrix-blas-zge) (p matrix-blas-dge))
(let* ((n (size evals))
(evec (cnew 0 n n)))
(do ((col 0 (incf col)))
@@ -66,11 +66,11 @@
(defun combine-ri-vectors (wr wi)
(let* ((n (size wr))
- (wr2 (make-instance 'blas-real :rows n :cols 1 :size n :store wr))
- (wi2 (make-instance 'blas-real :rows n :cols 1 :size n :store wi)))
+ (wr2 (make-instance 'matrix-dge :rows n :cols 1 :store wr))
+ (wi2 (make-instance 'matrix-dge :rows n :cols 1 :store wi)))
(if (.= wi2 0)
wr2
- (.+ wr2 (.* %i (convert wi2 'blas-complex))))))
+ (.+ wr2 (.* %i (convert wi2 'matrix-zge))))))
(defun dgeev-workspace-size (n lv? rv?)
;; Ask geev how much space it wants for the work array
@@ -112,7 +112,7 @@
(if vl-mat "V" "N") ; JOBVL
(if vr-mat "V" "N") ; JOBVR
n ; N
- (store a) ; A
+ (matrix-store a) ; A
n ; LDA
wr ; WR
wi ; WI
@@ -128,12 +128,12 @@
(vr-mat2 (rearrange-eigenvector-matrix evec vr-mat)))
(list evec vl-mat2 vr-mat2)))))
-(defmethod eigenvectors ((a blas-complex))
+(defmethod eigenvectors ((a matrix-zge))
(destructuring-bind (evals vl-mat vr-mat)
(zgeev (copy a) nil (create a 0))
(list evals vr-mat)))
-(defmethod eigenvalues ((a blas-complex))
+(defmethod eigenvalues ((a matrix-zge))
(destructuring-bind (evals vl-mat vr-mat)
(zgeev (copy a) nil nil)
evals))
@@ -166,8 +166,8 @@
(2n (* 2 n))
(xxx (allocate-real-store 2))
(w (cnew 0 n 1))
- (vl (if vl-mat (store vl-mat) xxx))
- (vr (if vr-mat (store vr-mat) xxx))
+ (vl (if vl-mat (matrix-store vl-mat) xxx))
+ (vr (if vr-mat (matrix-store vr-mat) xxx))
(lwork (zgeev-workspace-size n (if vl-mat t nil) (if vr-mat t nil)))
(work (allocate-real-store lwork))
(rwork (allocate-real-store 2n)))
@@ -176,9 +176,9 @@
(if vl-mat "V" "N") ;; JOBVL
(if vr-mat "V" "N") ;; JOBVR
n ;; N
- (store a) ;; A
+ (matrix-store a) ;; A
n ;; LDA
- (store w) ;; W
+ (matrix-store w) ;; W
vl ;; VL
(if vl-mat n 1) ;; LDVL
vr ;; VR
Modified: src/matlisp/inv.lisp
==============================================================================
--- src/matlisp/inv.lisp (original)
+++ src/matlisp/inv.lisp Sun May 17 15:02:03 2009
@@ -19,44 +19,44 @@
(in-package :lisplab)
-(defmethod minv! ((a blas-real))
+(defmethod minv! ((a matrix-blas-dge))
(let* ((N (rows a))
(ipiv (make-array N :element-type '(unsigned-byte 32)))
(msg "argument A given to minv is singular to working machine precision"))
(multiple-value-bind (_ ipiv info)
- (f77::dgetrf N N (store a) N ipiv 0)
+ (f77::dgetrf N N (matrix-store a) N ipiv 0)
(declare (ignore _))
(unless (zerop info)
(error msg))
(let ((work (make-array N :element-type 'double-float)))
(multiple-value-bind (_ __ info)
- (f77::dgetri N (store a) N ipiv work N 0)
+ (f77::dgetri N (matrix-store a) N ipiv work N 0)
(declare (ignore _ __))
(unless (zerop info)
(error msg))
a)))))
-(defmethod minv ((a blas-real))
+(defmethod minv ((a matrix-blas-dge))
(minv! (copy a)))
-(defmethod minv! ((a blas-complex))
+(defmethod minv! ((a matrix-blas-zge))
(let* ((N (rows a))
(ipiv (make-array N :element-type '(unsigned-byte 32)))
(msg "argument A given to mdiv is singular to working machine precision"))
(multiple-value-bind (_ ipiv info)
- (f77::zgetrf N N (store a) N ipiv 0)
+ (f77::zgetrf N N (matrix-store a) N ipiv 0)
(declare (ignore _))
(unless (zerop info)
(error msg ))
(let ((work (make-array (* 2 N) :element-type 'double-float)))
(multiple-value-bind (_ __ info)
- (f77::zgetri N (store a) N ipiv work N 0)
+ (f77::zgetri N (matrix-store a) N ipiv work N 0)
(declare (ignore _ __))
(unless (zerop info)
(error msg))
a)))))
-(defmethod minv ((a blas-complex))
+(defmethod minv ((a matrix-blas-zge))
(minv! (copy a)))
Modified: src/matlisp/mul.lisp
==============================================================================
--- src/matlisp/mul.lisp (original)
+++ src/matlisp/mul.lisp Sun May 17 15:02:03 2009
@@ -19,18 +19,18 @@
(in-package :lisplab)
-(defmethod m* ((a blas-real) (b blas-real))
+(defmethod m* ((a matrix-blas-dge) (b matrix-blas-dge))
(let* ((m (rows a))
(n (cols b))
(k (cols a))
- (c (new 'blas-real (list m n) t 0.0)))
- (f77::dgemm "N" "N" m n k 1.0 (store a) m (store b) k 0.0 (store c) m)
+ (c (mcreate a 0 (list m n))))
+ (f77::dgemm "N" "N" m n k 1.0 (matrix-store a) m (matrix-store b) k 0.0 (matrix-store c) m)
c))
-(defmethod m* ((a blas-complex) (b blas-complex))
+(defmethod m* ((a matrix-blas-zge) (b matrix-blas-zge))
(let* ((m (rows a))
(n (cols b))
(k (cols a))
- (c (new 'blas-complex (list m n) t 0.0)))
- (f77::zgemm "N" "N" m n k #C(1.0 0.0) (store a) m (store b) k #C(0.0 0.0) (store c) m)
+ (c (mcreate a 0 (list m n))))
+ (f77::zgemm "N" "N" m n k #C(1.0 0.0) (matrix-store a) m (matrix-store b) k #C(0.0 0.0) (matrix-store c) m)
c))
Modified: src/matrix/level1-classes.lisp
==============================================================================
--- src/matrix/level1-classes.lisp (original)
+++ src/matrix/level1-classes.lisp Sun May 17 15:02:03 2009
@@ -100,7 +100,7 @@
(defclass matrix-lisp-dge (matrix-implementation-lisp matrix-base-dge) ())
-(defclass matrix-dge (matrix-implementation-blas matrix-lisp-dge) ())
+(defclass matrix-blas-dge (matrix-implementation-blas matrix-lisp-dge) ())
(defclass matrix-dge (matrix-blas-dge) ()
(:documentation "General matrix with double float elements."))
@@ -144,7 +144,31 @@
(matrix-structure-diagonal matrix-element-complex-double-float matrix-implementation-base)
())
+;;; Function matrices (matrices without a store)
+(defclass function-matrix
+ (matrix-structure-general matrix-element-base matrix-implementation-base)
+ ((mref
+ :initarg :mref
+ :initform (constantly 0)
+ :accessor function-matrix-mref
+ :type function)
+ (set-mref
+ :initarg :set-mref
+ :initform (constantly nil)
+ :accessor function-matrix-set-mref
+ :type function)
+ (vref
+ :initarg :vref
+ :initform (constantly 0)
+ :accessor function-matrix-vref
+ :type function)
+ (set-vref
+ :initarg :set-vref
+ :initform (constantly nil)
+ :accessor function-matrix-set-vref
+ :type function))
+ (:documentation "Matrix without a store."))
Modified: src/matrix/level1-constructors.lisp
==============================================================================
--- src/matrix/level1-constructors.lisp (original)
+++ src/matrix/level1-constructors.lisp Sun May 17 15:02:03 2009
@@ -19,14 +19,14 @@
(in-package :lisplab)
-(export '(mat new col row))
+#+nil (export '(mat new col row))
-(export '(rmat rnew rcol rrow))
+(export '(funmat
+ rmat rnew rcol rrow
+ cmat cnew ccol crow))
-(export '(cmat cnew ccol crow))
-
-(defmacro mat (type &body args)
+#+nil (defmacro mat (type &body args)
"Creates a matrics"
`(convert
,(cons 'list (mapcar (lambda (x)
@@ -34,11 +34,11 @@
args))
,type))
-(defun col (type &rest args)
+#+nil (defun col (type &rest args)
"Creates a column matrix"
(convert (mapcar 'list args) type))
-(defun row (type &rest args)
+#+nil (defun row (type &rest args)
"Creates a row matrix"
(convert args type))
@@ -62,7 +62,7 @@
(defun rnew (value rows &optional (cols 1))
"Creates a new blas-real matrix"
- (new 'matrix-dge (list rows cols) t value))
+ (mnew 'matrix-dge value rows cols))
(defmacro cmat (&body args)
"Creates a blas-complex matrics"
@@ -83,4 +83,27 @@
(defun cnew (value rows &optional (cols 1))
"Create a new blas-complex matrix"
- (new 'matrix-zge (list rows cols) t value))
\ No newline at end of file
+ (mnew 'matrix-zge value rows cols))
+
+
+;;; Function matrix
+
+(defmacro funmat (rows cols args &body body)
+ "Creates a read only function matrix"
+ (let ((rows2 (gensym "rows"))
+ (cols2 (gensym "cols"))
+ (i (gensym))
+ (r (gensym))
+ (c (gensym)))
+ `(let ((,rows2 ,rows)
+ (,cols2 ,cols))
+ (make-instance 'function-matrix
+ :rows ,rows2
+ :cols ,cols2
+ :mref (lambda (self , at args)
+ (declare (muffle-conditions style-warning))
+ , at body)
+ :vref (lambda (self ,i)
+ ;; Default self vector reference in column major order
+ (multiple-value-bind (,r ,c) (floor ,i ,rows2)
+ (mref self ,r ,c)))))))
Modified: src/matrix/level1-generic.lisp
==============================================================================
--- src/matrix/level1-generic.lisp (original)
+++ src/matrix/level1-generic.lisp Sun May 17 15:02:03 2009
@@ -26,7 +26,7 @@
#-todo-remove(defmethod new (class dim &optional (element-type t) (value 0))
;;; TODO get rid of this default that calls the new constructor
- (mnew class dim element-type value))
+ (mnew class value (car dim) (cadr dim)))
#+todo-remove(defmethod convert (obj type)
(if (not (or (vector? obj) (matrix? obj)))
Modified: src/matrix/level1-interface.lisp
==============================================================================
--- src/matrix/level1-interface.lisp (original)
+++ src/matrix/level1-interface.lisp Sun May 17 15:02:03 2009
@@ -39,7 +39,7 @@
(defgeneric new (class dim &optional element-type value)
(:documentation "Deprecated. Use mnew in stead. Creates a new matrix filled with numeric arguments."))
-(defgeneric mnew (class dim &optional element-type value)
+(defgeneric mnew (class value rows &optional cols)
(:documentation "General matrix constructor. Creates a new matrix filled with numeric arguments."))
(defgeneric ref (matrix &rest subscripts)
Modified: src/matrix/level1-matrix.lisp
==============================================================================
--- src/matrix/level1-matrix.lisp (original)
+++ src/matrix/level1-matrix.lisp Sun May 17 15:02:03 2009
@@ -20,7 +20,38 @@
(in-package :lisplab)
-;;; Generic methods
+;;; This is OK, but could be optimzied!
+(defmacro w/mat (a args &body body)
+ (let ((a2 (gensym))
+ (x (first args))
+ (i (second args))
+ (j (third args)))
+ `(let ((,a2 ,a))
+ (dotimes (,i (rows ,a2))
+ (dotimes (,j (cols ,a2))
+ (let ((,x (mref ,a2 ,i ,j)))
+ (setf (mref ,a2 ,i ,j)
+ , at body))))
+ ,a2)))
+
+
+;;; Generic methods and functions
+
+(defun convert-list-to-matrix (list type)
+ (let* ((rows (length list))
+ (cols (length (car list)))
+ (m (mnew type 0 rows cols)))
+ (fill-matrix-with-list m list)))
+
+(defun convert-matrix-to-matrix (m0 type)
+ (let* ((rows (rows m0))
+ (cols (cols m0))
+ (m (mnew type 0 rows cols)))
+ (dotimes (i rows)
+ (dotimes (j cols)
+ (setf (mref m i j) (mref m0 i j))))
+ m))
+
(defmethod scalar? ((x matrix-base)) nil)
@@ -65,11 +96,20 @@
(make-matrix-instance 'matrix-zge dim value)
(make-matrix-instance 'matrix-dge dim value)))
+
+
+
+
;;; Spcialized for blas-dge
-(defmethod mnew ((class (eql 'matrix-dge)) dim &optional (element-type t) (value 0))
- (declare (ignore element-type))
- (make-matrix-instance class dim value))
+(defmethod convert ((x cons) (type (eql 'matrix-dge)))
+ (convert-list-to-matrix x type))
+
+(defmethod convert ((x matrix-base) (type (eql 'matrix-dge)))
+ (convert-matrix-to-matrix x type))
+
+(defmethod mnew ((class (eql 'matrix-dge)) value rows &optional (cols 1))
+ (make-matrix-instance class (list rows cols) value))
(defmethod mref ((matrix matrix-base-dge) row col)
(aref (the type-blas-store (matrix-store matrix))
@@ -91,11 +131,18 @@
(setf (aref (the type-blas-store (matrix-store matrix)) idx)
(the double-float (coerce value 'double-float))))
+
+
;;; Spcialized for blas-zge
-(defmethod mnew ((class (eql 'matrix-zge)) dim &optional (element-type t) (value 0))
- (declare (ignore element-type))
- (make-matrix-instance class dim value))
+(defmethod convert ((x cons) (type (eql 'matrix-zge)))
+ (convert-list-to-matrix x type))
+
+(defmethod convert ((x matrix-base) (type (eql 'matrix-zge)))
+ (convert-matrix-to-matrix x type))
+
+(defmethod mnew ((class (eql 'matrix-zge)) value rows &optional (cols 1))
+ (make-matrix-instance class (list rows cols) value))
(defmethod mref ((matrix matrix-base-zge) row col)
(ref-blas-complex-store (matrix-store matrix)
@@ -110,11 +157,24 @@
value)
(defmethod vref ((matrix matrix-base-zge) i)
- (ref-blas-complex-store (store matrix) i 0 1))
+ (ref-blas-complex-store (matrix-store matrix) i 0 1))
(defmethod (setf vref) (value (matrix matrix-base-zge) i)
(setf (ref-blas-complex-store (matrix-store matrix) i 0 1)
(coerce value '(complex double-float)))
value)
+;;; Function matrix
+
+(defmethod mref ((f function-matrix) row col)
+ (funcall (function-matrix-mref f) f row col))
+
+(defmethod (setf mref) (value (f function-matrix) row col)
+ (funcall (function-matrix-set-mref f) value f row col))
+
+(defmethod vref ((f function-matrix) idx)
+ (funcall (function-matrix-vref f) f idx))
+
+(defmethod (setf vref) (value (f function-matrix) idx)
+ (funcall (function-matrix-set-vref f) value f idx))
Modified: src/matrix/level1-util.lisp
==============================================================================
--- src/matrix/level1-util.lisp (original)
+++ src/matrix/level1-util.lisp Sun May 17 15:02:03 2009
@@ -29,6 +29,18 @@
:rows rows
:cols cols)))
+(defun fill-matrix-with-list (m x)
+ (let* ((rows (rows m))
+ (cols (cols m)))
+ (do ((xx x (cdr xx))
+ (i 0 (1+ i)))
+ ((= i rows))
+ (do ((yy (car xx) (cdr yy))
+ (j 0 (1+ j)))
+ ((= j cols))
+ (setf (mref m i j) (car yy))))
+ m))
+
(deftype type-blas-store ()
'(simple-array double-float (*)))
Modified: src/matrix/level2-generic.lisp
==============================================================================
--- src/matrix/level2-generic.lisp (original)
+++ src/matrix/level2-generic.lisp Sun May 17 15:02:03 2009
@@ -17,6 +17,8 @@
;;; with this program; if not, write to the Free Software Foundation, Inc.,
;;; 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+;;; TODO clean this up
+
(in-package :lisplab)
(defmethod square-matrix? (x)
Modified: src/matrix/level2-matrix-dge.lisp
==============================================================================
--- src/matrix/level2-matrix-dge.lisp (original)
+++ src/matrix/level2-matrix-dge.lisp Sun May 17 15:02:03 2009
@@ -19,6 +19,11 @@
(in-package :lisplab)
+(defmethod fill! ((a matrix-dge) value)
+ (let ((x (coerce value 'double-float))
+ (store (matrix-store a)))
+ (fill store x)))
+
(defmethod copy ((matrix matrix-base-dge))
(make-instance (class-name (class-of matrix))
:store (copy-seq (matrix-store matrix))
@@ -94,6 +99,6 @@
(matrix-store b)
(lambda (&rest args)
(coerce (apply f args) 'double-float))
- (matrix-store a) (mapcar #'store args))
+ (matrix-store a) (mapcar #'matrix-store args))
b))
Modified: src/matrix/level2-matrix-zge.lisp
==============================================================================
--- src/matrix/level2-matrix-zge.lisp (original)
+++ src/matrix/level2-matrix-zge.lisp Sun May 17 15:02:03 2009
@@ -18,6 +18,15 @@
(in-package :lisplab)
+(defmethod fill! ((a matrix-zge) value)
+ (let ((rx (coerce (realpart value) 'double-float))
+ (cx (coerce (imagpart value) 'double-float))
+ (store (matrix-store a)))
+ (loop for i from 0 below (length store) by 2 do
+ (setf (aref store i) rx
+ (aref store (1+ i)) cx))))
+
+
(defmethod copy ((matrix matrix-base-zge))
(make-instance (class-name (class-of matrix))
:store (copy-seq (matrix-store matrix))
Modified: system/lisplab.asd
==============================================================================
--- system/lisplab.asd (original)
+++ system/lisplab.asd Sun May 17 15:02:03 2009
@@ -41,11 +41,12 @@
(:file "level1-constructors")
(:file "level2-interface")
-
+ (:file "level2-generic")
+ (:file "level2-array-functions")
(:file "level2-matrix-dge")
(:file "level2-matrix-zge")
- (:file "level2-array-functions")
+
; (:file "level1-interface")
; (:file "level1-util")
@@ -82,7 +83,7 @@
;;
;; Linear algebra lisp implementation (Level 3)
;;
- #+nil (:module :linalg-native
+ (:module :linalg-native
:depends-on (:matrix :linalg-interface)
:pathname "../src/linalg/"
:serial t
@@ -94,14 +95,14 @@
;;
;; Fast Fourier transform (Level 3)
;;
- #+nil (:module :fft
+ (:module :fft
:depends-on (:matrix)
:pathname "../src/fft/"
:serial t
:components
(
(:file "level3-fft-interface")
- (:file "level3-fft-generic")
+ #+nil (:file "level3-fft-generic")
(:file "level3-fft-blas")))
;;
@@ -119,7 +120,7 @@
;;
;; Blas and Lapack implmentations (Level 3)
;;
- #+nil (:module :matlisp
+ (:module :matlisp
:depends-on (:matrix :linalg-interface)
:pathname "../src/matlisp/"
:serial t
More information about the lisplab-cvs
mailing list