;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;
;;; LUSH Lisp Universal Shell
;;;   Copyright (C) 2002 Leon Bottou, Yann Le Cun, AT&T Corp, NECI.
;;; Includes parts of TL3:
;;;   Copyright (C) 1987-1999 Leon Bottou and Neuristique.
;;; Includes selected parts of SN3.2:
;;;   Copyright (C) 1991-2001 AT&T Corp.
;;;
;;; 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., 59 Temple Place, Suite 330, Boston, MA 02111, USA
;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;; $Id: linear-equation.lsh,v 1.3 2003/06/04 18:40:44 leonb Exp $
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;; by Fu Jie Huang, May 2003

;; use the routine "dgesv" to solve a linear equation Ax=b
;;  [ 2   -1   1] [ ] [-1]
;;  [ 3    3   9]*[x]=[ 0]
;;  [ 3    3   5] [ ] [ 4]
;; obviously, 
;;  [ 1]
;;x=[ 2] 
;;  [-1]

(libload "lapack/lapack-d")
(setq A [d [2 3 3][-1 3 3][1 9 5]])  ;; row-wise storage
(setq b [d [-1 0 4]])                ;; row-wise storage
(setq ipiv (int-matrix 3))

(let ((n    ((int-matrix) 3))
      (nrhs ((int-matrix) 1))
      (lda  ((int-matrix) 3))
      (ldb  ((int-matrix) 3))
      (info ((int-matrix) 0)))
  (dgesv n nrhs A lda ipiv b ldb info))

(print "solution x=")
(pretty b)        ;; the output "x" is stored in "b"

;; try to peek into what happened during the procedure
(print "A=P*L*U:")
(pretty A)

(print "pivot indices:")
(pretty ipiv)

;; A is LU decomposed with pivoting:   A = P * L * U
;;     [  1  0  0]        [3   3   9]
;; L = [2/3  1  0]    U = [0  -3  -5]   
;;     [  1  0  1]        [0   0  -4]
;; L and U are packed inside A
;; the permutation matrix P can be derived from ipiv
