# Re: MCL linear systems

• To: info-mcl@cambridge.apple.com
• Subject: Re: MCL linear systems
• From: Dimitri Simos <dim@lissys.demon.co.uk>
• Date: Tue, 7 Jun 94 05:19:22 -0400

```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.,

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
(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
(list (aref ab i np)))))

(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

```