Permalink
Browse files

cmm4, small fixes

  • Loading branch information...
1 parent 3ac5166 commit 00ec6e56e11bcf49aff3e4ddde425f60f320caea @allchemist committed Jun 1, 2011
Showing with 43 additions and 18 deletions.
  1. +35 −12 applications/matrix-product.lisp
  2. +5 −3 applications/svd.lisp
  3. +3 −3 lapacke/highlevel.lisp
@@ -1,26 +1,24 @@
(in-package :sb-math)
-#|
+(export '(mm cmm orthogonal-bais-transform general-basis-transform))
+
;; generic matrix product
-(defgeneric matrix-product (m1 m2 &rest keys &key dest &allow-other-keys))
+(defgeneric mm (m1 m2 &rest keys &key dest &allow-other-keys))
-(defmethod matrix-product ((m1 vector) (m2 vector)
+(defmethod mm ((m1 vector) (m2 vector)
&rest keys &key dest &allow-other-keys)
- (print 'vector-to-vector)
(apply #'ger m1 m2 :dest dest keys))
-(defmethod matrix-product ((m1 array) (m2 vector)
+(defmethod mm ((m1 array) (m2 vector)
&rest keys &key dest &allow-other-keys)
- (print 'matrix-to-vector)
- (apply #'gemv m1 m2 :dest dest keys))
+ (apply #'gemv m1 m2 :dest dest keys))
-(defmethod matrix-product ((m1 array) (m2 array)
+(defmethod mm ((m1 array) (m2 array)
&rest keys &key dest &allow-other-keys)
- (print 'matrix-to-matrix)
(apply #'gemm m1 m2 :dest dest keys))
-(defmethod matrix-product ((m1 triangular-matrix) (m2 vector)
+#|(defmethod matrix-product ((m1 triangular-matrix) (m2 vector)
&rest keys &key dest &allow-other-keys)
(print 'triangular-to-vector)
(apply #'tpmv m1 (if dest (copy m2 dest) (copy m2)) keys))
@@ -34,8 +32,7 @@
(defmethod matrix-product ((m1 hermitian-matrix) (m2 vector)
&rest keys &key dest &allow-other-keys)
(print 'hermitian-to-vector)
- (apply #'hpmv m1 m2 :dest dest keys))
-|#
+ (apply #'hpmv m1 m2 :dest dest keys))|#
;; basis rotation
@@ -47,10 +44,36 @@
;; chain matrix multiplication
+(defun cmm (&rest matrices)
+ (ecase (length matrices)
+ (0 nil)
+ (1 (first matrices))
+ (2 (gemm matrices))
+ (3 (cmm3 matrices))
+ (4 (cmm4 matrices))))
+
(defun cmm3 (m1 m2 m3)
(if (< (+ (* (dim0 m1) (dim1 m1) (dim1 m2))
(* (dim0 m1) (dim1 m2) (dim1 m3)))
(+ (* (dim0 m2) (dim1 m2) (dim1 m3))
(* (dim0 m1) (dim1 m1) (dim1 m3))))
(gemm (gemm m1 m2) m3)
(gemm m1 (gemm m2 m3))))
+
+(defun cmm4 (m1 m2 m3 m4)
+ (let ((r1 (dim0 m1)) (c1 (dim1 m1))
+ (r2 (dim0 m2)) (c2 (dim1 m2))
+ (r3 (dim0 m1)) (c3 (dim1 m1))
+ (r4 (dim0 m2)) (c4 (dim1 m2))
+ (vals (make-list 6)))
+ (setf (elt vals 0) (+ (* r3 r4 c4) (* r2 r3 c4) (* r1 r2 c4)))
+ (setf (elt vals 1) (+ (* r2 r3 c3) (* r2 r4 c4) (* r1 r2 c4)))
+ (setf (elt vals 2) (+ (* r2 r3 c3) (* r1 r2 c3) (* r1 r4 c4)))
+ (setf (elt vals 3) (+ (* r1 r2 c2) (* r1 r3 c3) (* r1 r4 c4)))
+ (setf (elt vals 4) (+ (* r1 r2 c2) (* r3 r4 c4) (* r1 r3 c4)))
+ (ecase (position (apply #'min vals) vals :test #'eq)
+ (0 (gemm m1 (gemm m2 (gemm m3 m4))))
+ (1 (gemm m1 (gemm (gemm m2 m3) m4)))
+ (2 (gemm (gemm m1 (gemm m2 m3)) m4))
+ (3 (gemm (gemm (gemm m1 m2) m3) m4))
+ (4 (gemm (gemm m1 m2) (gemm m3 m4))))))
View
@@ -1,9 +1,11 @@
(in-package :sb-math)
-(defun pseudo-inverse (A)
+(export '(pinv em-norm fm-norm cond-number llsq))
+
+(defun pinv (A)
(multiple-value-bind (S U VT)
(svd A :left :all :right :all)
- (gemm
+ (gemm
(gemm VT (setf (diag (make-matrix-like A)) (map-matrix-/ S)) :transa :trans :transb :trans)
U :transb :trans)))
@@ -20,5 +22,5 @@
(error "Matrix is unconditioned: last singular value is zero")
(/ (aref S 0) last-sv))))
-(defun linear-least-squares (A Y)
+(defun llsq (A Y)
(gemv (pseudo-inverse A) Y))
View
@@ -1,6 +1,6 @@
(in-package :sb-math)
-(export '(transpose lu lu-inverse lu-solve svd))
+(export '(transpose lu inv solve svd))
(defun transpose (matrix &optional dest)
(let ((element-type (array-element-type matrix))
@@ -32,7 +32,7 @@
(assert (zerop info) nil "LU decomposition returned non-zero code: ~A" info)
(values %A ipiv))))
-(defun lu-inverse (A &key ipiv (pure t))
+(defun inv (A &key ipiv (pure t))
(let ((dim (dim0 A))
(element-type (array-element-type A))
(info -1)
@@ -54,7 +54,7 @@
(assert (zerop info) nil "LU inversion returned non-zero code: ~A" info))))
%A))
-(defun lu-solve (A B &key (trans :notrans) (pure t))
+(defun solve (A B &key (trans :notrans) (pure t))
(let* ((dim (dim0 A))
(dim1 (if (vectorp B) 1 (dim1 B)))
(element-type (array-element-type A))

0 comments on commit 00ec6e5

Please sign in to comment.