;;; -*- Mode: lisp; Syntax: ansi-common-lisp; Package: :lisplab; Base: 10 -*-

;;; Copyright (C) 2009 Knut S. Gjerden
;;; 
;;; 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.

;;; $Log: potrf.lisp,v $
;;; Revision 1.2   14.01.2010 13:28:40 knutgj
;;; Revision 1.1   06.08.2009 12:40:40 knutgj
;;; o Initial revision.


(in-package :lisplab)

(defgeneric potrf! (a &key uplo)
  (:documentation
"
  Syntax
  ======
  (POTRF a [:u :L])

  Purpose
  =======
  DPOTRF computes the Cholesky factorization of a real symmetric
  positive definite matrix A.

  The factorization has the form
     A = U**T * U,  if UPLO = 'U', or
     A = L  * L**T,  if UPLO = 'L',
  where U is an upper triangular matrix and L is lower triangular.

  This is the block version of the algorithm, calling Level 3 BLAS.

  Return Values
  =============
  [1] The factor U or L from the Cholesky
          factorization A = U**T*U or A = L*L**T.
  [2] INFO = T: successful
             i:  U(i,i) is exactly zero. 
"))



(defmethod potrf! ((a matrix-foreign-dge) &key uplo)
  (let* ((n (rows a))
	 (m (cols a)))

    (declare (type fixnum n m))
    (multiple-value-bind (new-a info)
	(f77::dpotrf (case uplo
                  (:L "L")
                  (:U "U")
                  (t "U")) ;; UPLO
		m          ;; N
		(vector-store a)  ;; A
		n          ;; LDA
		0)         ;; INFO
      (declare (ignore new-a))
      (values a (if (zerop info)
			 t
		       info)))))

(defmethod potrf! ((a matrix-foreign-zge) &key uplo)
  (let* ((n (rows a))
	 (m (cols a)))

    (declare (type fixnum n m))
    (multiple-value-bind (new-a info)
	(f77::zpotrf  (case uplo
                  (:L "L")
                  (:U "U")
                  (t "U")) ;; UPLO
		m          ;; N
		(vector-store a)  ;; A
		n          ;; LDA
		0)         ;; INFO
      (declare (ignore new-a))
      (values a (if (zerop info)
			 t
		       info)))))