Hyper-systolic algorithms for N-body computations and parallel level-3 BLAS libraries

Coded Distributed Matrix Multiplication (CDMM) is a distributed matrix multiplication (DMM) for large-scale matrices through a coding scheme such that any <italic xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink">R</i> worker nodes among all <italic xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink">N</i> worker nodes can recover the final product, where <italic xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink">N</i> corresponds to the length of the code and <italic xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink">R</i> ≤ <italic xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink">N</i> is called the recovery threshold. The state-of-art CDMM schemes, such as EP codes for Single DMM and GCSA codes for batch DMM, are defined over a Galois field GF(<italic xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink">q</i>) of size <italic xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink">q</i> ≥ <italic xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink">N</i>. These are inefficient for small Galois fields such as GF(2) and the integer residue ring Z<sub xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink"><i>p</i></sub><sup xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink">ℓ</sup> due to the lack of invertible elements for interpolation. DMM over Z<italic xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink">_p</i><sup xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink">ℓ</sup> (such as Z<sub xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink">2</sub><sup xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink">64</sup> ) is well-motivated in practice due to their direct compatibility with hardware. In this work, we construct efficient CDMM over the Galois ring GR(p<sup xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink">ℓ</sup>, <italic xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink">d</i>) which is an extension ring over Z<sub xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink"><i>p</i></sub><sup xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink">ℓ</sup> of degree <italic xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink">d</i>, particularly, GR(<italic xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink">p, d</i>) = GF(<italic xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink">p</i><sup xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink">d</sup>) is the Galois field and GR(<sub xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink"><i>p</i></sub><sup xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink">ℓ</sup>, 1) = Z<italic xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink">_p</i><sup xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink">ℓ</sup>. We first give a general CDMM framework for the batch of <italic xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink">n</i> matrix multiplications via the famous RMFE (Cascudo et al. Crypto’18). Compared with GCSA, our construction has a smaller recovery threshold by a factor of 1/<italic xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink">n</i>. Next, we optimize EP codes via batch preprocessing of the input matrices. We give two types of Single CDMM, which can achieve almost the same performance as EP codes over a Galois field with size <italic xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink">q</i> ≥ <italic xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink">N</i>. Finally, we present the experimental analysis of our CDMM on Galois rings.

We propose two coding schemes for distributed matrix multiplication in the presence of stragglers. These coding schemes are adaptations of Luby Transform (LT) codes and Raptor codes to distributed matrix multiplication and are termed <italic xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink">Factored LT (FLT) codes</i> and <italic xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink">Factored Raptor (FRT) codes</i> . We show that all nodes in the Tanner graph of a randomly sampled code have a tree-like neighborhood with high probability. This ensures that the density evolution analysis gives a reasonable estimate of the average recovery threshold of FLT codes. The recovery threshold of the proposed FLT codes is asymptotically optimal when the output degree distribution is Soliton. Empirically, we show that FRT codes have an excellent recovery threshold while the number of worker nodes is moderately large. In addition, using Azuma–Hoeffding inequality, we derive concentration results to show that the recovery threshold of a randomly chosen FLT code is close to the ensemble average. FLT and FRT codes have better recovery thresholds when compared to Product codes and they are expected to have better numerical stability when compared to Polynomial codes, while they can also be decoded with a low-complexity decoding algorithm. Finally, the proposed codes are better matched to the practically important case of sparse matrix-matrix multiplication as compared to many previous schemes.

This paper introduces a novel master-multi-SIMD on-chip multi-core architecture for embedded signal processing. The parallel architecture and its memory subsystem are described in this paper. We evaluate the large size matrix multiplication performance on this parallel architecture and compare it with a SIMD-extended data parallel architecture. We also examine how well the new architecture scales for different numbers of SIMD co-processors. The experimental results show that the ePUMA architecture’s memory subsystem can effectively hide the data access overhead. With its 8-way SIMD data path and multi-SIMD parallel execution, the ePUMA architecture improves the performance of matrix multiplication with a speedup of 45x from the conventional SIMD extension.

