; 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)