home *** CD-ROM | disk | FTP | other *** search
/ Eagles Nest BBS 8 / Eagles_Nest_Mac_Collection_Disc_8.TOAST / Developer Tools⁄Additions / MacScheme20 / Mathlib / matrix.scm < prev    next >
Encoding:
Text File  |  1989-04-27  |  6.2 KB  |  232 lines  |  [TEXT/????]
  1. ;;; $Header: matrix.scm,v 1.2 87/08/28 00:57:38 GMT gjs Exp $
  2. ;;;; Matrix Package
  3.  
  4. (if-mit
  5.  (declare (usual-integrations = + - * /
  6.                  zero? 1+ -1+
  7.                  ;; truncate round floor ceiling
  8.                  sqrt exp log sin cos)))
  9.  
  10. ;;; A matrix is represented as a vector of rows, where each row
  11. ;;; is a vector of field elements.
  12.  
  13. (define ((matrix-elementwise f) . matrices)
  14.   (generate-matrix
  15.     (num-rows (car matrices))
  16.     (num-cols (car matrices))
  17.     (lambda (i j)
  18.       (apply f (map (lambda (m) (matrix-ref m i j)) matrices)))))
  19.  
  20. (define (make-matrix n m . optional)
  21.   (let ((ans (make-vector n))
  22.         (filler (if optional (car optional) '())))
  23.     (let loop ((i 0))
  24.       (if (= i n)
  25.       ans
  26.       (begin (vector-set! ans i
  27.                               (make-vector m filler))
  28.          (loop (+ i 1)))))))
  29.  
  30. (define (matrix-by-rows M)
  31.   (apply vector (map list->vector M)))
  32.  
  33. (define (matrix-by-cols M)
  34.   (transpose (matrix-by-rows M)))
  35.  
  36. (define (generate-matrix rows cols proc)
  37.   (let ((ans (make-matrix rows cols)))
  38.     (let loop-rows ((i 0))
  39.       (if (= i rows)
  40.       ans
  41.       (begin (let loop-columns ((j 0))
  42.            (if (= j cols)
  43.                'done
  44.                (begin (matrix-set! ans i j (proc i j))
  45.                   (loop-columns (+ j 1)))))
  46.          (loop-rows (+ i 1)))))))
  47.  
  48. (define (num-rows m) (vector-length m))
  49. (define (num-cols m) (vector-length (vector-ref m 0)))
  50.  
  51. (define (matrix-ref m i j)
  52.   (vector-ref (vector-ref m i) j))
  53.  
  54. (define (matrix-set! m i j v)
  55.   (vector-set! (vector-ref m i) j v))
  56.  
  57. (define (vector->row-matrix v)
  58.   (vector v))
  59.  
  60. (define (vector->col-matrix v)
  61.   (transpose (vector->row-matrix v)))
  62.  
  63. (define (outer-product v1 v2)
  64.   (matrix*matrix (vector->col-matrix v1)
  65.                  (vector->row-matrix v2)))
  66.  
  67. (define (nth-row M n) ;return as a vector
  68.   (vector-ref M n))
  69.  
  70. (define (nth-col M n)
  71.   (nth-row (transpose M) n))
  72.  
  73. (define (copy-matrix m)
  74.   (generate-matrix (num-rows m) (num-cols m)
  75.     (lambda (i j) (matrix-ref m i j))))
  76.  
  77.  
  78. (define add-matrices
  79.   (vector-elementwise add-vectors))
  80.  
  81. (define sub-matrices
  82.   (vector-elementwise sub-vectors))
  83.  
  84. (define (mul-matrices . m)
  85.   (let loop ((remaining m) (prod (identity-matrix (num-rows (car m)))))
  86.     (if (null? remaining)
  87.         prod
  88.         (loop (cdr remaining) (matrix*matrix prod (car remaining))))))
  89.  
  90. (define (matrix*matrix m1 m2)
  91.   (let ((m1r (num-rows m1))
  92.     (m1c (num-cols m1))
  93.     (m2r (num-rows m2))
  94.     (m2c (num-cols m2)))
  95.     (if (not (and (= m1c m2r)))
  96.     (error "matrix sizes do not match -- MATRIX*MATRIX" m1 m2)
  97.     (let ((tm2 (transpose m2)))
  98.       (generate-matrix m1r
  99.                            m2c
  100.                        (lambda (i j)
  101.                              (dot-product (vector-ref m1 i)
  102.                       (vector-ref tm2 j))))))))
  103.  
  104.  
  105. (define (transpose m)
  106.   (generate-matrix (num-cols m)
  107.                    (num-rows m)
  108.                    (lambda (i j) (matrix-ref m j i))))
  109.  
  110. (define (matrix*scalar m k)
  111.   (generate-matrix (num-rows m) 
  112.                    (num-cols m)
  113.                    (lambda (i j) (* k (matrix-ref m i j)))))
  114.  
  115. (define (scale-matrix k)
  116.   (lambda (m) (matrix*scalar m k)))
  117.  
  118. (define matrix*vector
  119.   (lambda (mat vec)
  120.     (if (not (= (vector-length vec) (num-cols mat)))
  121.         (error "Not conformable for multiplication" mat vec)
  122.         (generate-vector (num-rows mat)
  123.                          (lambda (i)
  124.                            (dot-product (nth-row mat i) vec))))))
  125.  
  126. (define vector*matrix
  127.   (lambda (vec mat)
  128.     (if (not (= (vector-length vec) (num-rows mat)))
  129.         (error "Not conformable for multiplication" vec mat)
  130.         (generate-vector (num-cols mat)
  131.                          (lambda (i)
  132.                            (dot-product vec (nth-col mat i)))))))
  133.  
  134.  
  135. (define (zero-matrix n)
  136.   (make-vector n (make-vector n 0)))
  137.  
  138. (define (identity-matrix n)
  139.   (generate-vector
  140.     n
  141.     (lambda (i) (unit-vector n i))))
  142.  
  143. (define (general-determinant zero add sub mul)
  144.   (define (internal-determinant m)
  145.     (let ((size (num-cols m)))
  146.       (if (= size 1)
  147.       (matrix-ref m 0 0)
  148.       (let colloop ((i 0) (ans zero) (sgn 1))
  149.         (if (= i size)
  150.         ans
  151.         (colloop (1+ i)
  152.              (add ans
  153.                               (mul sgn
  154.                    (matrix-ref m 0 i)
  155.                    (internal-determinant (cofactor m 0 i))))
  156.                (sub zero sgn)))))))
  157.   internal-determinant)
  158.  
  159. (define (cofactor m i j)
  160.   (generate-matrix
  161.    (-1+ (num-rows m))
  162.    (-1+ (num-cols m))
  163.    (lambda (a b)
  164.      (matrix-ref
  165.       m
  166.       (if (< a i) a (1+ a))
  167.       (if (< b j) b (1+ b))))))
  168.  
  169. (define determinant (general-determinant 0 + - *))
  170.  
  171. (define (det M)
  172.   (let ((lu (lu-decompose M allow-zero-pivot))
  173.         (n (num-rows M)))
  174.     (let ((a (car lu))  ; the matrix part
  175.           (p (cadr lu)) ; the permutation, as a vector
  176.       (s (caddr lu))) ; sign of permutation
  177.       (if (< (lu-rank a) n)
  178.           0
  179.           (let loop ((i 0) (prod s))
  180.             (if (= i n)
  181.                 prod
  182.                 (loop (+ i 1) (* prod (matrix-ref a i i)))))))))
  183.  
  184. (define (diagonal m)
  185.   (if (not (= (num-rows m) (num-cols m)))
  186.       (error "diagonal of a non-square matrix" m))
  187.   (let loop ((i 0))
  188.     (if (= i (num-rows m))
  189.         '()
  190.         (cons (matrix-ref m i i)
  191.               (loop (1+ i))))))
  192.  
  193. (define (make-diagonal-matrix diag)
  194.   (let ((n (vector-length diag)))
  195.     (generate-matrix n n
  196.                      (lambda (i j)
  197.                        (if (not (= i j))
  198.                            0
  199.                            (vector-ref diag i))))))
  200.  
  201. (define (matrix-trace m)
  202.   (apply + (diagonal m)))
  203.  
  204.  
  205. (define (matinv a)
  206.   (let ((n (num-rows a))
  207.         (lu (lu-decompose a)))
  208.     (let ((lumat (car lu)) (luperm (cadr lu)))
  209.       (let ((m (generate-vector
  210.                  n
  211.                  (lambda (i)
  212.                    (let ((e (unit-vector n i)))
  213.                      (lu-backsubstitute lumat luperm e))))))
  214.         (transpose m)))))
  215.  
  216. (define (matrix^nth M p)
  217.   (let ((n (round p)))
  218.     (cond ((not (= n p))
  219.            (error "Only integer powers allowed -- MATRIX^NTH"))
  220.           ((< n 0) 
  221.            (matrix^nth (matinv M) (- n)))
  222.           ((zero? n)
  223.            (identity-matrix (num-rows M)))
  224.           (else
  225.            (let loop ((count n))
  226.              (cond ((= count 1) M)
  227.                    ((even? count) 
  228.                     (let ((a (loop (quotient count 2))))
  229.                       (matrix*matrix a a)))
  230.            (else
  231.             (matrix*matrix M (loop (-1+ count))))))))))
  232.