From PackedArray To NumericArray

Mark Sofroniou,
Wolfram Research

Abstract

PackedArray is a data structure that permits efficient machine precision computations in Mathematica. This talk gives an outline of the design, scope and implementation of PackedArrays. NumericArray is a new generalization of PackedArray that allows computations for a wider set of machine precision types. An overview of the goals and differences with PackedArrays will be given, together with some new algorithms that have been implemented. Integration into the new compiler will also be described.

Initialization

SetAttributes[SerialParallelTiming, HoldFirst];SerialParallelTiming[expr_]:=
Module[{opts ,pc,ts, tp},
opts = SystemOptions["ParallelOptions"];
pc = $ProcessorCount;
SetSystemOptions["ParallelOptions"->{"MKLThreadNumber"->1,"ParallelThreadNumber"->1}];
ts = AbsoluteTiming[expr];
SetSystemOptions["ParallelOptions"->{"MKLThreadNumber"->pc,"ParallelThreadNumber"->pc}];
tp = AbsoluteTiming[expr];
SetSystemOptions[opts];
{ts, tp}
];
$HistoryLength = 0;
Needs["Compile`"];

PackedArray

Started life in Compile in V3.0 as a data structure MTensor.

PackedArrays were added to the kernel in V4.0 using an expression wrapper for the MTensor data structure.

MTensor has three basic types

MachineInteger : largest signed machine integer (different for 32 or 64 bit machines)

MachineReal : 64 bit floating point number (double)

MachineComplex : complex pair of real 64 bit floating point numbers (double complex in C, complex<double> in C++)

MTensor types match the native machine number types in Mathematica

PackedArray functionality refined and extended over numerous release cycles.

PackedArray Advantages

PackedArrays are practically indisinguishable from Lists

Many functions generate PackedArrays automatically

pa = Range[10]

list = {1,2,3,4,5,6,7,8,9,10}

