home *** CD-ROM | disk | FTP | other *** search
/ Eagles Nest BBS 8 / Eagles_Nest_Mac_Collection_Disc_8.TOAST / Developer Tools⁄Additions / MacScheme20 / Mathlib / linear.scm < prev    next >
Encoding:
Text File  |  1989-04-27  |  8.8 KB  |  286 lines  |  [TEXT/????]
  1. ;;; $Header: linear.scm,v 1.2 87/08/28 00:55:51 GMT gjs Exp $
  2. ;;;; 6.003 Linear System Solver(s)
  3.  
  4. ;;; This file contains the following definitions:
  5. ;;;
  6. ;;; (solve-linear-system a-matrix b-vector) => x-vector
  7. ;;;
  8. ;;; (lu-backsubstitute lu-matrix lu-permutation b-vector) => x-vector
  9. ;;;
  10. ;;; (lu-null-vector <lu-matrix> <index>) ==> <vector>
  11. ;;;
  12. ;;; (lu-decompose matrix [barf-if-zero-pivot]) =>
  13. ;;;                             (lu-matrix  lu-permutation  sign)
  14. ;;;
  15. ;;; (lu-rank lu-matrix) => non-negative-integer
  16. ;;; (lu-nullity lu-matrix) => non-negative-integer
  17. ;;;
  18. ;;; (sign-of-permutation lu-permutation) => {+1, -1}
  19. ;;;
  20.  
  21. (if-mit
  22.  (declare (usual-integrations = + - * /
  23.                  zero? 1+ -1+
  24.                  ;; truncate round floor ceiling
  25.                  sqrt exp log sin cos)))
  26.  
  27. ;;; Solves an inhomogeneous system of linear equations, A*X=B, returning
  28. ;;; the vector X.
  29.  
  30. (define (solve-linear-system a b)
  31.   (let ((lu (lu-decompose a)))
  32.     (let ((lumatrix (car lu))
  33.       (lupermutations (cadr lu)))
  34.       (lu-backsubstitute lumatrix lupermutations b))))
  35.  
  36.  
  37. ;;; THE NEXT PROCEDURE ASSUMES THAT THE MATRIX IS IN LU FORM
  38. ;;; The rank is just the number of non-zero diagonal elements.
  39.  
  40. (define (lu-nullity m)
  41.   (count-elements heuristically-zero?
  42.                   (diagonal m)))
  43.  
  44. (define (lu-rank m)
  45.   (let ((n (lu-nullity m)))
  46.     (- (num-rows m) n)))         
  47.  
  48. (define (heuristically-zero? z)
  49.   (< (magnitude z) heuristic-zero-test-bugger-factor))
  50.  
  51. ;;; The following number was assigned by HAL. -- GJS & MH.
  52. (define heuristic-zero-test-bugger-factor 1e-9)
  53.  
  54. ;;; LU decomposition -- uses Crout's algorithm (see Press, section 2.3).
  55. ;;; The LU decomposition returns a LIST of the resulting matrix, the
  56. ;;; permutation vector, and the sign (+1 or -1) of the permutation.
  57. ;;; The only currently-defined flag controls error-trapping on detection
  58. ;;; of a zero-pivot. This flag, if omitted, defaults to true.
  59.  
  60. ;;; We define the flags first:
  61. (define-integrable barf-on-zero-pivot true)
  62. (define-integrable allow-zero-pivot false)
  63.  
  64. (define (lu-decompose m . flags)
  65.   (let* ((barf-if-zero-pivot? (if flags (car flags) true))
  66.          (n (num-rows m))
  67.      (perms (generate-vector n (lambda (i) i)))  ;row permutations
  68.          (sign 1)  ;1 for even permutations, -1 for odd
  69.          ;; We must copy the matrix, since Crout's algorithm clobbers it.
  70.          (m (copy-matrix m)))
  71.     (let jloop ((j 0))
  72.     (if (< j n)
  73.         (begin
  74.           (let iloop ((i 0)) ;compute elements above diagonal
  75.         (if (<= i j)
  76.             (begin (matrix-set! m i j (upper-eqn i j m))
  77.                (iloop (1+ i)))))
  78.           (let iloop ((i (1+ j))) ;compute elements below diagonal
  79.         (if (< i n)
  80.             (begin (matrix-set! m i j (lower-eqn i j m))
  81.                (iloop (1+ i)))))
  82.           (let* ((pivot-info (find-best-pivot m j n))
  83.                      (pivot (car pivot-info))
  84.              (pivot-index (cdr pivot-info)))
  85.         (if (< (abs pivot) tiny-pivot-bugger-factor)
  86.                     (if barf-if-zero-pivot?
  87.                 (error "Singular matrix - zero pivot")
  88.                 (jloop (1+ j)))
  89.             (let ((1/pivot (/ 1 pivot)))
  90.               (lu-row-swap m j pivot-index perms)
  91.               (if (not (= j pivot-index))
  92.                   (set! sign (- sign)))
  93.               (let iloop ((i (1+ j))) ;divide through by pivot
  94.             (if (= i n)
  95.                 'done
  96.                 (begin (matrix-set! m i j
  97.                     (* (matrix-ref m i j) 1/pivot))
  98.                    (iloop (1+ i)))))
  99.               (jloop (1+ j))))))))
  100.         (list m perms sign)))
  101.  
  102. ;;; As a practical matter, the following theoretical factor need not
  103. ;;; be adjusted if conditions of STP do not precisely apply.
  104. (define tiny-pivot-bugger-factor (/ 1 6.0238e23))
  105.  
  106. ;;; Note that we start looping at the jth row in searching for the
  107. ;;; next pivot.  We can make this somewhat better numerically by
  108. ;;; including scale factors for the columns, but i'm too lazy to do
  109. ;;; this now. -- HAL
  110.  
  111. (define (find-best-pivot m j n)
  112.   (let ((column (generate-vector n (lambda (i) (matrix-ref m i j)))))
  113.     (let iloop ((i j) (bestindex -1) (maxval -1) (pivot -1))
  114.       (if (= i n)
  115.           (cons pivot bestindex)
  116.       (let* ((p (vector-ref column i))
  117.          (absp (abs p)))
  118.         (if (> absp maxval)
  119.         (iloop (+ i 1) i absp p)
  120.         (iloop (+ i 1) bestindex maxval pivot)))))))
  121.  
  122. (define (upper-eqn i j m)
  123.   (if (= i 0)
  124.       (matrix-ref m i j)
  125.       (let kloop ((k 0) (sum 0))
  126.     (if (= k i)
  127.         (- (matrix-ref m i j) sum)
  128.         (kloop (+ k 1)
  129.            (+ sum (* (matrix-ref m i k) (matrix-ref m k j))))))))
  130.  
  131. (define (lower-eqn i j m)
  132.   (let kloop ((k 0) (sum 0))
  133.     (if (= k j)
  134.     (- (matrix-ref m i j) sum)
  135.     (kloop (+ k 1)
  136.            (+ sum (* (matrix-ref m i k) (matrix-ref m k j)))))))
  137.  
  138. (define (lu-row-swap  m i1 i2 perms)
  139.   (define (swap-elements vector i j)
  140.     (let ((temp (vector-ref vector i)))
  141.       (vector-set! vector i (vector-ref vector j))
  142.       (vector-set! vector j temp)))
  143.   (if (not (= i1 i2))
  144.       (begin (swap-elements perms i1 i2)
  145.          ;;uses fact that matrix is a vector of rows
  146.          (swap-elements m i1 i2))))
  147.  
  148. ;;; Back substitution (see Press, page 32)
  149.  
  150. (define (lu-backsubstitute m perm b)
  151.   (let* ((n (vector-length b))
  152.      (top (- n 1))
  153.      (y (make-vector n '()))
  154.      (x (make-vector n '())))
  155.     (let fdloop ((i 0))
  156.       (if (< i n)
  157.       (begin
  158.         (vector-set! y i
  159.           (- (vector-ref b (vector-ref perm i))
  160.          (let jloop ((j 0) (sum 0))
  161.            (if (= j i)
  162.                sum
  163.                (jloop (+ 1 j)
  164.                   (+ sum
  165.                  (* (vector-ref y j)
  166.                     (matrix-ref m i j))))))))
  167.           (fdloop (+ i 1)))))
  168.     (let bkloop ((i top))
  169.       (if (>= i 0)
  170.       (begin
  171.         (vector-set! x i
  172.           (/ (- (vector-ref y i)
  173.             (let jloop ((j (+ i 1)) (sum 0))
  174.               (if (= j n)
  175.               sum
  176.               (jloop (+ 1 j)
  177.                  (+ sum
  178.                     (* (vector-ref x j)
  179.                        (matrix-ref m i j)))))))
  180.          (matrix-ref m i i)))
  181.         (bkloop (- i 1)))))
  182.     x))
  183.  
  184. ;;; In the case of a homogeneous system we can solve if the matrix is 
  185. ;;;  singular.  For each non-negative integer index less than the nullity of
  186. ;;;  the lu-matrix we can obtain an independent solution vector (the entire
  187. ;;;  set spans the null space).
  188.  
  189. ;;; In this case, the upper part of the matrix contains the U part of
  190. ;;; the LU decomposition.  The L part is irrelevant to this
  191. ;;; process since we are solving a homogeneous equation.  The
  192. ;;; algorithm is to work backwards, from n-1 to 0, solving the
  193. ;;; equations  m[i,i]x[i] + m[i,i+1]x[i+1] + ... = 0.
  194. ;;; If m[i,i] is zero (to within the tolerance), then we can set x[i]
  195. ;;; to an arbitrary value.
  196. ;;; We set the kth arbitary unknown to 1 and the rest to 0 (0<=k<nullity).
  197. ;;; If a luser supplies a k>=nullity, he will get a zero vector.
  198.  
  199. (define (lu-null-vector m k)
  200.   (let* ((n (num-rows m))
  201.      (top (- n 1))
  202.      (x (make-vector n '())))
  203.     (let bkloop ((i top) (acount 0))
  204.       (if (>= i 0)
  205.       (let ((p (matrix-ref m i i)))
  206.         (if (heuristically-zero? p)
  207.                 (begin
  208.                   (if (= acount k)
  209.               (vector-set! x i 1)
  210.                       (vector-set! x i 0))
  211.               (bkloop (-1+ i) (1+ acount)))
  212.         (let ((s (let jloop ((j (+ i 1)) (sum 0))
  213.                (if (= j n)
  214.                    sum
  215.                    (jloop (+ 1 j)
  216.                       (+ sum
  217.                      (* (vector-ref x j)
  218.                         (matrix-ref m i j))))))))
  219.           (vector-set! x i (/ (- s) p))
  220.           (bkloop (- i 1) acount))))))
  221.     x))
  222.  
  223. ;;; The following procedure is a hack to get the sign of the permutation
  224. ;;;  generated by LU-DECOMPOSE.  This is unnecessary now, because LU-DECOMPOSE
  225. ;;;  has since been modified to return the permutation sign directly.
  226.  
  227. (define (sign-of-permutation perm) ;perm = #(p1 p2 ... pn)
  228.   (define (make-ordered-pairs L)   ;L is a list, not a vector
  229.     (define (pair-first-with-others L)
  230.       (if (< (length L) 2)
  231.           '()
  232.           (let ((a (car L)) (b (cdr L)))
  233.             (map (lambda (u) (cons a u)) b))))
  234.     (if (< (length L) 2)
  235.         '()
  236.         (append (pair-first-with-others L)
  237.                 (make-ordered-pairs (cdr L)))))
  238.   (let ((ordered-pairs (make-ordered-pairs (vector->list perm))))
  239.     (let loop ((L ordered-pairs) (inversion-count 0))
  240.       (if (null? L)
  241.           (if (even? inversion-count) 1 -1)
  242.           (let ((pair (car L)))
  243.             (if (> (car pair) (cdr pair))
  244.                 (loop (cdr L) (+ inversion-count 1))
  245.                 (loop (cdr L) inversion-count)))))))
  246.  
  247. ;;; This is stuff to test the LU-decomposition
  248.  
  249. (define (hilb n)
  250.   (generate-matrix n n (lambda (i j) (/ 1 (+ i j 2)))))
  251.  
  252. (define (upper-matrix m)
  253.   (let ((n (num-rows m)))
  254.     (generate-matrix n n
  255.       (lambda (i j)
  256.     (if (> i j)
  257.         0
  258.         (matrix-ref m i j))))))
  259.  
  260. (define (lower-matrix m)
  261.   (let ((n (num-rows m)))
  262.     (generate-matrix n n
  263.       (lambda (i j)
  264.     (cond ((< i j) 0)
  265.           ((= i j) 1)
  266.           (else (matrix-ref m i j)))))))
  267.  
  268.  
  269. (define (lower*upper m)
  270.   (matrix*matrix (lower-matrix m)
  271.          (upper-matrix m)))
  272.  
  273.  
  274. (define (checklu m)
  275.   (let ((ans (lu-decompose m allow-zero-pivot)))
  276.     (let ((ans-mat (car ans)) (ans-perm (cadr ans)))
  277.       (let ((prod (lower*upper ans-mat)))
  278.     (let ((n (num-rows m)))
  279.       (generate-matrix n n
  280.         (lambda (i j)
  281.               (- (matrix-ref prod i j)
  282.              (matrix-ref m
  283.                  (vector-ref ans-perm i)
  284.                  j)))))))))
  285.  
  286.