Re: Matrix or array operations library

[Matt Wette]

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