From: mj@cs.brown.edu (Mark Johnson)
Subject: Re: Numerical analysis routines, anyone?
Organization: Brown University Department of Computer Science
Date: Mon, 21 Dec 1992 15:46:53 GMT

This is certainly not complete, but it is at least a start.
I wrote this package for some statistical language processing
that I'm doing.  The matrix routines all use a generic
arithmetic interface because I hope at some stage to be
working with more complex numerical types (specifically,
a value and its first derivative).

Anyway, I hope this is useful,

Mark

;;; Matrix manipulation package
;;;
;;; requires a DOTIMES macro a la CommonLisp

;;; generic arithmetic interface

(define multiply *)
(define add +)
(define subtract -)
(define divide /)
(define zero 0)
(define =zero? zero?)
(define one 1)
(define (inverse x) (/ 1.0 x))

(define (sum n fn)
  (define (sum-helper i partial-sum)
    (if (zero? i)
        (add (fn 0) partial-sum)
        (sum-helper (1- i) (add (fn i) partial-sum))))
  (sum-helper (1- n) zero))
    
;;; generic matrix interface.  This version uses vectors of vectors,
;;; and uses the generic arithemtic above

(define (make-zero-matrix dim)
  (let ((mat (make-vector dim)))
    (dotimes (i dim)
      (vector-set! mat i (make-vector dim zero)))
    mat))

(define (make-matrix dim)
  (make-zero-matrix dim))

(define (matrix-ref mat i j)
  (vector-ref (vector-ref mat i) j))

(define (matrix-set! mat i j value)
  (vector-set! (vector-ref mat i) j value)
  mat)

(define (matrix-inc! mat i j inc)
  (matrix-set! mat i j (add (matrix-ref mat i j) inc)))

(define (matrix-dec! mat i j inc)
  (matrix-set! mat i j (subtract (matrix-ref mat i j) inc)))

(define (matrix-mult! mat i j inc)
  (matrix-set! mat i j (multiply (matrix-ref mat i j) inc)))

(define (matrix-div! mat i j inc)
  (matrix-set! mat i j (divide (matrix-ref mat i j)  inc)))

(define (matrix-dimension mat)
  (vector-length mat))

(define (make-unit-matrix dim)
  (let ((mat (make-zero-matrix dim)))
    (dotimes (i dim)
      (matrix-set! mat i i one))
    mat))

;(define (matrix-swap-elts mat i1 j1 i2 j2)
;  (let ((temp (matrix-ref mat i1 j1)))
;    (matrix-set! mat i1 j1 (matrix-ref mat i2 j2))
;    (matrix-set! mat i2 j2 temp)))

(define (matrix-swap-rows! mat row1 row2)
  	(let ((temp (vector-ref mat row1)))
     (vector-set! mat row1 (vector-ref mat row2))
     (vector-set! mat row2 temp)))

(define (copy-matrix mat)
  (let* ((n (matrix-dimension mat))
         (copy (make-matrix n)))
    (dotimes (i n)
      (dotimes (j n)
        (matrix-set! copy i j (matrix-ref mat i j))))
    copy))

;;; matrix multiplication

(define (matrix-multiply mat1 mat2)
  (let* ((n (matrix-dimension mat1))
         (ans (make-zero-matrix n)))
    (dotimes (i n)
      (dotimes (j n)
        (matrix-set! ans i j
                     (sum n (lambda (k)
                              (multiply (matrix-ref mat1 i k)
                                        (matrix-ref mat2 k j)))))))
    ans))

(define (matrix-subtract mat1 mat2)
  (let* ((n (matrix-dimension mat1))
         (ans (make-matrix n)))
    (dotimes (i n)
      (dotimes (j n)
        (matrix-set! ans i j (subtract (matrix-ref mat1 i j)
                                       (matrix-ref mat2 i j)))))
    ans))

;;; gauss-jordan does matrix inversion.  Note that this destructively changes 
;;; the argument matrix mat

(define (gauss-jordan mat)
  (let* ((n (matrix-dimension mat))
         (inv (make-unit-matrix n)))
    (dotimes (i n)
     	(do ((pi i (1+ pi)))
          ((or (= pi n)
               (not (=zero? (matrix-ref mat pi i))))
           (cond ((= pi n) (error "GAUSS-JORDAN: matrix to invert is singular"))
                 ((not (= i pi)) (matrix-swap-rows! mat i pi)
                                 (matrix-swap-rows! inv i pi)))))
      (let ((pivot (matrix-ref mat i i)))
        (dotimes (j n)
          (if (not (= j i))
            (let ((factor (divide (matrix-ref mat j i) pivot)))
              (matrix-set! mat j i zero)
              (do ((k (1+ i) (1+ k)))
                		((= k n))
             	  (matrix-dec! mat j k (multiply (matrix-ref mat i k) factor)))
              (dotimes (k n)
                (matrix-dec! inv j k (multiply (matrix-ref inv i k) factor))))))
        (matrix-set! mat i i one)
        (do ((j (1+ i) (1+ j)))
            ((= j n))
          (matrix-div! mat i j pivot))
        (dotimes (j n)
          (matrix-div! inv i j pivot))))
    inv))
            	  

Mark Johnson
Cognitive Science, Box 1978
Brown University
Xref: sparky comp.lang.scheme:2852 comp.lang.c:19132