Matrix API and layouts

Introduction

• APIs are implementations of algebras. The domain and operations of algebras should be reflected in the API but also the fundamental properties of the underlying hardware.
• Entities and operations:
  • Base-256 numbers, to 8 digits wide.
  • Contiguous blocks of equal-sized base-256 numbers, from 16 digits to 64 digits wide in total.
  • Registers and memory.
  • Memory addresses 4 to 8 bytes wide.
  • Memory offsets 4 bits to 8 bytes wide.

Hardware limitations

• Hardware properties:
  • No registers or stack for arbitrary size data. CPUs are fundamentally 3-operand (2 sources, 1 destination).
  • Random memory reads possible but very expensive.
  • Writes and random writes even more expensive.
  • Except for rather small (< 64KiB) contiguous blocks (caching)
  • Vectors of 4 or 8 elements hugely favoured.
  • Arithmetic 0.1-10 cycles, cache access 1-100 cycles, RAM access 1,000-10,000 cycles, SSD access 1,000,000-10,000,000 cycles, HDD access...
  • Write accesses to RAM and SSD much slower than reads.

Styles of operations

Side note on 3 styles of operations:

• Pure imperative (mul(&t1,b,c); mul(&t2,d,e); add(&a,t1,t2)): all operations do not return any values and all operands and results are parameters (in, out or change).
• Pure arithmetic (a = (b*c)+(d*e)): all operations return values and do not modify any operands. Much worse than pure imperative both for documentation and for efficiency: CPUs are not at all arithmetic, they are imperative and handle only limited size primitive types, even if there are exceptions. Also programming is not mathematics, not even applied mathematics.
• Both (addto(mul(&a,b,c),mul(&t1,d,e))): operations where possible both have in, out, change parameters and return values or references as appropriate. I prefer this.

API outlines

• Really important to have:
  • opaque types where implementations are hidden.
  • To have as many operands as possible.
  • To return a value or a reference that can be used as a further operand.
  • All allocations done if possible before the call (else use a mark-release separate stack allocator).
• Typical signatures (static or dynamic dispatching usually elides the _impl):
  • [v = | (void)] ro ref type_impl_GET[_field](&m,i, &v,j,k)
  • [m = | (void)] rw ref type_impl_PUT[_field](&m,i, v,j,k)
  • [m = | (void)] rw ref type_impl_OP[_field](&m, v)
  • [m = | (void)] rw ref type_impl_OP_N[_field](&m, n,v1,v2,...)
  • [m = | (void)] rw ref type_impl_OP_type_impl(&m, v)
  • [m = | (void)] rw ref type_impl_OP_N_type_impl(&m, n,v1,v2,...)
  • [m = | (void)] rw ref type_impl_FROM_type_impl(&m, v)
• Example:

  	    ref matrix matrix_rc_pow(ref matrix_rc m,const natural n)
  	    {
  	      if n > 1:
  		{ matrix mt; void) matrix_rc_transpose(mt,m); }
  	      for 2..n:
  		{ (void) matrix_rc_multiply_inner_transposed(m,mt); }
  	      return m;
  	    }
  	  

Layout fundamental choices

• Physical order:
  • An array of records
  • A record of arrays
  Like dict and list in Python.
• Space allocation:
  • Fixed size array
  • Resizable array
  • Resizable array with gap (usually only at the end).

Dense or sparse

Big question: dense or sparse?

• Random accesses very expensive, so sparse must be very sparse. Fortunately some important matrix types are very sparse.
• Most operations among sparse matrices result in dense matrices. Fortunately operations with sparse matrices and scalars result in sparse matrices.
• What about matrices that cannot fit in cache, memory or even one SSD or HDD? They are still in effect "sparse" in some sense.

Layout examples

• Ordering in memory:
  • Array of pointers to rows
  • Array of pointers to columns
  • Single array indexed row-wise
  • Single array indexed column-wise
  • Single array indexed submatrix-wise
  • Single array 4-8 wide row-strips
  • Single array 4-8 wide column-strips
  • Single array with pointer to closure of index procedure

  Pointers can mean relative indices.

• Ordering of elements:
  • Single array indexed by diagonal(s)
  • 3-diagonals array(s)
  • (Potentially recursive) tree of sub-matrices
  • Recursive tree of quadrants
  • Associative table key being indices.

Consequences

• Fitting subsets in cache, memory or SSD can give enormous speedups.
• If that is possible then pipeling can give huge additional speedups and to achieve that instead of doing conditionals inside loops it is better to be conditionals on entire loops if possible (using SIMD is a special short case). Arithmetic tricks can be used to avoid conditionals.
• In a game (3D triangles) application switching from row layout to columnn layout gave an overall speedup of 10%:
  • That is enormous for a single simple change.
  • The change was simple because all operations were shaped as in the above, so no applications code needed changing (really! totaly opaque matrix types).