Re: MCL linear systems

[Dimitri Simos <dim@lissys.demon.co.uk>]

In the MCL CD, look at:

MCL 2.0.1 CD:User Contributed Code:Math & Science:CL Math:matrix.lisp

I have not tried the above.

The following code is based on some fortran code in
"Applied Numerical Analysis", 3rd. Edn., 
Gerald & Wheatley, Addison Wesley.

It is ancient, horrible, a perfect example of how not to do things 
in Lisp, and will raise the blood presure of many a Lisper.
But it works. (I think. Full disclaimers).
And for integers, you get the EXACT solution (although you may
need more RAM than the number of atoms in the visible universe
for any serious matrices).

If you use this for anything serious, do let me know first.
This not shareware.

Dimitri.


(defun ith (i lst)    ;; like nth but with origin 1
  (nth (- i 1) lst))

(defun set-ith (i lst new)
  (setf (nth (- i 1) lst) new))

(defsetf ith set-ith)


(defun check-gauss-elim 
       (lhs solution)
  (let (ans)
    (dolist (sub lhs)
      (setq ans
            (append ans 
                    (list (apply #'+ 
                                 (mapcar #'* sub solution)))))) ans))


(defmacro dofor (var start end &rest body)
  (let ((result (gensym)))
    `(do ((,var ,start (+ 1 ,var))
          (,result nil))
         ((> ,var ,end) ,result)
       (setq ,result (progn ,@body)))))


(defun gauss-aux
       (ab n np)
  (prog* ((nm1 (- n 1)) ratio
          save
          value
          ipvt
          ip1
          nvbl
          l
          answer)
    (declare (optimize speed))
    (dofor i
           1
           nm1
           (setq ipvt i
                 ip1 (+ i 1))
           (dofor j
                  ip1
                  n
                  (if (< (abs (aref ab ipvt i))
                         (abs (aref ab j i)))
                    (setq ipvt j)))
           (if (< (abs (aref ab ipvt i))
                  1.0e-6)
             (error "Sorry, gauss-elimination fails"))
           (when (/= ipvt i)
             (dofor jcol
                    1
                    np
                    (setq save (aref ab i jcol))
                    (setf (aref ab i jcol)
                          (aref ab ipvt jcol)
                          (aref ab ipvt jcol)
                          save)))
           (block thirtytwo
             (dofor jrow
                    ip1
                    n
                    ;          (if (= (aref ab jrow i) 0.0)
                    ;              (return-from thirtytwo
                    ;                nil))
                    (setq ratio (/ (aref ab jrow i)
                                   (aref ab i i)))
                    (setf (aref ab jrow i) ratio)
                    (dofor kcol
                           ip1
                           np
                           (setf (aref ab jrow kcol)
                                 (- (aref ab jrow kcol)
                                    (* ratio
                                       (aref ab i kcol))))))))
    (if (< (abs (aref ab n n))
           1.0e-6)
      (error "Sorry, gauss-elimination fails"))
    (setf (aref ab n np)
          (/ (aref ab n np)
             (aref ab n n)))
    (dofor j
           2
           n
           (setq nvbl (- np j)
                 l (+ nvbl 1)
                 value (aref ab nvbl np))
           (dofor k
                  l
                  n
                  (setq value (- value
                                 (* (aref ab nvbl k)
                                    (aref ab k np)))))
           (setf (aref ab nvbl np)
                 (/ value (aref ab nvbl nvbl))))
    (dofor i
           1
           n
           (setq answer (append answer
                                (list (aref ab i np)))))
    (return answer)))


(defun gauss-elim
       (lhs-rows rhs-column)
  (let* ((n (length rhs-column))
         (np1 (+ n 1))
         (ab (make-array (list (+ 1 n) (+ 1 np1)))))
    (when (/= (length lhs-rows) n)
      (error "inconsistent input to gauss-elim."))
    (mapcar
     #'(lambda (l)
         (when (/= (length l) n)
           (error "inconsistent input to gauss-elim.")))
     lhs-rows)
    (dofor i
           1
           n
           (dofor j
                  1
                  n
                  (setf (aref ab j i)
                        (ith i (ith j lhs-rows)))))
    (dofor i
           1
           n
           (setf (aref ab i np1)
                 (ith i rhs-column)))
    (gauss-aux ab n np1)))



#|
;; example
(setq a '((45 66 23 89) 
          (72 27 14 53)
          (71 34 69 33)
          (55 15 72 41)))
(setq b '(2 4 6 8))
(setq c (gauss-elim a b))
(check-gauss-elim a c)

|#

________________________________________________________
Dimitri Simos, Lissys Ltd.  Tel&Fax: UK (+44) 509 235620