;;; 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-foreign-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"))
	(multiple-value-bind (_ ipiv info) 
	    (f77::dgetrf N N (vector-store a) N ipiv 0)
	  (declare (ignore _))
	  (unless (zerop  info)
	    (error msg))
	  ;; 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)))