[lisplab-cvs] r210 - in trunk/src: linalg matrix2
jivestgarden at common-lisp.net
jivestgarden at common-lisp.net
Fri Mar 23 19:19:58 UTC 2012
Author: jivestgarden
Date: Fri Mar 23 12:19:57 2012
New Revision: 210
Log:
Added level2 column operations
Modified:
trunk/src/linalg/level3-linalg-generic.lisp
trunk/src/matrix2/matrix2-dge.lisp
trunk/src/matrix2/matrix2-generic.lisp
trunk/src/matrix2/matrix2-interface.lisp
Modified: trunk/src/linalg/level3-linalg-generic.lisp
==============================================================================
--- trunk/src/linalg/level3-linalg-generic.lisp Sun Mar 18 08:40:45 2012 (r209)
+++ trunk/src/linalg/level3-linalg-generic.lisp Fri Mar 23 12:19:57 2012 (r210)
@@ -38,7 +38,7 @@
(defmethod mct ((a matrix-base))
(.conj (mtp a)))
-(defmethod m* ((a matrix-base) (b matrix-base))
+#+old (defmethod m* ((a matrix-base) (b matrix-base))
(let ((c (mcreate a 0 (list (rows a) (cols b)))))
(dotimes (i (rows c))
(dotimes (j (cols c))
@@ -47,6 +47,15 @@
(.* (mref a i k) (mref b k j)))))))
c))
+(defmethod m* ((a matrix-base) (b matrix-base))
+ (let ((c (mcreate a 0 (list (rows a) (cols b))))
+ (a (mtp a)))
+ (dotimes (i (rows c))
+ (dotimes (j (cols c))
+ (setf (mref c i j)
+ (col-col-mul-sum A i b j))))
+ c))
+
(defmethod m/ ((a matrix-base) (b matrix-base))
(m* a (minv b)))
@@ -64,27 +73,6 @@
(setf (vref col i) 1)
(LU-solve! LU col))))
A))
-
-#+nil (defmethod minv! ((a matrix-base))
- ;; Flawed. Does not work on when pivoting is needed
- "Brute force O(n^3) implementation of matrix inverse.
-Think I'll keep this for the general case since it works also
-when the elements cannot be ordered, unlike the LU-based version"
- (let* ((size (rows a))
- (temp 0))
- (dotimes (i size a)
- (setf temp (mref a i i))
- (dotimes (j size)
- (setf (mref a i j) (if (= i j)
- (./ (mref a i j))
- (./ (mref a i j) temp))))
- (dotimes (j size)
- (unless (= i j)
- (setf temp (mref a j i)
- (mref a j i) 0)
- (dotimes (k size)
- (setf (mref a j k)
- (.- (mref a j k) (.* temp (mref a i k))))))))))
(defmethod LU-factor! ((A matrix-base) p)
;; Translation from GSL.
@@ -182,8 +170,26 @@
-
-
-
+;;; Alternative code
+#+nil (defmethod minv! ((a matrix-base))
+ ;; Flawed. Does not work on when pivoting is needed
+ "Brute force O(n^3) implementation of matrix inverse.
+Think I'll keep this for the general case since it works also
+when the elements cannot be ordered, unlike the LU-based version"
+ (let* ((size (rows a))
+ (temp 0))
+ (dotimes (i size a)
+ (setf temp (mref a i i))
+ (dotimes (j size)
+ (setf (mref a i j) (if (= i j)
+ (./ (mref a i j))
+ (./ (mref a i j) temp))))
+ (dotimes (j size)
+ (unless (= i j)
+ (setf temp (mref a j i)
+ (mref a j i) 0)
+ (dotimes (k size)
+ (setf (mref a j k)
+ (.- (mref a j k) (.* temp (mref a i k))))))))))
Modified: trunk/src/matrix2/matrix2-dge.lisp
==============================================================================
--- trunk/src/matrix2/matrix2-dge.lisp Sun Mar 18 08:40:45 2012 (r209)
+++ trunk/src/matrix2/matrix2-dge.lisp Fri Mar 23 12:19:57 2012 (r210)
@@ -58,3 +58,62 @@
(- j dc)
rows)))))
B))
+
+;;;; The column operations
+
+(defmethod col-smul! ((A matrix-dge) i num)
+ (let* ((num (coerce num 'double-float))
+ (A-store (vector-store A))
+ (r (rows A))
+ (start (* r i))
+ (end (* r (1+ i))))
+ (declare (type type-blas-store A-store)
+ (type type-blas-idx start end))
+ (loop for k from start below end do
+ (setf (aref A-store k)
+ (* (aref A-store k) num))))
+ A)
+
+(defmethod col-swap! ((A matrix-dge) i j)
+ (let* ((A-store (vector-store A))
+ (r (rows A))
+ (tmp 0d0)
+ (ii (* i r))
+ (jj (* j r)))
+ (declare (type type-blas-store A-store)
+ (type type-blas-idx r ii jj)
+ (type double-float tmp))
+ (dotimes (k r)
+ (setf tmp (aref A-store ii)
+ (aref A-store ii) (aref A-store jj)
+ (aref A-store jj) tmp)
+ (incf ii)
+ (incf jj))
+ A))
+
+(defmethod col-sum ((A matrix-dge) i)
+ (let* ((A-store (vector-store A))
+ (r (rows A))
+ (start (* r i))
+ (end (* r (1+ i)))
+ (sum 0d0))
+ (declare (type type-blas-store A-store)
+ (type type-blas-idx i r start end)
+ (type double-float sum))
+ (loop for k from start below end do
+ (incf sum (aref A-store k)))
+ sum))
+
+(defmethod col-col-mul-sum ((A matrix-dge) i
+ (B matrix-dge) j)
+ (let ((A-store (vector-store A))
+ (B-store (vector-store B))
+ (r (rows A))
+ (sum 0d0))
+ (declare (type type-blas-store A-store B-store)
+ (type type-blas-idx i j r)
+ (type double-float sum))
+ (dotimes (k r)
+ (incf sum (* (aref A-store (column-major-idx k i r))
+ (aref B-store (column-major-idx k j r)))))
+ sum))
Modified: trunk/src/matrix2/matrix2-generic.lisp
==============================================================================
--- trunk/src/matrix2/matrix2-generic.lisp Sun Mar 18 08:40:45 2012 (r209)
+++ trunk/src/matrix2/matrix2-generic.lisp Fri Mar 23 12:19:57 2012 (r210)
@@ -141,19 +141,44 @@
(reshape a (list rows (/ (size a) rows) 1)))
-(defmethod row-swap! (A i j)
+(defmethod row-swap! ((A matrix-base) i j)
(dotimes (c (cols A))
(psetf (mref A i c) (mref A j c)
(mref A j c) (mref A i c)))
A)
-(defmethod row-mul! (A i num)
+(defmethod row-mul! ((A matrix-base) i num)
(dotimes (c (cols A))
(setf (mref A i c) (.* num (mref A i c))))
A)
-(defmethod row-add! (A i j num)
+(defmethod row-add! ((A matrix-base) i j num)
(dotimes (c (cols A))
(setf (mref A i c) (.+ (mref A i c) (.* num (mref A j c)))))
A)
+;;; The column operations
+
+(defmethod col-swap! ((A matrix-base) i j)
+ (dotimes (r (rows A))
+ (psetf (mref A r i) (mref A r j)
+ (mref A r j) (mref A r i)))
+ A)
+
+(defmethod col-smul! ((A matrix-base) i num)
+ (dotimes (r (rows A))
+ (setf (mref A r i) (.* num (mref A r i))))
+ A)
+
+(defmethod col-sum ((A matrix-base) i)
+ (let ((sum 0))
+ (dotimes (r (rows A))
+ (setf sum (.+ sum (mref A r i))))
+ sum))
+
+(defmethod col-col-mul-sum ((A matrix-base) i (B matrix-base) j)
+ (let ((sum 0))
+ (dotimes (r (rows A))
+ (setf sum (.+ sum (.* (mref A r i)
+ (mref B r j)))))
+ sum))
Modified: trunk/src/matrix2/matrix2-interface.lisp
==============================================================================
--- trunk/src/matrix2/matrix2-interface.lisp Sun Mar 18 08:40:45 2012 (r209)
+++ trunk/src/matrix2/matrix2-interface.lisp Fri Mar 23 12:19:57 2012 (r210)
@@ -65,7 +65,7 @@
(defgeneric get-col (matrix col)
(:documentation "Gets rows. Destructive"))
-;;; Row operations
+;;; Row operations. TODO remove these
(defgeneric row-swap! (matrix i j)
(:documentation "Swaps row i and j of matrix. Destructive."))
@@ -76,6 +76,21 @@
(defgeneric row-add! (matrix i j number)
(:documentation "Adds a multiplicum of row j to row i. A_ic=A_ic+number*A_jc. Destructive."))
+;;; Column operations.
+
+(defgeneric col-swap! (matrix i j)
+ (:documentation "Swaps row i and j of matrix. Destructive."))
+
+(defgeneric col-smul! (matrix i number)
+ (:documentation "Multiplies row i with a scalar. Destructive."))
+
+(defgeneric col-sum (matrix i)
+ (:documentation "Multiplies row i with a scalar. Destructive."))
+
+(defgeneric col-col-mul-sum (a i b j)
+ (:documentation "The dot product of column i of a and column j of b."))
+
+
;;;; Views
More information about the lisplab-cvs
mailing list