Map[Developer`PackedArrayQ, {pa, list}]

PackedArrays have better data locality and use less memory.

PackedArrays : ({{, ⋯, }})             Lists : ({{, ⋯, }})

Arithmetic and structural operations are much more efficient.

Map[ByteCount, {pa, list}]

Convert Lists to PackedArrays when there are benefits

Autocompile to run more efficiently above a certain threshold

PackedArray Disadvantages

PackedArrays are practically indisinguishable from Lists

Can be difficult to track when operations are packed/unpacked

Type overflow can lead to a significant degradation in performance

Can be useful to work with an invariant fixed-width data type

pv = ConstantArray[$MaxMachineNumber, 1];

Developer`PackedArrayQ[pv]

lv = 2 * pv

Developer`PackedArrayQ[lv]

NumericArray

There are currently 12 NumericArray types compared to 3 for PackedArray

Integer8, Integer16, Integer32, Integer64

UnsignedInteger8, UnsignedInteger16, UnsignedInteger32, UnsignedInteger64

Real32, Real64

ComplexReal32, ComplexReal64

Memory/speed benefits.

Useful in compilation, audio/video processing, machine learning etc

Construction and Coercion

Use NumericArray to construct/convert to a machine array of a specific type.

Various of conversion methods

Default method attempts to coerce source type to target type.

nai8 = NumericArray[{1,2,3}, "Integer8"]

NumericArrayQ[nai8]

Conversion from one NumericArray type to another  is possible.

nar32 =  NumericArray[nai8, "Real32"]

Construction and Coercion

Normal converts a NumericArray to a PackedArray when possible, or a list otherwise.

par64 = Normal[nar32]

Developer`PackedArrayQ[par64]

PackedArray  type conversion is a special case of NumericArray type conversion and now works in parallel.

pa = Range[1*^8];
SerialParallelTiming[N[pa];]

Distinguishing between NumericArrays

Can create NumericArrays with identical values but different types

na1 = NumericArray[{1,2,3}, "Integer8"]

na2 = NumericArray[{1,2,3}, "UnsignedInteger32"]

SameQ distinguishes types:

SameQ[na1, na2]

Equal compares values

Equal[na1, na2]

Unification

The PackedArray data structure (MTensor) is now identical to the NumericArray data structure (MNumericArray)

Advantage is a single (C++ template based) code base that is more compact and can be better optimized.

Completely replaces old MTensor code with extended type version for MNumericArray.

Reimplemented using modern paradigms (e.g. C++ Standard Template Library and parallel directives).

Functionality is implemented in WolframRTL is now linked to the kernel as a dynamic library.

TreePlot[{"WolframRTL"->"Kernel","WolframRTL"->"New Compiler","WolframRTL"->"Old Compiler","WolframRTL"->"Compiler\nStandalone Code","WolframRTL"->"LibraryLink",
"WolframRTL"->"WolframNTL?","WolframRTL"->"FrontEnd?"},Center,"WolframRTL",VertexLabels->None,VertexShapeFunction->({LightBlue,EdgeForm[Black],Disk[#,.25],Black,Text[#2,#1]}&),PlotLabel->"Uses of WolframRTL"]

Thread safe memory allocator 4 times faster than malloc for small memory blocks

System Integration

V11.3 allows machine underflow and subnormal machine floating point numbers

Now allow smallest negative signed integer to be represented as a machine integer

The largest positive signed machine integer on a 64 bit machine

Developer`MachineIntegerQ[(2^63)-1]

The smallest negative signed machine integer on a 64 bit machine

Developer`MachineIntegerQ[-(2^63)]

The entire range of machine integers can now be represented in PackedArrays and NumericArrays

NumericArray[{-(2^63), (2^63)-1}, "Integer64"]

Negation of the smallest negative signed machine integer is converted to a big integer, otherwise x < 0 ⇏ -x > 0.

x = -(2^63);
mx = - x;
Map[Developer`MachineIntegerQ, {x, mx}]

Ordering and Sort

Mathematica has always used merge sort for ordering.

Merge sort is a stable sort where the relative order of ties is maintained.

Merge sort has a bounded sub-quadratic worst case complexity.

Merge sort has linear complexity for sorted arrays.

Merge sort uses fewer comparisons than quick sort, better when comparison is expensive.

Since a PackedArray is an implicit list representation, Ordering must use a stable sort to match the List case.

Ordering[Developer`ToPackedArray[{1,2,1,3,1}]]

Ordering[{1,2,1,3,1}]

Sort

Sorting a contiguous machine array in PackedArray or NumericArray does not have the same stability restriction on the data as a generic List.

The ordering of ties is unimportant since the data is identical.

Sort[Developer`ToPackedArray[{1,2,1,3,1}]]

Sort[{1,2,1,3,1}]

Alternative sorting methods (quick sort) are more efficient for contiguous arrays e.g. better data locality.

Parallel implementations for sorting machine arrays is becoming more widely available: both quick-sort and stable-sort

Parallel Sort Example

Sort an array in reverse order.

data = Range[1.*^8, 1., -1.];

SerialParallelTiming[Sort[data];]

SerialParallelTiming[Ordering[data];]

For comparison V12.0 takes: 5.15 seconds for Sort and 4.55 seconds for Ordering.

V12.0 implementation of Sort and Ordering:

Around 2500 lines of C code to support the 3 PackedArray types

New implementation:

About 500 lines of C++ code, essentially an interface to the C++ STL std::sort and std::stable_sort.

Supports all 12 NumericArray types

Works in parallel

Vector Parallelization

Want to avoid a naive threshold above which parallelization is used.

Could have many small blocks with large communication overhead just above the threshold

Instead

Vectors of length n are divided into p ≥ 2 blocks of size k

(({{(Subscript[V, < 1 >])}, {⋮}, {(Subscript[V, < p >])}}))    with (p - 1) k < n ≤ p k

Block-size chosen according to operation cost and data locality (e.g. fit into L1 cache)

Gradually use more threads as length increaes to keep communication overhead constant.

Parallel Functionality

Parallel functionality for NumericArray and PackedArray.

When possible use C++ STL functions for base case, with higher level parallelism.

Can remove parallelism layer when compilers fully support all aspects of the C++17 standard.

Constructors and copy/convert functions

ConstantArray, Range

Multiple Part assignment, Take, Drop

Flatten, ArrayReshape

NumericArray, ToPackedArray

Arithmetic

78 unary operations

21 binary operations

5 ternary operations e.g. fused-multiply-add

Map-reduction: Total, Max, MaxAbs, Min, MinMax, vector-Dot

Other parallel operations

Reverse, RotateLeft, RotateRight

Differences, Ratios

Sort, NumericalOrder, Order, Ordering, OrderedQ

Early termination

Machine array comparison is also now implemented in parallel

SameQ, Equal, Less etc and floating point exception checking

Parallel thread termination (cancelling) can be expensive

Blocks allow ‘early-termination’ for predicates and arithmetic exceptions

Block computation has a pretest for early termination

n = 5*^8;
pa1 = ConstantArray[1., n];
pa2 = ConstantArray[1., n];

Worst case for serial comparison: all elements the same

SerialParallelTiming[SameQ[pa1, pa2]]

Worst case for parallel comparison: first element different

pa2[[1]] = 2.;
SerialParallelTiming[SameQ[pa1, pa2]]

Parallel Transpose

Matrix transposition is not part of LAPACK.

Now have new parallel implementation of matrix transposition (together with higher rank arrays)

Matrices are explicitly divided into blocks chosen for data locality

Out-of-place transpose: straightforward

In-place-transpose: many different cases to consider

n × n : straightforward

n × m : with n | m or m | n efficient methods (n = 2 m used in complex FFT)

n × m : with n ∤ m and m ∤ n

Matrix Transposition (out-of-place)

pv = Range[1., 1.*^8,1.];

pm = ArrayReshape[pv, {1*^4, 1*^4}];

SerialParallelTiming[Transpose[pm];]

The serial case is around 20% faster than V12.0 because of better data locality.

Partition into p × q blocks of size k × k (typically k =  64)

(({{(Subscript[A, < 1 1 >]), ⋯, (Subscript[A, < 1 q >])}, {⋮, ⋱, ⋮}, {(Subscript[A, < p 1 >]), ⋯, (Subscript[A, <  p q >])}}))    with (p - 1) k < n ≤ p k and  (q - 1) k < m ≤ q k

Matrix Transposition (in-place)

Cycle-based algorithm (used in MKL, our implementation is slightly faster)

Compute permutation cycles to successively swap pairs of elements

Permutation cycles depend on matrix dimensions

0, 7, 14 are cycles of length 1 so their positions are unchanged in the transpose

({{A,  , (A^T)}, {(({{0, 5, 10}, {1, 6, 11}, {2, 7, 12}, {3, 8, 13}, {4, 9, 14}})),  ,  (({{0, 1, 2, 3, 4}, {5, 6, 7, 8, 9}, {10, 11, 12, 13, 14}}))}})

Matrix Transposition (in-place)

(1, 3, 9, 13, 11, 5) is a cycle

({{(({{0, 5, 10}, {1, 6, 11}, {2, 7, 12}, {3, 8, 13}, {4, 9, 14}})), (({{0, 5, 10}, {3, 6, 11}, {2, 7, 12}, {1, 8, 13}, {4, 9, 14}})), (({{0, 5, 10}, {3, 6, 11}, {2, 7, 12}, {9, 8, 13}, {4, 1, 14}})), (({{0, 5, 10}, {3, 6, 11}, {2, 7, 12}, {9, 8, 1}, {4, 13, 14}})), (({{0, 5, 10}, {3, 6, 1}, {2, 7, 12}, {9, 8, 11}, {4, 13, 14}})), (({{0, 1, 10}, {3, 6, 5}, {2, 7, 12}, {9, 8, 11}, {4, 13, 14}}))}})

(2, 6, 4, 12, 8, 10) is a cycle

({{(({{0, 1, 10}, {3, 6, 5}, {2, 7, 12}, {9, 8, 11}, {4, 13, 14}})), (({{0, 1, 10}, {3, 2, 5}, {6, 7, 12}, {9, 8, 11}, {4, 13, 14}})), (({{0, 1, 10}, {3, 4, 5}, {6, 7, 12}, {9, 8, 11}, {2, 13, 14}})), (({{0, 1, 10}, {3, 4, 5}, {6, 7, 2}, {9, 8, 11}, {12, 13, 14}})), (({{0, 1, 10}, {3, 4, 5}, {6, 7, 8}, {9, 2, 11}, {12, 13, 14}})), (({{0, 1, 2}, {3, 4, 5}, {6, 7, 8}, {9, 10, 11}, {12, 13, 14}}))}})

Swap the dimensions:

({{(({{0, 1, 2}, {3, 4, 5}, {6, 7, 8}, {9, 10, 11}, {12, 13, 14}})), ↦,  (({{0, 1, 2, 3, 4}, {5, 6, 7, 8, 9}, {10, 11, 12, 13, 14}}))}})

Matrix Transposition (in-place)

Cycle-based algorithm:

Elegant and simple to implement

Does not use cache efficiently because of irregular memory access patterns

Inherently sequential

Matrix Transposition (in-place)

New algorithm from NVIDIA researchers

A Decomposition for In-place Array Transposition, Catanzaro Keller & Garland, 2014

Consider data as 2-D array instead of 1-D array

Use row and column shuffles to reorder

In this example gcd(8, 4) = 2 so column rotation is required in the first phase with blocks of size 2

Matrix Transposition (in-place)

Implementation

Reduce expensive modular arithmetic (Arithmetic Strength Reduction).

Factor of 4 speedup for serial case.

Serial case becomes more efficient than cycle-based method.

Reduce auxiliary workspace and data movement

Arithmetic

Experimental kernel branch for comparison of single and double precision arithmetic

pv = Range[1., 5.*^8, 1.];

ByteCount[pv]

SerialParallelTiming[Sqrt[pv];]

nv = NumericArray[pv, "Real32"];

ByteCount[nv]

SerialParallelTiming[Sqrt[nv];]

MIxed types is a serious problem for binary operands:

na1 = NumericArray[{3,4,5}, "UnsignedInteger8"];
na2 = NumericArray[{1,2,3}, "Integer8"];
na1 + na2

Kernel vs Compiler

Proper support of NumericArray would require introduction of typed numbers

Example

Example of type information loss in kernel:

Part[NumericArray[{1,2,3}, "Integer8"], 1]

Example

Example of kernel semantics:

na1 = NumericArray[{3,4,5}, "UnsignedInteger8"];
na2 = NumericArray[{1,2,3}, "UnsignedInteger8"];

Expression na1 - na2 as subtraction is perfectly valid.

Kernel converts to Plus[na1, Minus[na2]] which involves negation of an unsigned type

Kernel vs Compiler

Example

In the following should convert integral constants to Typed[x, "UnsignedInteger8"]

na = NumericArray[{1,2,3}, "UnsignedInteger8"];
1 + 2 (na^2)

Example

Can do constant type markup and dead code elimination in compiler

fun = Function[Typed[x,"UnsignedInteger8"],
 If [x < Typed[0,TypeOf[x]],
x+ Typed[1,TypeOf[x]],
 x+ Typed[2,TypeOf[x]]
]
]

cir=CompileToIR[fun];
cir["toString"]

Compiler

PackedArray interface for new compiler:

pafun=Function[Typed[arg,TypeSpecifier["PackedArray"]["Real64",1]],Total[arg]];
cpafun = FunctionCompile[pafun]

cpafun[Range[1., 10., 1.]]

Extension of compiler interface for NumericArray:

nafun=Function[Typed[arg,TypeSpecifier["NumericArray"]["Real32",1]],Total[arg]];

Conclusions

NumericArray is an extension of PackedArray including more native machine number types

PackedArray data structure is now equivalent to the NumericArray data structure.

Arithmetic and structural operations are performed in WolframRTL

OpenMP primitives used for parallelization

PackedArray operations are becoming faster irrespective of the compiler.

WolframRTL ongoing rewrite using modern C++

WolframRTL code base over time:

({{, }})

Immediate Future

Integration of NumericArray operations into the new compiler

More structural operations

Typed random number generation

Longer term

Dense linear algebra

Sparse linear algebra

Integration into kernel? Need to interact with:

Numbers

PackedArrays

Lists

SparseArrays

StructuredArrays