[lisplab-cvs] r134 - trunk/src/matlisp

[Jørn Inge Vestgården]

Author: jivestgarden
Date: Wed Mar 10 14:52:37 2010
New Revision: 134

Log:
lapack lu not yet ok

Added:
   trunk/src/matlisp/lu.lisp
Modified:
   trunk/src/matlisp/inv.lisp

Modified: trunk/src/matlisp/inv.lisp
==============================================================================
--- trunk/src/matlisp/inv.lisp	(original)
+++ trunk/src/matlisp/inv.lisp	Wed Mar 10 14:52:37 2010
@@ -1,4 +1,4 @@
-;;; Lisplab, matliap/div.lisp
+;;; Lisplab, matlisp/div.lisp
 ;;; Lapack-based matrix inversion
 
 ;;; Copyright (C) 2009 Joern Inge Vestgaarden
@@ -23,7 +23,7 @@
   (if cl-user::*lisplab-liblapack-path*
       (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"))
+	     (msg "Argument A given to minv is singular to working machine precision"))
 	(multiple-value-bind (_ ipiv info) 
 	    (f77::dgetrf N N (matrix-store a) N ipiv 0)
 	  (declare (ignore _))
@@ -46,7 +46,7 @@
   (if cl-user::*lisplab-liblapack-path*
       (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"))
+	     (msg "Argument A given to mdiv is singular to working machine precision"))
 	(multiple-value-bind (_ ipiv info) 
 	    (f77::zgetrf N N (matrix-store a) N ipiv 0)
 	  (declare (ignore _))

Added: trunk/src/matlisp/lu.lisp
==============================================================================
--- (empty file)
+++ trunk/src/matlisp/lu.lisp	Wed Mar 10 14:52:37 2010
@@ -0,0 +1,51 @@
+;;; Lisplab, matlisp/lu.lisp
+;;; LU-factorization
+
+;;; Copyright (C) 2009 Joern Inge Vestgaarden
+;;;
+;;; This program is free software; you can redistribute it and/or modify
+;;; it under the terms of the GNU General Public License as published by
+;;; the Free Software Foundation; either version 2 of the License, or
+;;; (at your option) any later version.
+;;;
+;;; This program is distributed in the hope that it will be useful,
+;;; but WITHOUT ANY WARRANTY; without even the implied warranty of
+;;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+;;; GNU General Public License for more details.
+;;;
+;;; You should have received a copy of the GNU General Public License along
+;;; with this program; if not, write to the Free Software Foundation, Inc.,
+;;; 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+;;; TODO: maybe speed up.
+
+(in-package :lisplab)
+
+(defmethod LU-factor! ((A matrix-blas-dge) p)
+  (if cl-user::*lisplab-liblapack-path*
+      (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))))
+      (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