;;; Lisplab, matlisp/div.lisp
;;; Lapack-based matrix inversion

;;; 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.

(in-package :lisplab)

(defmethod minv! ((a matrix-foreign-dge))
  (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))
	  (let ((work (make-array N :element-type 'double-float)))
	    (multiple-value-bind (_ __ info)	  
		(f77::dgetri N (vector-store a) N ipiv work N 0)
	      (declare (ignore _ __))
	      (unless (zerop info)
		(error msg))
	      a))))
      ;; Othervise, call native lisp implementation
      (call-next-method)))
			
(defmethod minv ((a matrix-foreign-dge))
  (minv! (copy a)))

(defmethod minv! ((a matrix-foreign-zge))
  (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"))
	(multiple-value-bind (_ ipiv info) 
	    (f77::zgetrf N N (vector-store a) N ipiv 0)
	  (declare (ignore _))
	  (unless (zerop  info)
	    (error msg  ))
	  (let ((work (make-array (* 2 N) :element-type 'double-float)))
	    (multiple-value-bind (_ __ info)
		(f77::zgetri N (vector-store a) N ipiv work N 0)
	      (declare (ignore _ __))
	      (unless (zerop info)
		(error msg))
	      a))))
      ;; Othervise, call native lisp implementation
      (call-next-method)))
			
(defmethod minv ((a matrix-foreign-zge))
  (minv! (copy a)))