guile-user
[Top][All Lists]

## Re: Matrix or array operations library

 From: Matt Wette Subject: Re: Matrix or array operations library Date: Sun, 27 Jan 2019 10:41:57 -0800 User-agent: Mozilla/5.0 (X11; Linux x86_64; rv:60.0) Gecko/20100101 Thunderbird/60.4.0

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

```

reply via email to