[lisplab-cvs] r168 - in trunk/src: linalg test

[Jørn Inge Vestgården]

Author: jivestgarden
Date: Mon May 17 14:31:01 2010
New Revision: 168

Log:
bugfix lin-solve

Modified:
   trunk/src/linalg/level3-linalg-dge.lisp
   trunk/src/linalg/level3-linalg-generic.lisp
   trunk/src/linalg/level3-linalg-interface.lisp
   trunk/src/test/test-methods.lisp

Modified: trunk/src/linalg/level3-linalg-dge.lisp
==============================================================================
--- trunk/src/linalg/level3-linalg-dge.lisp	(original)
+++ trunk/src/linalg/level3-linalg-dge.lisp	Mon May 17 14:31:01 2010
@@ -32,10 +32,10 @@
 	(setf (mref b i j) (mref a j i))))
     b))
 
-(defmethod mct ((a matrix-lisp-dge))
+(defmethod mct ((a matrix-base-dge))
   (mtp a))
 
-(defmethod m* ((a matrix-lisp-dge) (b matrix-lisp-dge))
+(defmethod m* ((a matrix-base-dge) (b matrix-base-dge))
   (let* ((N (rows a))
 	 (M (cols b))
 	 (S (rows b))
@@ -57,7 +57,7 @@
 	    (setf (refc i j) cij)))) 
       c)))
 
-(defmethod LU-factor! ((A matrix-lisp-dge) p)
+(defmethod LU-factor! ((A matrix-base-dge) p)
   ;; Translation from GSL. 
   ;; Destructive LU factorization. The outout is PA=LU,
   ;; stored in one matrix, where the diagonal elements belong
@@ -100,7 +100,7 @@
 	  
 (defun L-solve!-blas-real (L x col)
   ;; Solve Lx=b
-  (declare (matrix-lisp-dge L x))
+  (declare (matrix-base-dge L x))
   (let ((L2 (matrix-store L))
 	(x2 (matrix-store x))
 	(N (rows x)))
@@ -117,7 +117,7 @@
   x)
 
 (defun U-solve!-blas-real (U x col)
-  (declare (matrix-lisp-dge U x))
+  (declare (matrix-base-dge U x))
   (let* ((U2 (matrix-store U))
 	 (x2 (matrix-store x))
 	 (N (rows x))
@@ -141,8 +141,9 @@
   (U-solve!-blas-real LU x col)
   x)
 
-(defmethod lin-solve ((A matrix-lisp-dge) (b matrix-lisp-dge))
-  (destructuring-bind (LU pvec sign) (LU-factor A)
+(defmethod lin-solve ((A matrix-base-dge) (b matrix-base-dge))
+  (destructuring-bind (LU pvec sign) 
+      (LU-factor! (copy A) (make-permutation-vector (size b)))
     (let ((b2 (copy b)))
       (dotimes (i (rows A))
 	(setf (vref b2 (vref pvec i)) (vref b i)))
@@ -159,10 +160,10 @@
 	(LU-solve!-blas-real LU A  (vref p i)))))
   A)
       
-(defmethod minv! ((A matrix-lisp-dge))
+(defmethod minv! ((A matrix-base-dge))
   (minv!-blas-real A))
 
-(defmethod minv ((A matrix-lisp-dge))
+(defmethod minv ((A matrix-base-dge))
   (minv! (copy A)))
 
 

Modified: trunk/src/linalg/level3-linalg-generic.lisp
==============================================================================
--- trunk/src/linalg/level3-linalg-generic.lisp	(original)
+++ trunk/src/linalg/level3-linalg-generic.lisp	Mon May 17 14:31:01 2010
@@ -160,17 +160,18 @@
 				(mref U i i)))))
       x))
 
-(defun LU-solve! (LU x)
+(defmethod LU-solve! ((LU matrix-base) (x matrix-base))
   (L-solve! LU x)
   (U-solve! LU x)
   x)
 
 (defmethod lin-solve ((A matrix-base) (b matrix-base))
-  (destructuring-bind (LU pvec sign) (LU-factor A)
+  (destructuring-bind (LU pvec det) 
+      (LU-factor! (copy A) (make-permutation-vector (rows A)))    
     (let ((b2 (copy b)))
-      (dotimes (i (rows A))
-	(setf (vref b2 (vref pvec i)) (vref b i)))
-      (LU-solve! LU b2))))
+    (dotimes (i (size b))
+      (setf (vref b2 i) (vref b (vref pvec i)))) 
+    (LU-solve! LU b2))))
   
 (defmethod mdet ((A matrix-base))
   (destructuring-bind (LU pvec det) 

Modified: trunk/src/linalg/level3-linalg-interface.lisp
==============================================================================
--- trunk/src/linalg/level3-linalg-interface.lisp	(original)
+++ trunk/src/linalg/level3-linalg-interface.lisp	Mon May 17 14:31:01 2010
@@ -71,6 +71,10 @@
 L is low diagonal with unity at diagnoals, U is upper diagnoal and 
 P is an permutation matrix, so that A = P L U."))
 
+(defgeneric LU-solve! (LU x)
+  (:documentation "Solves the linear equation system LUx=b, where LU is an output 
+from LU-factor!."))
+
 (defgeneric lin-solve (A b)
   (:documentation  "Solves the linear system of equations Ax=b."))
 

Modified: trunk/src/test/test-methods.lisp
==============================================================================
--- trunk/src/test/test-methods.lisp	(original)
+++ trunk/src/test/test-methods.lisp	Mon May 17 14:31:01 2010
@@ -100,7 +100,11 @@
 	  (c #md((1 2 -1) (3 4 9) (1 1 1)))
 	  (d #mz((1 2 2.1) (3 5 %i) (-%i %i -%i)))
 	  (x #mm((1 2 2.1) (3 5 %i) (-%i %i -%i)))
-	  (args (list a b c d x)))
+	  (args (list a b c d x))
+	  (vec1 (mcol 'matrix-ge 1 2 3))
+	  (vec2 (dcol 1 2 3)))
+     (simple-non-nil-check #'lin-solve (list x vec1))
+     (simple-non-nil-check #'lin-solve (list c vec2))
      (mapc (lambda (x) (simple-non-nil-check #'mtp (list x))) args)
      (mapc (lambda (x) (simple-non-nil-check #'mct (list x))) args)
      (mapc (lambda (x) (simple-non-nil-check #'minv (list x))) args)