The distributed matrix multiplication problem with unknown number of stragglers is considered, where the goal is to allow a master to efficiently and flexibly obtain the product of two massive matrices by distributing the computation across <tex xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink">$N$</tex> servers. We assume there are at most <tex xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink">$N-R$</tex> stragglers but the exact number is not known a priori. Motivated by reducing the latency, a flexible solution is proposed to fully utilize the computation capability of available servers. The computing job for each server is separated into 2 layers, constructed based on Entangled Polynomial (EP) codes by Yu el al. The final results can be obtained when a larger number of servers complete the task from the first layer or a smaller number of servers complete the tasks from both 2 layers. The required finite field size of the proposed solution is less than <tex xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink">$2N$</tex> . Moreover, the optimal partitioning of the input matrices is discussed. Our constructions can also be generalized to batch matrix multiplication.

In this paper, a unified algebraic transformation approach is presented for designing parallel recursive and adaptive digital filters and singular value decomposition (SVD) algorithms. The approach is based on the explorations of some algebraic properties of the target algorithms' representations. Several typical modern digital signal processing examples are presented to illustrate the applications of the technique. They include the cascaded orthogonal recursive digital filter, the Givens rotation-based adaptive inverse QR algorithm for channel equalization, and the QR decomposition-based SVD algorithms. All three examples exhibit similar throughput constraints. There exist long feedback loops in the algorithms' signal flow graph representation, and the critical path is proportional to the size of the problem. Applying the proposed algebraic transformation techniques, parallel architectures are obtained for all three examples. For cascade orthogonal recursive filter, retiming transformation and orthogonal matrix decompositions (or pseudo-commutativity) are applied to obtain parallel filter architectures with critical path of five Givens rotations. For adaptive inverse QR algorithm, the commutativity and associativity of the matrix multiplications are applied to obtain parallel architectures with critical path of either four Givens rotations or three Givens rotations plus two multiply-add operations, whichever turns out to be larger. For SVD algorithms, retiming and associativity of the matrix multiplications are applied to derive parallel architectures with critical path of eight Givens rotations. The critical paths of all parallel architectures are independent of the problem size as compared with being proportional to the problem size in the original sequential algorithms. Parallelism is achieved at the expense of slight increase (or the same for the SVD case) in the algorithms' computational complexity.

In distributed matrix multiplication, a common scenario is to assign each worker a fraction of the multiplication task, by partitioning the input matrices into smaller submatrices. In particular, by dividing two input matrices into m-by-p and p-by-n subblocks, a single multiplication task can be viewed as computing linear combinations of pmn submatrix products, which can be assigned to pmn workers. Such block-partitioning based designs have been widely studied under the topics of secure, private, and batch computation, where the state of the arts all require computing at least "cubic" (pmn) number of submatrix multiplications. Entangled polynomial codes, first presented for straggler mitigation, provides a powerful method for breaking the cubic barrier. It achieves a subcubic recovery threshold, meaning that the final product can be recovered from any subset of multiplication results with a size order-wise smaller than pmn. In this work, we show that entangled polynomial codes can be further extended to also include these three important settings, and provide a unified framework that order-wise reduces the total computational costs upon the state of the arts by achieving subcubic recovery thresholds.

In this paper, we propose a parallel processing architecture, based on two TMS320C44 VME-bus DSP boards, for real time implementation of Cell Averaging and Clutter Map Constant False Alarm Rate (CFAR) detectors, and a data extractor based on a centroidal interpolation of the position of the detected target. The optimal processing speed has been achieved by fully exploiting the capacities of the dasiaC44 processor. The implemented system is well adapted for two dimension radars. The overall processing scheme interconnections and the real time implementation results are presented and discussed.

The goal of the paper is to present a novel approach for designing massively parallel architectures for image recognition systems. Massive parallel hardware produces high speed image recognition. Hardware architectures work much faster and are more effective for application specific processing than a software-based image recognition system. The concept of designing optimal image recognition systems with high speed, low complexity, good portability, and low cost has been present for a long time. The authors have considered the basic design factors, simple and regular design, concurrency and communication, and balancing computation with I/O. The technique presented provides a new method of designing an image recognition system by concentrating on modularity and massive parallelism. This is accomplished through the design of a processing element layer, local communicator layer, and network of new image layer. The processing element (PE) layer performs the basic recognition functions on individual pixels through logical operations. Both the local communicator and network of new image layers institute a mesh network topology, which has proven to be the fastest and most efficient. The local communicator layer has the most important function, it communicates with the other two layers, combines results from the PE layer, and houses the main control unit. All layers work simultaneously, or in other words, in a parallel fashion. The flow diagram of a CAD designer algorithm, based on this new design technique, is presented.

