codepad

Language:

; solving systems of linear equations

(define-syntax fold-of
  (syntax-rules (range in is)
    ((_ "z" f b e) (set! b (f b e)))
    ((_ "z" f b e (v range fst pst stp) c ...)
      (let* ((x fst) (p pst) (s stp)
             (le? (if (positive? s) <= >=)))
        (do ((v x (+ v s))) ((le? p v) b)
          (fold-of "z" f b e c ...))))
    ((_ "z" f b e (v range fst pst) c ...)
      (let* ((x fst) (p pst) (s (if (< x p) 1 -1)))
        (fold-of "z" f b e (v range x p s) c ...)))
    ((_ "z" f b e (v range pst) c ...)
      (fold-of "z" f b e (v range 0 pst) c ...))
    ((_ "z" f b e (x in xs) c ...)
      (do ((t xs (cdr t))) ((null? t) b)
        (let ((x (car t)))
          (fold-of "z" f b e c ...))))
    ((_ "z" f b e (x is y) c ...)
      (let ((x y)) (fold-of "z" f b e c ...)))
    ((_ "z" f b e p? c ...)
      (if p? (fold-of "z" f b e c ...)))
    ((_ f i e c ...)
      (let ((b i)) (fold-of "z" f b e c ...)))))

(define-syntax sum-of (syntax-rules ()
  ((_ arg ...) (fold-of + 0 arg ...))))

(define (make-matrix rows columns . value)
  (do ((m (make-vector rows)) (i 0 (+ i 1)))
      ((= i rows) m)
    (if (null? value)
        (vector-set! m i (make-vector columns))
        (vector-set! m i (make-vector columns (car value))))))

(define (matrix-rows x) (vector-length x))

(define (matrix-cols x) (vector-length (vector-ref x 0)))

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

(define (matrix-set! m i j x) (vector-set! (vector-ref m i) j x))

(define-syntax for
  (syntax-rules ()
    ((for (var first past step) body ...)
      (let ((ge? (if (< first past) >= <=)))
        (do ((var first (+ var step)))
            ((ge? var past))
          body ...)))
    ((for (var first past) body ...)
      (let* ((f first) (p past) (s (if (< first past) 1 -1)))
        (for (var f p s) body ...)))
    ((for (var past) body ...)
      (let* ((p past)) (for (var 0 p) body ...)))))

(define (lu-decomposition a)
  (let* ((n (matrix-rows a))
         (l (make-matrix n n 0))
         (u (make-matrix n n 0)))
    (for (i n) (matrix-set! l i i 1))
    (for (k n)
      (matrix-set! u k k (matrix-ref a k k))
      (for (i (+ k 1) n)
        (matrix-set! l i k (/ (matrix-ref a i k) (matrix-ref u k k)))
        (matrix-set! u k i (matrix-ref a k i)))
      (for (i (+ k 1) n)
        (for (j (+ k 1) n)
          (matrix-set! a i j
            (- (matrix-ref a i j)
               (* (matrix-ref l i k)
                  (matrix-ref u k j)))))))
    (values l u)))

(define a
  #( #( 2  3  1  5 )
     #( 6 13  5 19 )
     #( 2 19 10 23 )
     #( 4 10 11 31 )))

(call-with-values
  (lambda () (lu-decomposition a))
  (lambda (l u) (display l) (newline)
                (display u) (newline)))
(newline)

(define (vector-swap! v a b)
  (let ((t (vector-ref v a)))
    (vector-set! v a (vector-ref v b))
    (vector-set! v b t)))

(define (matrix-swap! m ar ac br bc)
  (let ((t (matrix-ref m ar ac)))
    (matrix-set! m ar ac (matrix-ref m br bc))
    (matrix-set! m br bc t)))

(define (lup-decomposition a)
  (let* ((n (matrix-rows a)) (pi (make-vector n 0)))
    (for (i n) (vector-set! pi i i))
    (for (k n)
      (let ((p 0) (k-prime 0))
        (for (i k n)
          (let ((x (abs (matrix-ref a i k))))
            (when (< p x) (set! p x) (set! k-prime i))))
        (when (zero? p)
          (error 'lup-decomposition "singular matrix"))
        (vector-swap! pi k k-prime)
        (for (i n) (matrix-swap! a k i k-prime i))
        (for (i (+ k 1) n)
          (matrix-set! a i k
                       (/ (matrix-ref a i k)
                          (matrix-ref a k k)))
          (for (j (+ k 1) n)
            (matrix-set! a i j
                         (- (matrix-ref a i j)
                            (* (matrix-ref a i k)
                               (matrix-ref a k j))))))))
    (values a pi)))

(define a
  #( #(  2    0    2  3/5 )
     #(  3    3    4   -2 )
     #(  5    5    4    2 )
     #( -1   -2 17/5   -1 )))

(call-with-values
  (lambda () (lup-decomposition a))
  (lambda (a pi) (display a) (newline)
                 (display pi) (newline)))
(newline)

(define (lower lu i j)
  (cond ((< i j) 0) ((= i j) 1)
    (else (matrix-ref lu i j))))

(define (upper lu i j)
  (if (< j i) 0
    (matrix-ref lu i j)))

(define (lup-solve lu pi b)
  (let* ((n (matrix-rows lu))
         (y (make-vector n))
         (x (make-vector n)))
    (for (i n)
      (vector-set! y i
                   (- (vector-ref b (vector-ref pi i))
                      (sum-of (* (lower lu i j) (vector-ref y j))
                        (j range 0 i)))))
    (for (i (- n 1) -1)
      (vector-set! x i
                   (/ (- (vector-ref y i)
                         (sum-of (* (upper lu i j) (vector-ref x j))
                           (j range (+ i 1) n)))
                      (upper lu i i))))
    x))

(define a
  #( #(1 2 0)
     #(3 5 4)
     #(5 6 3)))

(define b #(1/10 25/2 103/10))

(call-with-values
  (lambda () (lup-decomposition a))
  (lambda (lu p) (display (lup-solve lu p b))))
(newline)

Private [?]

Run code