[lisplab-cvs] r138 - in trunk: . src/linalg src/matlisp

[Jørn Inge Vestgården]

Author: jivestgarden
Date: Sat Mar 20 11:25:57 2010
New Revision: 138

Log:
added lu-factorization from lapack

Modified:
   trunk/lisplab.asd
   trunk/package.lisp
   trunk/src/linalg/level3-linalg-generic.lisp
   trunk/src/matlisp/lu.lisp

Modified: trunk/lisplab.asd
==============================================================================
--- trunk/lisplab.asd	(original)
+++ trunk/lisplab.asd	Sat Mar 20 11:25:57 2010
@@ -179,6 +179,7 @@
      (:file "mul")
      (:file "inv")
      (:file "geev")
+     (:file "lu")
      (:file "tridiag")))))
 
 (defsystem :lisplab-fftw

Modified: trunk/package.lisp
==============================================================================
--- trunk/package.lisp	(original)
+++ trunk/package.lisp	Sat Mar 20 11:25:57 2010
@@ -60,9 +60,8 @@
    ;; Some general methods
    "COPY"
    "CONVERT"
-   "SCALAR?"
-   "VECTOR?"
-   "MATRIX?"
+   "VECTOR-P"
+   "MATRIX-P"
 
    ;; Basic methods (The dotted algebra)
    ".+"
@@ -193,6 +192,9 @@
    "MMAX" 
    "MABSMIN" 
    "MABSMAX"
+   "ROW-SWAP!"
+   "ROW-MUL!"
+   "ROW-ADD!"
    "SUB-MATRIX" ; To level3 ?
    "CIRC-SHIFT"
    "PAD-SHIFT"

Modified: trunk/src/linalg/level3-linalg-generic.lisp
==============================================================================
--- trunk/src/linalg/level3-linalg-generic.lisp	(original)
+++ trunk/src/linalg/level3-linalg-generic.lisp	Sat Mar 20 11:25:57 2010
@@ -109,13 +109,11 @@
 		       i-pivot i))))
 	(unless (= i-pivot j)
 	  (dotimes (i N)
-	    (let ((tmp (mref A j i)))
-	      (setf (mref A j i) (mref A i-pivot i) 
-		    (mref A i-pivot i) tmp)))
-	  (let ((tmp (vref p j))) 
-	    (setf (vref p j) (vref p i-pivot) 
-		  (vref p i-pivot) tmp)
-	    (setf sign (- sign)))))
+	    (psetf (mref A j i) (mref A i-pivot i) 
+		   (mref A i-pivot i) (mref A j i)))
+	  (psetf (vref p j) (vref p i-pivot) 
+		 (vref p i-pivot) (vref p j))
+	  (setf sign (- sign))))
       (unless (zerop (mref A j j))
 	(loop for i from (1+ j) below N do
 		(let ((Aij (./ (mref A i j) (mref A j j))))

Modified: trunk/src/matlisp/lu.lisp
==============================================================================
--- trunk/src/matlisp/lu.lisp	(original)
+++ trunk/src/matlisp/lu.lisp	Sat Mar 20 11:25:57 2010
@@ -26,26 +26,20 @@
       (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"))
-	(format t "Hei~%")
 	(multiple-value-bind (_ ipiv info) 
 	    (f77::dgetrf N N (matrix-store a) N ipiv 0)
 	  (declare (ignore _))
 	  (unless (zerop  info)
 	    (error msg))
-
-	  ;; TOOD must change ipiv to a an actual permutation vector !!!!!
-
-
-	  ;; Change from 1 based to zero based index
-	  (dotimes (i (length ipiv))
-	    (setf (aref ipiv i) (1- (aref ipiv i))))		    
-	  (list A ipiv (getrf-get-ipiv-det ipiv))))
+	  ;; Now convert the interchange vector ipiv to the permutation vector p
+	  ;; Also convert from one to zero-indexed.
+	  (let ((det 1)
+		(p (make-permutation-vector (size ipiv))))
+	    (dotimes (i (length ipiv))
+	      (unless (= (1+ i) (aref ipiv i))
+		(setf det (* det -1))
+		(psetf (vref p i) (vref p (1- (vref ipiv i)))
+		       (vref p (1- (vref ipiv i))) (vref p i))))
+	    (list A p det))))
       (call-next-method)))
 
-(defun getrf-get-ipiv-det (ipiv)
-  (let ((det 1))
-    ;; TODO maybe speed up
-    (dotimes (i (length ipiv))
-      (unless (= i (aref ipiv i))
-	(setf det (* det -1))))
-    det))
\ No newline at end of file