;;; -*- 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: potrs.lisp,v $ 
;;; Revision 1.2   14.01.2010 13:29:27 knutgj
;;; Revision 1.1   06.08.2009 12:19:27 knutgj
;;; o Initial revision.


(in-package :lisplab)

(defgeneric potrs! (a b &key uplo)
  (:documentation
   "
  Syntax
  ======
  (POTRS! a b [:U :L])

  Purpose
  =======
  Solves a system of linear equations
      A * X = B  or  A' * X = B
  with a general N-by-N matrix A using the Cholesky LU factorization computed
  by POTRF.  A and are the results from POTRF, UPLO specifies
  the form of the system of equations: 
           = 'U':   A = U**T*U 
           = 'L':   A = L*L**T

  Return Values
  =============
  [1] The NxM matrix X. (overwriting B)
  [4] INFO = T: successful
             i:  U(i,i) is exactly zero.  The LU factorization
                 used in the computation has been completed, 
                 but the factor U is exactly singular.
                 Solution could not be computed.
"))

(defgeneric potrs (a b &key uplo)
  (:documentation
 "
  Syntax
  ======
  (POTRS a b [:U :L])

  Purpose
  =======
  Same as POTRS! except that B is not overwritten.
"))

(defmethod potrs! :before ((a matrix-foreign-dge) (b matrix-foreign-dge) &key uplo)
  (declare (ignore uplo))
  (let ((n-a (rows a))
        (m-a (cols a))
        (n-b (rows b)))
    (if (not (= n-a m-a n-b))
        (error "Dimensions of A,B given to POTRS! do not match"))))

(defmethod potrs! ((a matrix-foreign-dge) (b matrix-foreign-dge) &key uplo)
  (let* ((n (rows a))
         (m (cols b)))
    (multiple-value-bind (x info)
        (f77::dpotrs (case uplo
                  (:L "L")
                  (:U "U")
                  (t "U"))  ;;UPLO
                n           ;;N
                m           ;;NRHS
                (vector-store a)   ;;A
                n           ;;LDA
                (vector-store b)   ;;B
                n           ;;LDB
                0)          ;;INFO
        (values 
         (make-instance 'matrix-dge :dim  (list n m) :store x)
         (if (zerop info)
             t
           info)))))

(defmethod potrs! ((a matrix-foreign-zge) (b matrix-foreign-zge) &key uplo)

  (let* ((n (rows a))
         (m (cols b)))
    (multiple-value-bind (x info)
        (f77::zpotrs (case uplo
                  (:L "L")
                  (:U "U")
                  (t "U"))
                n
                m
                (vector-store a)
                n
                (vector-store b)
                n
                0)        
        (values 
         (make-instance 'matrix-zge :dim (list n m) :store x) 
         (if (zerop info)
             t
           info)))))

(defmethod potrs ((a matrix-foreign-dge) (b matrix-foreign-dge) &key uplo)
  (potrs! a (copy b) :uplo uplo))

(defmethod potrs ((a matrix-foreign-zge) (b matrix-foreign-zge) &key uplo)
  (potrs! a (copy b) :uplo uplo))