Free text retrieval is an important problem which can significantly benefit from a parallel architecture. Signature methods have been proposed to answer text retrieval queries in parallel machines [Sta88, LF92], under the assumption that the main memory is sufficient to hold the entire signature file. We propose the use of a Parallel Bit-Sliced Signature File method on a SIMD machine architecture when the size of the signature file exceeds the available memory. We propose that we need not examine all the bit slices; instead we use a partial fetch slice swapping algorithm. This method achieves graceful performance degradation according to the database size. We provide formulae for the optimal number of signature slices to fetch and match with the query signature. Arithmetic examples show that our method can handle a 128GB database with a 2sec response time on a machine with the characteristics of the Connection Machine.

Highly parallel computer architectures and algorithms for speech analysis operations are surveyed. The classes of architectures considered include SIMD, skewed-SIMD, MIMD, and data flow machines, associative processors, and pipelined systolic and wavefront arrays. Parallel algorithms for digital filtering, FFTs, and linear predictive coding are summarized. System requirements are analyzed in terms of number and complexity of processors, memory requirements, and nature and complexity of inter-processor communications.

The kinematics, dynamics, Jacobian, and their corresponding inverses are six major computational tasks in the real-time control of robot manipulators. The parallel algorithms for these computations are examined and analyzed. They are characterized based on six well-defined features that have greatest effects on the execution of parallel algorithms. These features include type of parallelism, degree of parallelism (granularity), uniformity of operations, fundamental operations, data dependency, and communication requirement. It is found that the inverse dynamics, the forward dynamics, the forward kinematics and the forward Jacobian computations possess highly regular properties and they are all in homogeneous linear recursive form. The inverse Jacobian is essentially the problem of solving a system of linear equations. The closed-form solution of the inverse kinematics problem is obviously non-uniform and robot dependent. The iterative solution for the inverse kinematics problem seems uniform and the parallel portions of the algorithm involve the forward kinematics, the forward Jacobian, and the inverse Jacobian computations. Suitable algorithms for the six basic robotics computations are selected and parallelized to make use of their common features. The characterization of the six basic robotics algorithms is tabulated for discussion and the results can be used to design better parallel architectures or a common architecture for the computation of these robotics algorithms.

An optical packet switch (OPS) is called performance guaranteed if it can achieve 100% throughput with bounded packet delay. Presently, high speedup requirement and large packet delay are two main disadvantages in designing performance guaranteed OPS. Survivability is another important issue that must be considered for real OPS implementations. In this paper, we propose a two-layer parallel OPS architecture together with an efficient scheduling scheme to address all the above issues. The tradeoff between speedup and packet delay under this new parallel architecture is also formulated to provide more design flexibility. Compared to the single-layer OPS, our proposed solution can simultaneously reduce both speedup and packet delay. For example, a delay of 4δN slots can be achieved with a speedup of 2 in our solution (where N is the switch size and δ is the switch reconfiguration overhead), whereas the single-layer OPS needs a speedup of 6 for a delay of 7δN slots. We show that this significant improvement benefits from a careful overall design rather than simply adding an extra switching layer. © 2005 IEEE.

A nonlinear, nonrecursive, two-dimensional filter has been developed for edge-preserving noise smoothing of images. It can also be beneficial as a pre-processor for data compression purposes. The filter is applied over a 3x3 window of pixels and uses an adaptive threshold. As with all image filters, the computational requirements lead to the utillization of special purpose VLSI hardware for real time implementation. The algorithm is being implemented in VLSI using a systolic array structure and parallel architecture, in which each processing element is assigned to one pixel. The structure of the processing element is very simple: except for some control gates and registers, it contains only a one-bit full adder as the main functional element. The most important feature of the proposed architecture is that the nonlinear operation has been implemented using no more than the simple tools and the silicon area required for a linear convolutional algorithm. Other properties of the above design are: simple interconnection scheme between processing elements and common control circuitry for all the processing elements on a single chip. Due to the topological simplicity of the processing element, arrays of size 16x16 cells per chip are believed to be feasible, using standard CMOS technology.

Special issue on <i>specialized</i> computer architecture simulators that see the present and may hold the future

This paper compares two general library routines for performing parallel distributed matrix multiplication. The PUMMA algorithm utilities block scattered data layout, whereas BiMMeR utilizes virtual 2-D torus wrap. The algorithmic differences resulting from these different layouts are discussed us well as the general issues associated with different data layouts for library routines. Results on the Intel Delta for the two matrix multiplication algorithms are presented. >