Arrays proposal: Order of operations in array-indexes-set! and array-slice-set!?

General
1

SRFI 231 has array-assign!, which seems similar to math/array's array-indexes-set! and array-slice-set!. I'm wondering whether the specification of array-assign! should be changed to determine the order in which elements are assigned, and I'm studying array-indexes-set! and array-slice-set! for ideas.

Common example algorithms for array libraries are the Jacobi and Gauss-Seidel iterations for approximately solving Laplace's equation, see, e.g., Gauss–Seidel method - Wikipedia . The basic iteration is the same for both:

(define (jacobi-iteration! OUT IN)
  ;; If IN and OUT are different, this is a Jacobi iteration for Laplace's equation.
  ;; If IN and OUT are the same, this is a Gauss-Seidel iteration.
  (let ((interior (interval-dilate (array-domain IN) '#(1 1)  '#(-1 -1))))
    (array-assign!
     (array-extract OUT interior)
     (apply array-map
            (lambda (left right up down)
              (fl* 0.25 (fl+ left right up down)))
            (map (lambda (translation)
                   (array-extract (array-translate IN translation) interior))
                 '(#(0 -1) #(0 1) #(-1 0) #(1 0)))))))

For Jacobi's method, the order that elements are computed and assigned to the interior of OUT doesn't matter, because IN and OUT don't share elements, but for the Gauss-Seidel method the ordering makes a difference.

Often, algorithms assume that elements are computed and assigned to the interior of OUT in lexicographical (row-major) order, but that isn't (yet) guaranteed by array-assign!. For example, for so-called symmetric Gauss-Seidel, each iteration reverses the order in which you compute new elements, as in

(define (gauss-seidel-2d X iterations tolerance exact-answer #!optional (symmetric? #t))
  (let* ((X_0 (array-copy X)))
    (let loop ((n 0)
               (X_k X_0))
      (jacobi-iteration! X_k X_k)                   ;; forward gauss-seidel
      (if symmetric?
          (let ((reversed (array-reverse X_k)))
            (jacobi-iteration! reversed reversed))) ;; backward gauss-seidel
      (if (and (fxzero? (fxremainder n 10))
               (or (fx>= n iterations)
                   (fl< (frobenius exact-answer X_k) tolerance)))
          (begin
            (pretty-print (list gauss-seidel-iterations: n error: (frobenius exact-answer X_k) symmetric?: symmetric?))
            X_k)
          (loop (fx+ n 1) X_k)))))

In searching the math/array documentation, I see for array-indexes-set! and array-slice-set! (which appear to be similar to SRFI 231's array-assign!, please correct me if I'm wrong) the only mention of order of computation/assignment is something like

So I'm wondering whether the order of assignment is specified by the two math/array routines.

As for array-assign! in the SRFI followup, I can see the following considerations:

  1. Not specifying the order in which elements are assigned to the destination array allows some hypothetical theoretical future implementation to be faster, by exploiting parallelism, the run-time layout of arrays, etc.
  2. Specifying the order in which elements are assigned to the destination array simplifies the specification of some "bulk" array algorithms, where you don't have to drop down to index loops of the underlying interval operations to achieve your goals.

The sample implementation is slow compared to hand-written loops in C or Fortran, and I think it's unlikely that anyone will write a high-performance implementation in, e.g., the style of PetaLisp GitHub - marcoheisig/Petalisp: Elegant High-Performance Computing (which I admire very much, but it appears that Marco has gone on to other things).

So the question seems to be: Should the library be designed to allow precise specification of algorithms without dropping down into loops (or the low-level interval-for-each interface), limiting future hypothetical high-performance implementations, or should it have imprecise specifications of the order of some computations, which makes specifying some algorithms more difficult, but allows more flexibility in implementations?

Brad

1 Like
2

After a few day's thought, it's clear that yes, array-assign! should, by default, assign elements one by one, in lexicographical order, or one couldn't reason at all about algorithms operating on overlapping arrays.

For example, when moving a one-dimensional subarray of an array with lower bound zero to another subarray, one could write (clumsily) :

(define S_1 (make-specialized-array-from-data (string-copy "0123456789ABCDEF") char-storage-class #t))
(define S_2 (make-specialized-array-from-data (string-copy "0123456789ABCDEF") char-storage-class #t))

(define (non-overlapping-move! i  ;; index of left-most char in substring
                               j  ;; signed offset of how far to move substring
                               n  ;; number of characters to move
                               A) ;; array of characters

  (define (negative-case i j n A)
    (array-assign! (array-translate (array-extract A (make-interval (vector (+ i j))
                                                                    (vector (+ i j n))))
                                    (vector (- j)))
                   (array-extract A (make-interval (vector i) (vector (+ i n))))))
  
  (cond ((negative? j)
         (negative-case i j n A))
        ((positive? j)
         (let ((L (interval-upper-bound (array-domain A) 0))
               (reversed-A (array-reverse A)))
           (negative-case (- L (+ i n)) (- j) n reversed-A)))
        (else
         A)))

(pretty-print (array-body S_1))
(non-overlapping-move! 5 -3 5 S_1)
(pretty-print (array-body S_1))
(newline)
(pretty-print (array-body S_2))
(non-overlapping-move! 5 3 5 S_2)
(pretty-print (array-body S_2))

which gives the result

"0123456789ABCDEF"
"0156789789ABCDEF"

"0123456789ABCDEF"
"0123456756789DEF"

but only if array-assign! moves elements in a fixed, row-major, order.

Perhaps I'll add an option to allow moving elements out of order, but the sample implementation won't take it into account.

1 Like