Re: Matrix or array operations library

[Matt Wette <matt.wette@gmail.com>]

So I tried an example using Guile to perform Gaussian elimination
(in situ, no column pivoting so unstable) using shared-arrays and
using explicit computation of offsets.  The explicit approach is
much more grungy but I'm guessing it does a lot less memory allocation.

Also, is it possible to get a pointer to the first element of an array
(i.e., it's #f64() root vector), for the purpose of using with FFI?

(use-modules ((srfi srfi-1) #:select (fold)))
(use-modules (srfi srfi-11))

;; === method 1: use share arrays

(define (axpy1 a x y)
   (let ((n (car (array-dimensions x))))
     (if (not (= n (car (array-dimensions y)))) (throw 'error))
     (let loop ((i 0))
       (when (< i n)
	(array-set! y (+ (* a (array-ref x i)) (array-ref y i)) i)
	(loop (1+ i))))))

;; modify A, b in-situ so that A is upper triangular
(define (reduce1 A b n)
   (let ((n-1 (1- n)))
     (let loop ((kl (iota n)))
       (when (pair? kl)
	(let* ((k (car kl))
	       (alpha (- (/ 1.0 (array-ref A k k))))
	       (x (make-shared-array A (lambda (j) (list k j)) n)))
	  (let loop1 ((i (1+ k)))
	    (when (< i n)
	      (let ((a (* alpha (array-ref A i k)))
		    (y (make-shared-array A (lambda (j) (list i j)) n)))
		(axpy1 a x y)
		(array-set! b (+ (* a (array-ref b k)) (array-ref b i)) i))
	      (loop1 (1+ i)))))
	(loop (cdr kl))))))
   
;; solve A*x = b for x where A is upper triangular
(define (bksolv1 A b x n)
   (let loop ((k (1- n)))
     (unless (negative? k)
       (let loop1 ((r 0.0) (j (1- n)))
	(cond
	 ((< k j)
	  (loop1 (+ r (* (array-ref A k j) (array-ref x j))) (1- j)))
	 (else
	  (array-set! x (/ (- (array-ref b k) r) (array-ref A k k)) k))))
       (loop (1- k)))))


;; === method 2: use root vector and computing offsets explicitly

;; based on root vectors
;; x, y are root vectors (from shared-array-root)
;; x0, y0 are starting indices
;; xi, yi are increments
(define (axpy2 n a x x0 xi y y0 yi)
   (let loop ((ix x0) (iy y0) (i 0))
     (when (< i n)
       (array-set! y (+ (* a (array-ref x ix)) (array-ref y iy)) iy)
       (loop (+ ix xi) (+ iy yi) (1+ i)))))

;; @deffn {Procedure} array-vec-model a vec-idx . idxs
;; Generate params for accessing a vector in an array.
;; Returned expression is
;; (values rv off inc)
;; where @code{rv} is the root-vector, @code{off} is the offset of the
;; first element and @code{inc} is the increment.
;; @end deffn
(define (array-vec-model a vec-idx . idxs)
   (let ((ic (shared-array-increments a)))
     (values (shared-array-root a)
	    (fold (lambda (bd ic ix of) (+ of (* ic (- ix (car bd)))))
		  (shared-array-offset a) (array-shape a) ic idxs)
	    (list-ref ic vec-idx))))

;; @deffn {Procedure} array-mat-model a row-idx col-idx . idxs
;; Generate params for accessing a vector in an array.
;; Returned expression is
;; (values rv off inc)
;; where @code{rv} is the root-vector, @code{off} is the offset of the
;; first element and @code{inc} is the increment.
;; @end deffn
(define (array-mat-model a row-idx col-idx . idxs)
   (let ((ic (shared-array-increments a)))
     (values (shared-array-root a)
	    (fold (lambda (bd ic ix of) (+ of (* ic (- ix (car bd)))))
		(shared-array-offset a) (array-shape a) ic idxs)
	    (list-ref ic row-idx)
	    (list-ref ic col-idx))))

;; modify A, b in-situ so that A is upper triangular
(define (reduce2 A b n)
   (let-values
       (((av a0 arinc acinc) (array-mat-model A 0 1 0 0))
        ((bv b0 binc) (array-vec-model b 0 0)))
     (let loop ((kl (iota n)) (akk a0) (bk b0))
       (when (pair? kl)
	(let* ((k (car kl))
	       (alpha (- (/ 1.0 (array-ref av akk)))))
	  (let loop1 ((i (1+ k)) (aik (+ akk arinc)) (bi (+ bk binc)))
	    (when (< i n)
	      (let ((a (* alpha (array-ref av aik))))
		(axpy2 (- n k) a av akk acinc av aik acinc)
		(array-set! bv (+ (* a (array-ref bv bk)) (array-ref bv bi)) bi))
	      (loop1 (1+ i) (+ aik arinc) (+ bi binc)))))
	(loop (cdr kl) (+ akk arinc acinc) (+ bk binc))))))

;; solve A*x = b for x where A is upper triangular
(define (bksolv2 A b x n)
   (let*-values
       (((av a00 arinc acinc) (array-mat-model A 0 1 0 0))
        ((bv b0 binc) (array-vec-model b 0 0))
        ((xv x0 xinc) (array-vec-model x 0 0))
        ((bn) (+ b0 (* (1- n) binc))))
     (let loop ((k (1- n))
	       (akn (+ a00 (* (1- n) arinc) (* (1- n) acinc)))
	       (xk (+ x0 (* (1- n) xinc))))
       (unless (negative? k)
	(let loop1 ((r 0.0) (j (1- n)) (akj akn) (bj bn))
	  (cond
	   ((< k j)
	    (loop1 (+ r (* (array-ref av akj) (array-ref x j)))
		   (1- j) (- akj acinc) (- bj binc)))
	   (else
	    (array-set! xv (/ (- (array-ref b j) r) (array-ref av akj)) xk))))
	(loop (1- k) (- akn arinc) (- xk xinc))))))


;; === example problem

(define A
   (list->typed-array
    'f64 2 '((11.0 1.2 1.3 1.4)
	    (2.1 22.0 2.3 2.4)
	    (3.1 3.2 33.0 3.4)
	    (4.1 4.2 4.3 44.0))))
(define b
   (list->typed-array
    'f64 1 '(4.0 3.0 2.0 1.0)))

(define x (make-typed-array 'f64 0.0 4))

;; expected result
(define x-soln
   (list->typed-array
    'f64 1 '(0.352856 0.103010 0.019728 -0.021913)))

(display "expect:\n") (vec-disp x-soln) (newline)

(when #f
   (reduce1 A b 4)
   ;;(display "A:\n") (mat-disp A)
   ;;(display "b:\n") (vec-disp b)
   (bksolv1 A b x 4)
   (display "method1:\n") (vec-disp x)
   )

(when #t
   (reduce2 A b 4)
   ;;(display "A:\n") (mat-disp A)
   ;;(display "b:\n") (vec-disp b)
   (bksolv2 A b x 4)
   (display "method2:\n") (vec-disp x)
   )