Parallel algorithms development for programmable logic devices

TL;DR

This paper presents a formal methodology for designing massively parallel algorithms on programmable logic devices, using functional specifications, refinement techniques, and CSP-inspired operators, demonstrated through matrix multiplication implementations on FPGA hardware.

Contribution

It extends formal hardware synthesis methodologies with a step-wise refinement approach using functional programming and CSP-inspired operators for efficient parallel algorithm design on PLDs.

Findings

Multiple hardware implementations with different performance trade-offs.

Successful compilation and testing on FPGA hardware.

Performance analysis of various algorithm refinements.

Abstract

Programmable Logic Devices (PLDs) continue to grow in size and currently contain several millions of gates. At the same time, research effort is going into higher-level hardware synthesis methodologies for reconfigurable computing that can exploit PLD technology. In this paper, we explore the effectiveness and extend one such formal methodology in the design of massively parallel algorithms. We take a step-wise refinement approach to the development of correct reconfigurable hardware circuits from formal specifications. A functional programming notation is used for specifying algorithms and for reasoning about them. The specifications are realised through the use of a combination of function decomposition strategies, data refinement techniques, and off-the-shelf refinements based upon higher-order functions. The off-the-shelf refinements are inspired by the operators of Communicating Sequential Processes (CSP) and map easily to programs in Handel-C (a hardware description language). The Handel-C descriptions are directly compiled into reconfigurable hardware. The practical realisation of this methodology is evidenced by a case studying the matrix multiplication algorithm as it is relatively simple and well known. In this paper, we obtain several hardware implementations with different performance characteristics by applying different refinements to the algorithm. The developed designs are compiled and tested under Celoxica's RC-1000 reconfigurable computer with its 2 million gates Virtex-E FPGA. Performance analysis and evaluation of these implementations are included.

Equations

$$\offinterlineskip\halign to=433.62pt{\strut$\@zlign#$\hfil\cr PRD \penalty 10000\ (Item \penalty 10000\ a) = out.element.channel \penalty 10000\ ! \penalty 10000\ a \rightarrow SKIP \crcr}$$

$$\offinterlineskip\halign to=433.62pt{\strut$\@zlign#$\hfil\cr\begin{array}[]{@{}lllll}PRD \penalty 10000\ (\langle{s}\rangle) = ((;)_{i=1}^{i=length(s)} (PRD \penalty 10000\ s_{i}) [out.elements.channel / out]);\\ out.eotChannel \penalty 10000\ ! \penalty 10000\ eot \rightarrow SKIP\end{array}\crcr}$$

$$\offinterlineskip\halign to=433.62pt{\strut$\@zlign#$\hfil\cr\begin{array}[]{@{}lllll}PRD \penalty 10000\ (\lfloor v \rfloor_{n})&=&\mathrel{\mathord{\mid}\mathord{\mid}\mathord{\mid}}_{i=1}^{i=n} (PRD \penalty 10000\ v_{i}) [ out.elements_{i}.channel / out ]\end{array}\crcr}$$

$$\offinterlineskip\halign to=433.62pt{\strut$\@zlign#$\hfil\cr P \rhd Q = (P [mid / out] \penalty 10000\ || \penalty 10000\ Q [mid / in]) \verb"\" \{ mid \} \crcr}$$

$$\offinterlineskip\halign to=433.62pt{\strut$\@zlign#$\hfil\cr f :: A \rightarrow B \crcr}$$

$$\offinterlineskip\halign to=433.62pt{\strut$\@zlign#$\hfil\cr f \sqsubseteq F \crcr}$$

$$\offinterlineskip\halign to=433.62pt{\strut$\@zlign#$\hfil\cr SMAP(F) =\crcr\mathop{\mathstrut{\mathchar 278\relax}}\nolimits X \penalty 10000\ \bullet in.eotChannel \penalty 10000\ ? \penalty 10000\ eot \rightarrow out.eotChannel \penalty 10000\ ! \penalty 10000\ eot \rightarrow SKIP \crcr\mathrel{\Box}\crcr F [in.elements.channel / in, out.elements.channel / out] ; X \crcr}$$

$$\offinterlineskip\halign to=433.62pt{\strut$\@zlign#$\hfil\cr\begin{array}[]{@{}lllll}VMAP_{n}(F)&=&\mathrel{\mathord{\mid}\mathord{\mid}\mathord{\mid}}_{i=1}^{i=n} F [ in_{i} / in, out_{i} / out]\end{array}\crcr}$$

$$\offinterlineskip\halign to=433.62pt{\strut$\@zlign#$\hfil\cr zipWith :: (A\rightarrow B\rightarrow C) \rightarrow[A]\rightarrow[B]\rightarrow[C] \crcr}$$

$$\offinterlineskip\halign to=433.62pt{\strut$\@zlign#$\hfil\cr zipWith \penalty 10000\ (\oplus) \penalty 10000\ [x_{1}, x_{2}, ... x_{n}] [y_{1}, y_{2}, ... y_{n}] = [x_{1} \oplus y_{1}, x_{2} \oplus y_{2},..., x_{n} \oplus y_{n}] \crcr}$$

$$\offinterlineskip\halign to=433.62pt{\strut$\@zlign#$\hfil\cr\begin{array}[]{lcl}SZIPWITH(F) =\\ \mathop{\mathstrut{\mathchar 278\relax}}\nolimits X \penalty 10000\ \bullet in.eotChannel \penalty 10000\ ? \penalty 10000\ eot \rightarrow out.eotChannel \penalty 10000\ ! \penalty 10000\ eot \rightarrow SKIP\\ \mathrel{\Box}\\ F [in_{1}.elements.channel / in_{1}, in_{2}.elements.channel / in_{2},\\ out.elements.channel / out] ; X\end{array} \crcr}$$

$$\offinterlineskip\halign to=433.62pt{\strut$\@zlign#$\hfil\cr vzipWith \penalty 10000\ (\oplus):: \lfloor A \rfloor_{n} \rightarrow,\lfloor B \rfloor_{n} \rightarrow\lfloor C \rfloor_{n} \crcr}$$

$$\offinterlineskip\halign to=433.62pt{\strut$\@zlign#$\hfil\cr VZIPWITH \penalty 10000\ (\oplus) = \mathrel{\mathord{\mid}\mathord{\mid}\mathord{\mid}}_{i=1}^{i=n}F[out_{i}/out,c_{i}/in_{1},d_{i}/in_{2}] \crcr}$$

Full text

1 Introduction

The rapid progress and advancement in electronic chips technology provides a variety of new implementation options for system engineers. The choice varies between the flexible programs running on a general purpose processor (GPP) and the fixed hardware implementation using an application specific integrated circuit (ASIC). Many other implementation options present, for instance, a system with a RISC processor and a DSP core. Other options include graphics processors and microcontrollers. Specialist processors certainly improve performance over general-purpose ones, but this comes as a quid pro quo for flexibility. Combining the flexibility of GPPs and the high performance of ASICs leads to the introduction of reconfigurable computing (RC) as a new implementation option with a balance between versatility and speed.

Generally, reconfigurable computing is computer processing with highly flexible computing fabrics. The principal difference when compared to using ordinary microprocessors is the ability to make substantial changes to the data path itself in addition to the control flow. In the last decade, there was a renaissance in the area of reconfigurable computing research with many proposed reconfigurable architectures developed both in industry and academia such as, Matrix, Garp, RAW, DPGA, RaPiD, PRISM, Pleiades, and Morphosys [1]. Such designs were feasible due to the relentless progress of silicon technology that allowed complex designs to be implemented on a single chip.

Field Programmable Gate Arrays (FPGAs), nowadays are important components of RC-systems, have shown a dramatic increase in their density over the last few years. For example, companies like Xilinx [2] and Altera [3] have enabled the production of FPGAs with several millions of gates, such as in Virtex-II Pro and Stratix-II FPGAs. The versatility of FPGAs, opened up completely new avenues in high-performance computing. These reconfigurable digital electronic hardware circuits can be combined with high-level software and design methodologies to form a powerful paradigm for computing.

The traditional implementation of a function on an FPGA is done using logic synthesis based on VHDL, Verilog or a similar HDL (hardware description langauge). These discrete event simulation languages are rather different from languages, such as C, C++ or JAVA. Many FPGA implementation tools are primarily HDL-based and not well integrated with high-level software tools. Furthermore, these HDL-based IP (intellectual property) cores are expensive and they have complex licensing schemes [4]. These obstacles had caused some blockage to the infiltration of FPGAs as the main platform solution for hardware engineers. An interesting step towards more success in hardware compilation is to grant a higher-level of abstraction from the point of view of programmer. Designer productivity can be improved and time-to-market can be reduces by making hardware design more like programming in a high-level langauge. Recently, vendors have initiated the use of high-level languages dependent tools like Handel-C [5, 6, 7, 8], Forge [9], Nimble [10, 11], SystemC [12] and Viva [13] (an object-oriented graphical development environment for programming FPGAs).

With the availability of powerful high-level tools accompanying the emergence of multi-million FPGA chips, more emphasis should be placed on affording an even higher level of abstraction in programming reconfigurable hardware. Building on these research motivations, in the work in hand, we extend and examine a methodology whose main objective is to allow for a higher-level correct synthesis of massively parallel algorithms and to map (compile) them onto reconfigurable hardware. Our main concern is with behavioural refinement, in particular the derivation of parallel algorithms. The presented methodology systematically transforms functional specifications of algorithms into parallel hardware implementations. It builds on the work of Abdallah and Hawkins [14, 15, 16, 17] extending their treatment of data and process refinement.

This paper is divided so that the following section introduces the adopted development methodology. Section 3 presents the theoretical background. In Section 4, we put some emphasis on the approach to develop different implementations of the matrix multiplication algorithm. The following section details the development steps. Section 7 demonstrates selected implementations. In Section 8, we analyze and evaluate the performance of the suggested implementations. Finally, Section 10 concludes the paper.

2 The Development Method

Although compilers can expose parallelism through data flow analysis [18], imperative languages are perhaps not ideal as a starting point. This is because imperative programs already incorporate design decisions (concerning control flows and data structures), preconditions (that can be assumed), post-conditions (that must be achieved), and invariants (that must be maintained). The direct manipulation of state makes it both difficult to prove that any two pieces of code are equivalent, and to perform substitutions, modify and rewrite the algorithm. Functional languages [19], such as Haskell [20], however, do not manipulate state directly, and as such gain the property of referential transparency. Any sub-expression of an algorithm can be substituted for any other that is provably equivalent. This is aided by an effective set of laws given to us by such reasoning frameworks as Bird-Merteen Formalism (BMF) [21], along with a wealth of other work in the functional programming and parallel processing fields [22, 23, 24, 25, 26, 27].

Although, many hardware development methods still use the powerful data flow analysis, such as Viva [13], the attractions for using the functional paradigm has incited many researchers. This triggered many investigations in this area, such as Lava [28], Hawk [29, 30], Hydra [31], HML [32], MHDL [33], DDD system [34], SAFL [35], MuFP [36], Ruby [37], and Form [38].

The suggested development model adopts the transformational programming approach for deriving massively parallel algorithms from functional specifications (See Figure 1). The functional notation is used for specifying algorithms and for reasoning about them. This is usually done by carefully combining a small number of higher-order functions that serve as the basic building blocks for writing high-level programs. The systematic methods for massive parallelisation of algorithms work by carefully composing an ”off-the-shelf” massively parallel implementation of each of the building blocks involved in the algorithm. The underlying parallelisation techniques are based on both pipelining and data parallelism.

Higher-order functions, such as map, filter, foldl, and foldr, provide a high degree of abstraction in functional programs [20]. Not only they do allow clear and succinct specifications for a large class of algorithms, but they also are ideal starting points for generating efficient implementations by a process of mathematical calculation using BMF. Over the past decade, there have been attempts to apply BMF for generating data parallel programs from abstract specifications using the skeleton approach [24, 15]. The main attraction of this approach is the potential for increasing reusability of parallel programs without sacrificing too much performance. The essence of this approach is to design a generic solution once, and to use instances of the design many times for various applications. Accordingly, this approach allows portability by implementing the design on different parallel architectures.

In order to develop generic solutions for general parallel architectures, it is necessary to formulate the design within a concurrency framework such as CSP [15, 8]. Often parallel functional programs show peculiar behaviours which are only understandable in the terms of concurrency rather than relying on hidden implementation details. The formalisation in CSP (of the parallel behaviour) leads to better understanding and allows for analysis of performance issues. The establishment of refinement concepts between functional and concurrent behaviours may allow systematic generation of parallel implementations for various architectures. This gives the ability to exploit well-established functional programming (FP) paradigms and transformation techniques in order to develop efficient CSP processes. These systematic refinement rules refine the specification to what we call the CSP implementation stage, where parallelism will be described using Hoare’s CSP [39]. Again, this allows issues of immense practical importance, such as, the careful reasoning about data distribution, network topology, and locality of communications.

The previous stages of development require a back-end stage for realising the developed designs. We note at this point that the Handel-C language relies on the parallel constructs in CSP to model concurrent hardware resources. Mostly, algorithms described with CSP could be implemented with Handel-C. Accordingly, this langauge is suggested as the final reconfigurable hardware realisation stage in the proposed methodology. It is noted that, for the desired hardware realisation, Handel-C enables the integration with VHDL and EDIF (Electronic Design Interchange Format) and thus various synthesis and place-and-route tools.

3 Background

Abdallah and Hawkins defined in [16] some constructs used in the development model. Their investigation looked in some depth at data refinement; which is the means of expressing structures in the specification as communication behaviour in the implementation.

3.1 Data Refinement

In the following we present some datatypes used for refinement, these are stream, vector, and combined forms.

The stream is a purely sequential method of communicating a group of values. It comprises a sequence of messages on a channel, with each message representing a value. Values are communicated one after the other. Assuming the stream is finite, after the last value has been communicated, the end of transmission (EOT) on a different channel will be signaled. Given some type A, a stream containing values of type A is denoted as $\langle A\rangle$.

Each item to be communicated by the vector will be dealt with independently in parallel. A vector refinement of a simple list of items will communicate the entire structure in a single. Given some type A, a vector of length n, containing values of type A, is denoted as $\lfloor A\rfloor_{n}$.

Whenever dealing with multi-dimensional data structures, for example, lists of lists, implementation options arise from differing compositions of our primitive data refinements - streams and vectors. Examples of the combined forms are the Stream of Streams, Streams of Vectors, Vectors of streams, and Vectors of Vectors. These forms are denoted by: $\langle S_{1},S_{2},...,S_{n}\rangle$ , $\langle V_{1},V_{2},...,V_{n}\rangle$, $\lfloor S_{1},S_{2},...,S_{n}\rfloor$, and $\lfloor V_{1},V_{2},...,V_{n}\rfloor$.

3.2 Process Refinement

The refinement of the formally specified functions to processes is the key step towards understanding possible parallel behaviour of an implementation. In this section, the interest is in presenting refinements of a subset of functions - some of which are higher-order. A bigger refined set of these functions is discussed in [15].

Generally, These highly reusable building blocks can be refined to CSP in different ways. This depends on the setting in which these functions are used (i.e. with streams, vectors etc.), and leads to implementations with different degrees of parallelism. Note that we don’t use CSP in a totally formal way, but we use it in a way that facilitates the Handel-C coding stage later. Recall for the following subsections that values are communicated through as an elements channel, while a single bit is communicated through another eotChannel channel to signal the end of transmission (EOT).

3.2.1 Produce

The produce process (PRD) is fundamental to process refinement. It is used to produce values on the channels of a certain communication construct (Item, Stream, Vector, and so on). These values are to be received and manipulated by another processes.

Items

For simple, single item types (int, char, bool, etc.), the produce process is very simple. This is depicted in Figure 2. Here the output is just a single channel.

The definition in CSP notation is very straightforward:

[TABLE]

Streams

The produce process for streams is depicted in Figure 3. As already noted, the output in this case is a pair of two other channels. One channel will produce the values of the stream, and the other will be a simple channel used to signal EOT.

In a more general case, the structure of the values which the stream is carrying is not necessarily known. These may be simple items, but may also be streams or vectors. Generally, producing a stream could be described as:

[TABLE]

Vectors

For vectors of size $n$, $n$ instances of the produce process are composed in parallel, one for each item in the vector. The output here is an array of channels. This is depicted in Figure 4. A general definition is given below:

[TABLE]

A process STORE stores a communication construct in a variable. We use this process to store items, vectors, streams, or combinations. A subscript letter is used with the processes PRD and STORE to indicate the type of communication. We sometimes omit this subscript if the communication structure is clear from context.

3.2.2 Feeding Processes

The feed operator in CSP models function application. The feed operator is written . The feed operator takes two processes, composes them together in parallel, and renames both the output of the first and the input of the second to a new name, which is then hidden. Given the lifted concepts of CSP channel renaming and hiding, the definition can remain the same regardless of the type of the communicating construct (Item, Stream, Vector or any combination).

[TABLE]

3.2.3 Formal Process Refinement

Given the definition of a feed operator that operates on processes, a formal definition of process refinement could be delivered. Consider a function $f$, which takes in values of type $A$ and returns values of type $B$. Assume that the data refinement step has already been performed, such that $A$ and $B$ are both types of some transmission value:

[TABLE]

Then, consider a potential refinement for $f$, a process $F$. The operator $\sqsubseteq$ denotes a process refinement, where the left hand side is a function, and the right hand side is a process. To state that $f$ is refined to $F$, or in other words, the process $F$ is a valid refinement of the function $f$, the following may be used:

[TABLE]

The rules of refinement were proven once [15], and in this paper we use them systematically to refine the functional specification into a network of communicating processes.

3.2.4 MAP the Process Refinement of the Higher-order Function map

Now the attention is turned to the refinement of the widely used higher-order function map [16] . Employing this function in stream and vector settings is presented. The refinement for combined structures is to be made in a similar way.

Streams

A process implementing the functionality of $map\,f$ in stream terms should input a stream of values, and output a stream of values with the function $f$ applied (See Figure 5).

In general, the handling of the EOT channels will be the same. However, the handling of the value will vary depending on the type of the elements of the input and output stream.

[TABLE]

Vectors

In functional terms, the functionality of $map\,f$ in a list setting is modelled by $vmap\,f$ in the vector setting. Consider $F$ as a valid refinement of the function $f$. The implementation of $VMAP$ can then proceed by composing $n$ instances of $F$ in parallel, and directing an item from the input vector to each instance for processing (See Figure 6). In CSP we have:

[TABLE]

3.2.5 ZIPWITH the Process Refinement of the Higher-order Function zipWith

Recall another higher-order function, namely zipWith. This function is used to zip two lists (taking one element from each list) with a certain operation. Formally:

[TABLE]

[TABLE]

3.2.6 Streams

The process implementation of ($zipWith\,f)$ in stream terms should input two streams of values, and output a stream of values with the function $f$ applied (See Figure 7).

Again, the handling of the EOT channel will be the same. Nevertheless, the handling of the value will vary depending on the type of the input and output streams elements.

[TABLE]

3.2.7 Vectors

To implement the data parallel version of this higher-order function, we refine it to a process VZIPWITH that takes two vectors as input and zips the two lists with a process F; F is a refined process from the function ($\oplus$). This is depicted as in Figure 8.

[TABLE]

[TABLE]

3.3 Decomposition of Higher-Order Functions for Parallel Programs Derivation

Intrinsic richness and usefulness of higher-order functions could be made clear by recalling some of previous work presented in [40]. This work concentrated on providing systematic decomposition methods for exploiting pipelined parallelism in instances of the higher-order function foldr. In this section, the decomposition for the higher-order function map is shown. This decomposition rule is used in the forthcoming applications.

The following decomposition rule is recalled along with its corresponding CSP implementation. This rule decomposes a specification of the form (map (h m)), where h is a function and m is a given list of values. The CSP network SPEC which refines this specification is shown in Figure 9.

$spec::A\rightarrow B\mathchar 24635\relax\;h::[T]\rightarrow A\rightarrow B\mathchar 24635\relax\;m::[T]\mathchar 24635\relax\;e::B$ $f::T\rightarrow A\rightarrow B\rightarrow B$ $spec=map$  $(h$   $m)$ $h$  $[]$  $x$  $=$  $e$ $h$  $(a$  $:$  $s)$   $x$  $=$  $f$  $a$  $x$  $(h$  $s$  $x)$

This could be decomposed to:

$spec=(\circ)/[final^{{}_{\ast}}]$  ${++}$  $(map\circ f^{\prime}\ast m)$  ${++}$  $[initial^{{}_{\ast}}]$ $f^{\prime}$  $a$  $\langle x,y\rangle$  $=$  $\langle x,$  $f$  $a$  $x$  $y\rangle$ $initial$  $x$ $=\langle x,$  $y\rangle$ $final$   $\langle x,$  $y\rangle=y$

The pipelined network of CSP processes SPEC, which refines the functional specification spec is synthesised as follows:

$SPEC=(\gg)/[MAP(initial)]\mathrel{{+}\mkern-8.0mu{+}}((MAP\circ f^{\prime})\ast(reverse$  $m))\mathrel{{+}\mkern-8.0mu{+}}[MAP(final)]$ $MAP(initial)=\mathop{\mathstrut{\mathchar 278\relax}}\nolimits X$  $\bullet$    $left?eot\rightarrow right!eot\rightarrow SKIP$

$|$

$left?x\rightarrow right!\langle x,e\rangle\rightarrow X$

$MAP(f^{\prime}$  $as)=\mathop{\mathstrut{\mathchar 278\relax}}\nolimits X$  $\bullet$   $left?eot\rightarrow right!eot\rightarrow SKIP$

$|$

$left?\langle x,y\rangle\rightarrow right!\langle x,f$  $a$  $x$  $y\rangle\rightarrow X$

$MAP(final)=\mathop{\mathstrut{\mathchar 278\relax}}\nolimits X$  $\bullet$    $left?eot\rightarrow right!eot\rightarrow SKIP$

$|$

$left?\langle x,y\rangle\rightarrow right!y\rightarrow X$

It is important not jump to the conclusion that every parallel algorithm can be deve1oped this way. There are two limitations to this approach. First, it only deals with systems which can be specified functionally. Second, it may not be possible to develop some parallel algorithms which use multi-directional communications using this method. This second point will be practically assessed in later section while designing a multilevel-pipelined parallel program.

3.4 Handel-C as a Stage in the Development Model

Based on datatype refinement and the skeleton afforded by process refinement, the desired reconfigurable circuits are built. Circuit realisation is done using Handel-C, as it is based on the theories of CSP [39] and Occam [41].

From a practical standpoint, each refined datatype is defined as a structure in Handel-C, while each process is implemented as a macro procedure. The constructs corresponding to the CSP stage are divided into 2 main categories for organisation purposes. The first category represents the definitions of the refined datatypes. The second category implements the refined processes. The refined processes are divided into different groups. The utility processes group contains macros responsible for producing, storing, sinking, broadcasting data and etc. The basic processes group contains macros that correspond to simple arithmetic and logical operations. These basic processes could be simple addition, multiplication, etc. The higher-order processes group contains the macros realising the CSP implementations corresponding to the higher-order functions. A separate group contains the macros that handle the FPGA card setup and general functionality. The reusable macros found in these groups serves as building blocks used for constructing a certain specified algorithm. This organisation is visualised in Figure 10.

3.4.1 Datatypes Definitions

The datatypes definitions are implemented using structures. This method supports recursive as well as simple types. The definition for an Item of a type Msgtype is a structure that contains a communicating channel of that type.

#define Item(Name, Msgtype) struct { chan Msgtype channel; Msgtype message; } Name

For generality in implementing processes the type of the communicating structure is to be determined at compile time. This is done using the typeof type operator, which allows the type of an object to be determined at compile time. For this reason, in each structure we declare a message variable of type Msgtype.

A stream of items, called StreamOfItems, is a structure with three declarations a communicating channel, an EOT channel, and a message variable [16]:

#define StreamOfItems(Name, Msgtype) struct { Msgtype message; chan Msgtype channel; chan Bool eotChannel; } Name

A vector of items, called VectorOfItems, is a structure with a variable message and another array of sub-structure elements [16].

#define VectorOfItems(Name, n, Msgtype) struct { struct { chan Msgtype channel; } elements[n]; Msgtype message; } Name

Other definitions are possible, but it affects the way a channel is called using the structure member operator (.). Examples of different extended definitions are as follows (the first definition reuses the Item structure, while the second one employs channel arrays supported in Handel-C):

#define VectorOfItems(Name, n, Msgtype) struct { struct { Item(element, MsgType); } elements[n]; } Name

#define VectorOfItems(Name, n, Msgtype) struct { chan Msgtype channel[n]; Msgtype messages; } Name

In general, there are certain limitations in the Handel-C language, which make the expression of a number of useful constructs either difficult or impossible. An example of an impossible to implement declaration is as follows:

StreamOfItems(Name, VectorOfItems(Name, n, Int16));

A simple preprocessor would facilitate a higher level of Handel-C generic definitions, and allow them to flow much more freely from our functional specifications. The implementation of such a preprocessor is being investigated within our research group.

3.4.2 Utilities Macros

The utility processes used in the implementation are related to the employed datatypes. The Handel-C implementation of these processes relies on their corresponding CSP implementation. In the following, we present an instance of these utility macros.

macro proc ProduceItem(Item, x){ Item.channel ! x;}

macro proc StoreItem(Item, x){ Item.channel ? x;}

3.4.3 Basic Processes Macros

This group of macros represents the fine-grained processes. A sample basic macro procedure Addition is included as an example.

macro proc Addition(xItem, yItem, output){ typeof (xItem.message) x,y; xItem.channel ? x; yItem.channel ? y; output.channel ! (x + y);}

3.4.4 Higher-Order Processes Macros

An example for an implementation in Handel-C of the CSP refinement of a higher-order function (map) is done as follows. The process hinges around a loop which terminates when the variable eot is set to true. At each step of the loop, the process enters a wait state until either the EOT or the value channel of the input stream is willing to communicate. If the EOT channel is willing to communicate, the input is consumed from it and stored in the variable eot, then output an EOT message for the output stream. If the value channel of the input stream is willing to communicate, the value is consumed then $F$ is applied to it giving the result on the output stream channel [16].

macro proc SMAP (streamin, streamout, F){ Bool eot; eot = False; do{ prialt{ case streamin.eotChannel ? eot: streamout.eotChannel ! True; break; default: F (streamin.elements,streamout.elements); break; }} while (!eot)}

We turn the attention to providing a definition in Handel-C for the behaviour of the process VMAP. Here we can employ Handel-C’s enumerated par construct to place $n$ instances of the process $F$ in parallel. Each instance is passed to the corresponding channels from both the input and output conduits [16].

macro proc VMAP (n, vectorin, vectorout, F) { typeof (n) c; par (c = 0 ; c < n ; c++){ F(vectorin.elements[c], vectorout.elements[c]);}}

The implementations of the stream and vector settings of the remaining high-order functions is done in a similar manner.

3.4.5 The RC-1000 System Control Macros

The Celoxica RC-x000 boards provide high-performance, real-time processing capabilities and are optimised for the Celoxica DK design suite. The RC-1000 is a standard PCI bus card, with four onboard banks of SRAM, equipped with a Xilinx Virtex with up to 2 million system gates[6].

According to the characteristics of the used system, some reusable macro procedures were implemented to be employed in the development model. For instance, reading or storing an Item from (in) a bank could be implemented as in the following macros:

macro proc ReadItemFromBank1(r){ Int temp;

PP1000ReadBank1(0, temp);
r.channel ! temp;}

macro proc StoreItemToBank1(r) { Int temp; unsigned int 21 count;

r.channel? temp;
PP1000WriteBank1(0, temp);}

3.5 Evaluation Tools and Performance Metrics

Different tools are used to measure the performance metrics used for the analysis. These tools include the design suite (DK) from Celoxica, where we get the number of NAND gates for the design as compiled to the Electronic Design Interchange Format (EDIF). The DK also affords the number of cycles taken by a design using its simulator. Accordingly, the speed of a design could be calculated depending on the expected maximum frequency of the design. The maximum frequency could be determined by the timing analyzer. Accordingly, the time to execute 1 cycle could be determined (Period). The execution time of a design is then the Period multiplied by the number of cycles. Thereat, the throughput is calculated depending on the amount of data processed in that execution time.

To get the practical execution time as observed from the host computer, the C++ high-precision performance counter is used. The counter probes the execution of the design after loading the image of the design into the FPGA till termination.

The information about the hardware area occupied by a design, i.e. number of Slices used after placing and routing the compiled code, is determined by the ISE place and route tool. Using the same tool we are able to get more detailed statistics about our compilation. In the current investigation the only used metrics are the number of Slices and the Total Equivalent Gate Count for a design.

4 A New Approach for developing a Matrix Multiplication Algorithm

Many parallel implementations of matrix multiplication have been investigated in the literature. Although this algorithm is simple and it has a long history, the continuous advancement in computer architectures made the study of this algorithm very interesting. Many matrix multiplication algorithms were suggested for parallel implementation. Horowitz and Zorat in [42] and Hake in [43] suggested a recursive divide-and-conquer solution. Fox et al. in [44] and Canon in [45] presented different ways this algorithm could be developed for a mesh topology. Other implementations were discussed in [46, 47, 48, 49].

An important requirement for any parallel hardware computation is the proper use of available computing resources. This is done either to minimise overall computation time, to minimise the chip area, or to compromise between these two goals. With the development model in hand, design flexibility is one of the main advantages granted. Accordingly, five refined designs from the functional specification of the standard matrix multiplication algorithm are presented. These designs vary between high-speed implementation with expensive use of resources, and lower speed implementation with less use of resources. The development of the matrix multiplication algorithm is presented in the following sections.

The development will start by formalizing the matrix multiplication algorithm. The functional specification will favor the use of the predefined high-order functions. This functional specification will be the source for different refinements with different degrees of parallelism described using CSP notation. The created parallel designs will be used in the section Section 7 as the basis of the code written in Handel-C.

5 Formal Functional Specification

An informal definition of the problem considers that the multiplication of two matrices ass and bss produces the matrix css whose elements, $c_{ij}(0\leq i<n,0\leq j<k)$ are computed as follows:

$c_{ij}=\sum_{t=0}^{m-1}a_{i,t}b_{t,j}$

where ass is an $n\times m$ matrix and bss is an $m\times k$ matrix. Items a and b correspond to elements from matrices ass and bss, respectively.

Generally, partitioning can be done very easily with matrix multiplication, where each matrix is divided into sub-matrices that can be manipulated as if they were a single matrix element [44]. This method is used to divide matrices with large dimensions to suit the expected limited capability of the available computer.

Turning our attention to the formalisation of the algorithm. A functional specification of matrix multiplication is formulated as a function mmult that takes two matrices as inputs and returns a matrix as a result. In this definition, we assume the first matrix is represented as a list of rows and the second matrix is represented as a list of columns.

mmult :: [[Int]] -> [[Int]] -> [[Int]] mmult ass bss = map (vmmult ass) bss

vmmult :: [[Int]] -> [Int] -> [Int] vmmult ass bs = map (scalarp bs) ass

scalarp :: [Int] -> [Int]-> Int scalarp as bs = sum (zipwithmul as bs)

sum :: [Int] -> Int sum rs = foldrl (+) rs

zipwithmul :: [Int] -> [Int] -> [Int] zipwithmul as bs = zipwith (*) as bs

The suggested algorithm for multiplying two matrices is done by mapping (using the higher-order function map) the function (vmmult ass) to all vectors in bss. This function is the multiplication of a vector with a matrix. It takes two inputs a matrix (list of lists) ass and a vector (list) bs and returns a list cs (a column in the resulting matrix).

In turn, vmmult maps the function (scalarp bs) over the list of lists ass. The function scalarp defines the scalar product of two vectors. The inputs to this function are two lists as and bs. The higher-order function zipWith is used in the function zipwithmul to zip the inputs with multiplication, then the function sum is employed to fold the already zipped lists with addition. The output of this composition is an element from the resultant matrix.

According to this specification, the implementation under HUGs98 Haskell compiler is tested at the unit, component and integration levels.

6 Algorithm Refinements

Clearly, parallelism is not a part of the starting specification of the stated problem. Typically, parallelism and communications are introduced at this stage in the development for the sole purpose of capturing functionally equivalent, but parallel, designs.

Applying the provably correct refinement rules, the previous specification is refined to CSP as a middle stage towards hardware realisation. The capability of doing different data refinements implies the availability of various designs. Whereby, each design would have different characteristics and levels of parallelism.

In the following refinements, five designs are presented. The first design is a data-parallel design, while, the second design is a stream-based design. The third and fourth designs addresses refinement to pipelined parallelism using function decomposition strategy. The last design is a 2D pipelined design with an extension leading to a systolic implementation of the problem.

For more clarification of the used terms we recall the following informal definitions. A data-parallel design replicates the same processes in order to compute for different data inputs. Commonly, a Single Program Multiple Data (SPMD) approach is used in data-parallel models, where data are distributed across processors. In pipelined computations, a program is divided into a series of tasks that have to be completed consecutively. Accordingly, these tasks are executed by separate pipeline stages. The pipe stages stream data from stage to stage to form the required computation. A stream-based design eliminates some replication (in data-parallel designs) or some stages (in pipelined designs) and it processes the input and output as streams of data. Systolic arrays are another parallel computing architecture. It is best described by analogy with the regular pumping of blood by the heart. In systolic arrays, processors are arranged in an array where data flow across the array elements between neighbours. For instance, a process firstly takes in data from one or more neighbours (North and West). Secondly, the process manipulates the input data. Finally, the process outputs results in the opposite direction (South and East).

6.1 First Design - Data Parallelism

Recalling the high-level specification of mmult:

mmult :: [[Int]] -> [[Int]] -> [[Int]] mmult ass bss = map(vmmult ass) bss

In this design we consider the refinement of the input bss and the output css as vectors of items of size k, where each item is a list.

$mmult(ass)::\lfloor[Int]\rfloor_{k}\rightarrow\lfloor[Int]\rfloor_{k}$

The CSP implementation of the functional specification of the matrix multiplication algorithm mmult realises this function as a process MMULT. The CSP process MMULT is the parallel execution of k-copies of VMMULT (the refinement of vmmult). This description is implemented using VMAP the vector setting of the higher-order function map. Therefore, it is the interleaving with renaming of the process VMMULT for all columns of bss:

$mmult(ass)\sqsubseteq MMULT(ass)$ $MMULT(ass)=VMAP_{k}(VMMULT(ass))$

The list ass is passed as an argument to each of the processes VMMUL(ass) in the above design. This design can be pictured as in Figure 11.

The list ass could be explicitly passed to the process VMMULT by exploiting the following algebraic identity:

$VMMULT(ass)=PRD(ass)\rhd VMMULT$

The effect of applying this step to the previous design is visualised in Figure 12. In this version, the list ass is locally produced and fed to each process VMMULT in the vector. The effect of having k parallel copies of PRD(ass) communicating with k instances of VMMULT can be achieved by factorising the process PRD(ass) and broadcasting its output to the relevant processing elements in the network. Applying this rule will result in a semantically equivalent version of MMULT which has a different layout, this is shown in Figure 13.

Now we turn our attention to the refinement of the function vmmult.

vmmult::[Int] -> [[Int]] -> [Int] vmmult bs ass = map (scalarp bs) ass

Clearly, vmmult (bs) is a map pattern. Since map has two different implementations, we will consider them in turn in this and the next designs. In this design, the CSP implementation realises the function vmmult as a process VMMULT with a vector of items ass=$\lfloor\textit{as}_{1},\textit{as}_{2},...\textit{as}_{n}\rfloor$ as input and a vector of items cs=$\lfloor\textit{c}_{1},\textit{c}_{2},...\textit{c}_{n}\rfloor$ as output. By refining scalarp into VSCALARP, the CSP implementation of vmmult (bs) is again the off the shelf refinement of a vector map:

$vmmult(bs)::\lfloor[Int]\rfloor_{m}\rightarrow[Int]_{n}$

$VMMULT(bs)=VMAP_{n}(VSCALARP(bs))$

By appealing to the same technique already used in the refinement of VMMULT we get:

$VSCALARP(bs)=PRD(bs)\rhd VSCALARP$

This leads to a new design ofVMMULT(bs):

$VMAP_{n}(PRD(bs)\rhd VSCALARP)$

$BROADCAST_{n}(bs)\rhd VMAP_{n}(VSCALARP(bs))$

Figure 14 shows the process VMMULT. This step also shows clearly the replication of the process VSCALARP, which is an indicator for a later replication in the hardware implementation.

Figure 15 expands the main building block in Figure 13 by corresponding configuration in Figure 14. This gives a two dimensional visualisation of the process MMULT as a data parallel implementation.

The next step is to present the building block VSCALARP corresponding to a CSP refinement of the function:

scalarp :: [Int] -> [Int] -> Int scalarp as bs = sum (zipwithmul as bs)

This function can be refined as the piping of two processes VZIPWITH and VFOLD corresponding to refinements of the functions zipwithmul and sum respectively.

$scalarp::\lfloor Int\rfloor_{m}\rightarrow\lfloor Int\rfloor_{m}\rightarrow Int$

$VSCALARP=VZIPWITH_{m}(MUL)>>_{m}VFOLD_{m}(ADD)$

This description is depicted in Figure 16.

For completeness, the CSP implementations of the simple addition and multiplication functions are:

$ADD=(in_{1}?a\rightarrow SKIP\mathrel{\mathord{\mid}\mathord{\mid}\mathord{\mid}}in_{2}?b\rightarrow SKIP)\mathchar 24635\relax\;out!(a+b)$

$MUL=(in_{1}?a\rightarrow SKIP\mathrel{\mathord{\mid}\mathord{\mid}\mathord{\mid}}in_{2}?b\rightarrow SKIP)$;$out!(a\times b)$

6.2 Second Design - Streaming I/O

For this design the refinement of mmult is not changed. The change will appear in the refinement of the function vmmult. Recall the formal specification of this function:

vmmult::[Int] -> [[Int]] -> [Int] vmmult bs ass = map (scalarp bs) ass

At this point, the input list ass is viewed as a stream of values ass= $\langle\textit{as}_{1},\textit{as}_{2},...\textit{as}_{n}\rangle$ and the output list as a stream of values as well. The CSP refinement of* vmmult(bs)* is directly obtained from the off-the-shelf stream-based implementation of the higher-order function map:

$vmmult(bs)::\langle[Int]\rangle\rightarrow\langle Int\rangle$

$VMMULT(bs)=MAP(VSCALARP(bs))$

Figure 17 shows the new version of the process VMMULT. This step also shows clearly that there is no more replication of the process* VSCALARP*, which is an indicator for the later reduction in use of hardware resources.

Keeping the refinement of the remaining functions the same, MMULT process looks as in Figure 18.

6.3 Third Design - Pipelining

Demonstrating the refinement to pipelined parallelism is the purpose of this design. Generally, this kind of parallelism is a very effective means for achieving efficiency in numerous algorithms. Usually, pipelined parallelism is much harder to detect than data parallelism. Accordingly, the function decomposition strategy, found in [40] and recalled in Section 3.3 is used. This strategy aims at exhibiting pipelined parallelism in functional programs. According to the decomposition rule the definition of the function mmult is pipelined.

mmult ass bss = map (vmmult ass) bss

vmmult [] bs = [] vmmult (as : ass) bs = vscalarp as bs : vmmult ass bs

The recursive function in this case is vmmult. The value to be passed to the next stage of the pipe is a tuple. Its first is the input vector and its second is result of applying vscalarp on the input vector from matrix bss and the present argument from matrix ass. According to the decomposition rule, the efficient implementation of mmult as a pipelined network of CSP processes can be as follows:

$MMULT=PRD(bss)\rhd((\gg)/(MAP\circ f^{\prime})\ast[as_{n},as_{n-1},...,as_{0}])$

$MAP(initial)=\mathop{\mathstrut{\mathchar 278\relax}}\nolimits Z$  $\bullet$    $left?eot\rightarrow right!eot\rightarrow SKIP$

$|$

$left?x\rightarrow right!\langle bs,[]\rangle\rightarrow Z$

$MAP(f^{\prime}$  $as)=\mathop{\mathstrut{\mathchar 278\relax}}\nolimits Z$  $\bullet$   $left?eot\rightarrow right!eot\rightarrow SKIP$

$|$

$left?\langle bs,y\rangle\rightarrow right!\langle bs,y{++}(vscalarp$  $as$  $bs)\rangle\rightarrow Z$

$MAP(final)=\mathop{\mathstrut{\mathchar 278\relax}}\nolimits Z$  $\bullet$    $left?eot\rightarrow right!eot\rightarrow SKIP$

$|$

$left?\langle bs,y\rangle\rightarrow right!y\rightarrow Z$

The decomposed pipelined network is shown in Figure 19. In this design, the matrix bss is input to the network as a stream of vectors (columns) $\langle\textit{bs}_{1},\textit{bs}_{2},...\textit{bs}_{k}\rangle$. The matrix ass vectors (row by row) are produced in the pipe stages. The result is considered as a stream of streams $\langle\textit{cs}_{1},\textit{cs}_{2},...\textit{cs}_{k}\rangle$ The first result to appear from the network is the output stream (column) $\langle\textit{cs}_{k}\rangle$ corresponding to the first input vector bsk. This design is independent from the size of k a dimension of bss and css.

6.4 Fourth Design - Pipelined Turnout Stages

This design makes use of an optimisation of the previous design. In this case, the input matrix* bss* is refined as a stream of vectors $\langle\textit{bs}_{1},\textit{bs}_{2},...\textit{bs}_{k}\rangle$ and the matrix ass is refined as arguments in the pipeline stages. The output from each stage is turned out as a result, besides forwarding the input from bss to the next stage. Thus, the output from this pipeline is a vector of streams as shown in Figure 20. Note that, this design also doesn’t depend on the size of k - the dimension in bss and css. The CSP implementation is as follows:

$MMULT=PRD(Bss)((\gg)/(MAP\circ f^{\prime})\ast[as_{n},as_{n-1},...,as_{0}])\mathrel{\mathord{\parallel}}SINK$

$MAP(f^{\prime}$  $as)=\mathop{\mathstrut{\mathchar 278\relax}}\nolimits Z$  $\bullet$    $left?eot\rightarrow right!eot\rightarrow SKIP$

$|$

$left?bs\rightarrow down!(vscalarp$  $as$  $bs)\rightarrow right!bs\rightarrow Z$

6.5 Multilevel Pipelined Design

This design applies the function decomposition strategy for pipelined parallelism on two levels. The first level is decomposing the vector matrix multiplication into a pipeline of processes performing the scalar product of two vectors, this is similar to the fourth design. An addition is made by pipelining the scalar product routine creating a second level pipeline. The final structure of the suggested design is multilevel pipelines realising the matrix multiplication algorithm. The decomposition of the process VSCALARP is done in a similar manner (See Figure 21):

$VSCALARP=PRD(0)\rhd((\gg)/(MAP\circ f^{\prime})\ast[a_{m},a_{m-1},...,a_{0}])$

$MAP(f^{\prime}$  $a)=\mathop{\mathstrut{\mathchar 278\relax}}\nolimits Z$  $\bullet$    $left?eot\rightarrow right!eot\rightarrow SKIP$

$|$

$left?x\rightarrow up?b\rightarrow right!(x+(a\times b))\rightarrow Z$

Thus MMULT implementation is as follows:

$MMULT=PRD(Bss)((\gg)/(MAP\circ f^{\prime})\ast[as_{n},as_{n-1},...,as_{0}])\mathrel{\mathord{\parallel}}SINK$

$MAP(f^{\prime}$  $as)=\mathop{\mathstrut{\mathchar 278\relax}}\nolimits Z$  $\bullet$    $left?eot\rightarrow right!eot\rightarrow SKIP$

$|$

$left?bs\rightarrow(PRD(bs)\rhd VSCALARP(as)[down/right])\mathchar 24635\relax\;$

$right?bs\rightarrow Z$

An optimisation of this design would lead to a systolic solution. The main idea of the change is to enable the communication between parallel VSCALARP stages in MMULT. A VSCALARP is to be the parallel execution of basic cells, each responsible of forwarding down the upper input from bss. The cell is also responsible for doing its part for the scalar product computation. This part is done by outputting to the right the result of adding the left input to the multiplication of the upper input b (from bss) with the item a (from ass). The basic process is called CELL (See Figure 22) and defined as:

$CELL(a)=PRD(a)\rhd(up?u\rightarrow left?l\rightarrow right!(u*a+l)\rightarrow down!u)$

Then, VSCALARP is implemented as:

$VSCALARP(as)=PRD(0)\rhd$

$(\mathrel{\mathord{\parallel}}_{i=1}^{i=m}CELL(as[i])[d_{i}/left,d_{(i+1)}/right,e_{i}/up,e_{(i+1)}/down])$

Thus, MMULT implementation is as follows:

$MMULT=PRD(Bss)((\gg)/(MAP\circ f^{\prime})\ast[as_{n},as_{n-1},...,as_{0}])\mathrel{\mathord{\parallel}}SINK$

$MAP(f^{\prime}$  $as)=\mathop{\mathstrut{\mathchar 278\relax}}\nolimits Z$  $\bullet$    $up?eot\rightarrow down!eot\rightarrow SKIP$

$|$

$VSCALARP(as)[right/d_{m}]\mathchar 24635\relax\;Z$

This implementation is depicted in Figure 23.

For the sake of giving a similar design, an intuitive CSP implementation of the matrix multiplication algorithm shown in Figure 23 is as follows:

$CELL_{(i,j)}(ass[i,j])=\mathop{\mathstrut{\mathchar 278\relax}}\nolimits X\bullet PRD(ass[i,j])\rhd((up?u\rightarrow down!u\rightarrow(SKIP\lessdot u=eot\gtrdot left?l\rightarrow right!(u*ass[i,j]+l))\rightarrow SKIP\rightarrow X)$

The matrix multiplication process MMULT is then implemented as:

$MMULT=BROADCAST_{n}(0)[d/out]\rhd$

$(\mathrel{\mathord{\parallel}}_{i=1}^{i=n}(\mathrel{\mathord{\parallel}}_{j=1}^{j=m}CELL(i,j)[b_{ij}/lwft,b_{i(j+1)}/right,e_{ij}/up,e_{(i+1)j}/down]))$

Finally, these designs depend only on the dimensions n and m from ass.

7 Reconfigurable Hardware Implementation

As a stage in the development model, Handel-C code follows the refined CSP implementation of the presented designs. The targeted circuit implementation is to have the same topology of the communicating processes shown in the refinement section. In the following subsections, some pieces of code taken from the five different designs are presented.

7.1 First Design - Data Parallelism

From first design, we recall the CSP implementation of the process MMULT:

$MMULT(ass)=VMAP_{K}(VMMULT(ass))$

The code corresponding to the process MMULT is the macro MatrixMult with bss as an input and css as an output. This macro is vector mapping of the macro VectMatrixMult, where the vector of vectors ass is internally produced. Recall that the macro VMap is the implementation of the higher-order process VMAP. In this case, VMap works by distributing the matrix bss for each process VectMatrixMult, where a vector from css will be the output.

macro proc MatrixMult (bss, css, n){ VMap(bss, css, n, VectMatrixMult);}

Then, the macro VectMatrixMult implements the process VMMULT, refined as:

$VMMULT(ass)=PRD(ass)\rhd VMMULT$

This applies again VMap calling the macro VScalarP (the implementation of the process VSCALARP) for each vector in ass. At his point, bss is sinked as it will be later internally produced in VScalarP.

macro proc VectMatrixMult (bss, cs) { VectorOfVectorsOfItems(ass, n, m, Int); SinkVectorOfVectorsOfItems(bss, n, m, sink); par { ProduceVectorOfVectorsOfItems(ass, n, m, tempAss); VMap(ass, cs, n, VScalarP);}}

The macro VScalarP implements the CSP process VSCALARP:

$VSCALARP(bs)=PRD(bs)\rhd VSCALARP$

The internal production for a vector from bss is done according to an index. This allows producing a different vector bs from bss for each execution of VScalarP.

macro proc VScalarP(as, outputItem, index){ VectorOfItems(internalV, m, Int); VectorOfItems(bs, m, Int); par { ProduceVectorOfItems(bs, m, tempBss[index]); VZipWith(m, as, bs, internalV, Multiplication); VFoldR (internalV, outputItem, m, Addition, 0);}}

The macros Addition and Multiplication corresponding to the processes ADD and MUL are implemented as:

macro proc Multiplication(xItem, yItem, output) { Int x,y; xItem.Channel ? x; yItem.Channel ? y; output.Channel ! (x * y);}

macro proc Addition(xItem, yItem, output) { Int x,y; xItem.Channel ? x; yItem.Channel ? y; output.Channel ! (x + y);}

Running the above implementation is done practically by producing bss and storing css from/to a buffer. The RC-1000 board internal SRAMS are used as the input and output buffers. The main macro call that runs the code implementing the first design for the matrix multiplication algorithm is as follows:

LoadVectorOfVectorsOfItemsFromBank0(n, m, ass);

par{ ProduceVectorOfVectorsOfItems(bss, m, k, bssTemp); MatrixMult(bss, css, m); StoreVectorOfVectorsOfItems(css, n, k, c ssTemp);}

7.2 Second Design - Streaming I/O

The code implementation of the current design reflects the change as applied to the first design in the refinement to CSP change. The process VMMULT definition is recalled:

$VMMULT(ass)=PRD(ass)\rhd VMMULT$

The only modification to the implementation is done by refining ass to a stream of vectors produced internally within VectMatrixMult. To meet the change, this macro employs the macro SMap the sequential implementation of the process SMAP. Accordingly, the Handel-C code is:

macro proc VectMatrixMult (bss, cs) { StreamOfVectorsOfItems(ass, n, m, Int); SinkVectorOfVectorsOfItems(bss, n, m, sink); par { ProduceStreamOfVectorsOfItems(Ass, n, m, tempAss); SMap(ass, cs, n, VScalarP);}}

7.3 Third Design - Pipelining

For this design, the decomposed process for a single pipe stage is implemented as the macro PipeStage. This macro starts by inputting from the left a vector bs and a stream cs. These two inputs are the output of the left identical pipe stage. For the initial pipe stage one input is the produced stream of vectors bss and nothing on the stream channel. The PipeStage then computes for the scalar product, forward the new stream cs and the vector bs to the right identical process and finally waits for another input signal. When the end of transmission is reached, an EOT signal is sent to the right pipe stage. The output bss from the final stage is sinked, while the stream of streams output css is stored as the result. The code of a single pipe stage is as follows:

macro proc PipeStage(tupleIn, tupleOut, iAs){ . . . Item(outputItem, Int16); VectorOfItems(tempBs, m, Int16); VectorOfItems(tempAs, m, Int16);

do{
    prialt{
    case tupleIn.element1.elements[0].channel ? tempVbs[0]:
        par(j = 1; j < m; j++){
         sOVIn.elements[j].channel ? tempVbs[j];}
         StoreStreamOfItems(tupleIn.element2.elements, iAs, tCs);
        par{
         ProduceVectorOfItems(tempBs, m, tempVbs);
         ProduceVectorOfItems(tempAs, m, ass[iAs]);
         VScalarP(tempAs, tempBs, m, outputItem);
         StoreItem(outputItem, tCs[iAs]);}

        ProduceVectorOfItems(tupleOut.element1, m, tempVbs);
        ProduceStreamOfItems(tupleOut.element2, (iAs + 1), tCs);
        break;

    case TupleIn.element1.eotChannel ? eot:
        TupleIn.element2.eotChannel ? eot1;
        TupleOut.element1.eotChannel ! True;
        TupleIn.element2.eotChannel ! True;
        break;}} while (!eot);}

The general replicating macro that corresponds to the employed decomposition is implemented as in the following code section. This macro takes the advantage of using the ifselect and par Handel-C constructs. By using ifselect, whole statements can be selected or discarded at compile time, depending on the evaluation of the expression. Accordingly, the par statement selects only one macro execution for P according to the value of c for each replication. The parameter c corresponds to the pipe stage number. In this design this parameter is initialised to 1 instead of 0 since a different initial pipe stage is implemented to overcome a limitation in Handel-C. This limitation forbids the production of stream of streams with a size 0 needed for the first pipe stage. The final picture of this implementation is best depicted as in Figure 19.

macro proc Pipe (tIn, tOut, n, P){ typeof(tIn) cmids[n - 1];

    par (c = 1; c < n; c++){
        ifselect (c == 1)
           P(tIn, cmids[c], c);
        else ifselect (c < n - 1)
            P(cmids[c - 1], cmids[c], c);
        else
            P(cmids[c - 1], TOut, c);}}

The execution of the matrix multiplication is done as follows.

par { ProduceStreamOfVectorsOfItemsFromBank0(bssIn, m, bss); PipeStageInitial(bssIn, tIn, 0); Pipe (tIn, tOut, n, PipeStage); StoreStreamOfStreamsOfItemsInBank1(tOut.element2, n, tempOut1); SinkStreamOfVectorsOfItems(tOut.element1, m, tempOut2);}

7.4 Fourth Design - Pipelined Turnout Stages

In a similar way of implementing the third design, this design uses a modified version of the previous macros. In the new TurnoutPipeStage the output is a vector of streams vOSOUT. Each pipe stage outputs its own stream and issues its own termination signal through its own eot channel. This macro works by firstly inputting a vector from the stream bss, forwards it to the right stage, then computes for the scalar product turning out its result. These steps are repeated till the end of transmission of bss. At that point, local termination signals are issued from the stage. The final picture of this implementation is best depicted as in Figure 20.

macro proc TurnoutPipeStage(sOVIn, sOVOut, vOSOut, indexForAs) { . . . VectorOfItems(tempBs, m, Int16); VectorOfItems(tempAs, m, Int16); do{ prialt{ case sOVIn.elements[0].channel ? tempVbs[0]: par(j = 1; j < m; j++){ sOVIn.elements[j].channel ? tempVbs[j];} ProduceVectorOfItems(sOVOut, m, tempVbs); par{ ProduceVectorOfItems(tempBs, m, tempVbs); ProduceVectorOfItems(tempAs, m, ass[indexForAs]); VScalarP(tempAs, tempBs, m, vOSOut);} break; case sOVIn.eotChannel ? eot: sOVOut.eotChannel ! True; vOSOut.eotChannel ! True; break;}} while (!EOT);}

This pipe pattern is implemented as:

macro proc TurnoutPipe(in1, out1, out2, n, p) { typeof(in1) cmids1[n + 1];

par (c = 0; c < n; c++){
    ifselect (c == 0)
        p (in1, cmids1[c], out2.elements[c], c);
    else ifselect (c < n - 1)
        p (cmids1[c - 1], cmids1[c], out2.elements[c], c);
    else
        p(cmids1[c - 1], out1, out2.elements[c], c);}}

The execution of this design’s implementation is done as follows:

void main(void) { . . . StreamOfVectorsOfItems(bssIn, m, Int16); StreamOfVectorsOfItems(bssOut, m, Int16); VectorOfStreamsOfItems(cssOut, n, Int16); par { ProduceStreamOfVectorsOfItems(bssIn, n, m, ass); TurnOutPipe(bssIn, bssOut, cssOut, 3, TurnoutPipeStage); StoreVectorOfStreamsOfItems(cssOut, n, k, tempCss); SinkStreamOfVectorsOfItems(bssOut, m, tempBss);} StoreMatrixInBank2(tempBss, n, k);}

It is clear from this design that n-parallel output streams are employed. Consequently, a buffer with multi-concurrent access is needed. This kind of access is not allowed with the available single R/W onboard SRAMs. Besides, the number of available banks for concurrent access is only 4. This introduces a limitation to the practical implementation of this design:

The parameter k will appear as a static constant in the compilation. Thus, for any new stream with a certain length a new compilation is needed.

• 2.

The limited ability of the FPGA to store the results on its internal area, especially, for matrices with large dimensions.

the suggested general solution to this problem is storing the resultant VectorOfStreamsOfItems on the local FPGA memory, and then storing them back as a stream of values or a vector of streams with up to 4 parallel streams only.

7.5 Multilevel Pipelines Design

This design implementation uses two kinds of pipe replicating macros for the two needed pipelined levels. The first pipe implements the decomposition for the vector matrix multiplication process. Thus, we reuse the macro TurnoutPipe employed in the fourth design. Furthermore, the second pipe macro, called SystolicPipe, implements the scalar product process. This macro replicates pipe stages with two inputs and two outputs. Recall the CSP implementation for a the basic computation cell:

$CELL(a)=PRD(a)\rhd(up?u\rightarrow left?l\rightarrow right!(u*a+l)\rightarrow down!u)$

The code corresponding to CELL is:

macro proc SystolicPipeCell(li, upV, r, downV, i, j) { . . . Item(temp, Int16);

upV.Channel ? tempb;
result = tempb * ass[i][j];
par{
    Addition(temp, li, r);
    ProduceItem(temp, result);}
downV.Channel ! tempb;}

To implement the complete algorithm, these cells are composed to form a scalar product pipeline using the macro SystolicPipe.

macro proc SystolicPipe (in1l, in2up, out1r, out2d, n, i,P){ typeof(in1l) cmids1[n + 1];

par (c = 0; c < n; c++){ ifselect (c == 0) P(in1l,in2up.elements[c],cmids1[c],out2d.elements[c],i,c); else ifselect (c < n - 1) P(cmids1[c-1],in2up.elements[c],cmids1[c],out2d.elements[c],i,c); else P(cmids1[c-1],in2up.elements[c],out1r,out2d.elements[c],i,c);}}

This scalar product pipe is the core of the pipe stage needed for the pipeline implementing the matrix multiplication algorithm. A pipe stage is implemented as:

macro proc PipeStage(bssIn, bssOut, cssOut, i){ VectorOfItems(tempBs, m, Int16); Item(rO, Int16); Item(lI, Int16); . . . do{ prialt{ case bssIn.elements[0].channel ? temp[0]: par(j = 1; j < m; j++){ bssIn.elements[j].channel ? temp[j];} par{ ProduceItem(lI, 0); ProduceVectorOfItems(tempBs, m, temp); SystolicPipe(lI, tempBs, rO, cssOut, m, i, SystolicPipeCell); StoreItem(rO, tempItem);} break;

    case bssIn.eotChannel ? eot:
      bssOut.eotChannel ! eot;
      break;}} while (!eot);}

This pipe stage is then replicated to form a pipeline using the predefined macro TurnoutPipe. The main code section running the above implementation is:

void main(void) { . . . StreamOfVectorsOfItems(bssIn, m, Int16); StreamOfVectorsOfItems(bssOut, m, Int16); StreamOfVectorsOfItems(cssFinal, n, Int16);

par{
    ProduceStreamOfVectorsOfItemsFromBank1(bssIn, m);
    {TurnoutPipe(bssIn, bssOut, cssFinal, n, PipeStage);
    cssFinal.eotChannel !True;}
    SinkStreamOfVectorsOfItemsToBank3(bssOut, m);
    StoreStreamOfVectorsOfItemsInBank2(cssFinal,  n);}}

8 Performance Analysis and Evaluation

The development is originated from a specification stage, whose main key feature is its powerful higher-level of abstraction. During the specification, the isolation from parallel hardware implementation technicalities allowed for deep concentration on the specification details. Whereby, for the most part, the style of specification comes out in favor of using higher-order functions. Two other inherent advantages for using the functional paradigm are clarity and conciseness of the specification. This was reflected throughout all the presented studies. At this level of development, the correctness of the specification is insured by construction from the used correct building blocks. The implementation of the formalised specification is tested under Haskell by performing random tests for every level of the specification.

The correctness will be carried forward to the next stage of development by applying the provably correct rules of refinement. The available pool of refinement formal rules enables a high degree of flexibility in creating parallel designs. This includes the capacity to divide a problem into completely independent parts that can be executed simultaneously (pleasantly parallel). Conversely, in a nearly pleasantly parallel manner, the computations might require results to be distributed, collected and combined in some way. Remember at this point, that the refinement steps are systematic and done by combining off-the-shelf reusable instances of basic building blocks.

In this case study, we will measure the speed in Items per Second (ips). For instance, a $3\times 3$ matrix has 9 items.

In Table 1 the results for running the different designs are presented. The first design occupied an area of 564 Slices per Item running at a speed of 257.14 Kips for a network with dimensions of $(3\times 3\times 3)$. The second design occupied a smaller area, as expected, but the realised design runs with a speed of 280 Kips. The pipelined second and third designs achieved a better speed of 2.25 Mips with less areas of 227 and 237 Slices per Item. The smallest area ratio of 158.7 Slices per Item has been occupied when placing and routing the 2D pipelined design. The speed achieved by this design is 3.1 Mips.

The second design is found to be $8.9\%$ faster than the first design, also the Slices to Items ratio of the second design is 42.7$\%$ less than that of the first design. Thus, the second design can accommodate for a larger number of items with a better speed as compared to the first design. The modification done to the third pipelined design yielding the fourth pipelined design shows that there were no effect on the speed of execution. However, the fourth design occupied a 4.4$\%$ larger area. Thus, the modification didn’t leave a positive effect on the performance. The 2D pipelined design has shown a better performance than the other designs, for instance, it occupies a 30.1$\%$ less Slices per Item than the third pipelined design, also achieving a 38$\%$ higher speed.

The $(11\times 11)$ 2D pipelined cells design is independent from the third dimension $k$. In Table 2 we compare the execution time of running for different values of $k$ between the RC-1000 and two computer machines. These are a 1.2 GHz Athlon AMD machine with 512MB of RAM, and a 1.4 GHz P4 with 1GB of RAM. It is shown in Figures 25 and  24 that these machines will outperform the suggested design implementation on the RC-1000 when the value of $k$ is nearly at a value of $299$ items. We note here the possible effect of the bus connecting the memory and FPGA on the speed of execution. To cope with this limitation a suggestion could be proposed for adding a cache memory to handle the input and output streams of data.

9 Acknowledgement

I would like to thank Dr. Ali Abdallah, Prof. Mark Josephs, Prof. Wayne Luk, Dr. Sylvia Jennings, and Dr. John Hawkins for their insightful comments on the research which is partly presented in this paper.

10 Conclusion

Mapping parallel versions of algorithms onto hardware could enormously improve computational efficiency. Recent advances in the area of reconfigurable computing came in the form of FPGAs and their high-level HDLs such as Handel-C. In this paper, we build on these recent technological advances by presenting, demonstrating and examining a systematic approach of behavioural synthesis. This system creates a functional specification of an algorithm without defining parallelism. Correspondingly, an efficient parallel implementation is derived in the form of CSP network of processes. Accordingly, we create efficient parallel implementations in Handel-C. The presented work included theory and practices about the suggested methodology. This paper also presented a demonstration for using a proposed model to synthesise reconfigurable hardware for the matrix multiplication algorithm. The general functional specification is discussed firstly followed by the provably correct step-wise refitment to CSP. Many possible designs were engineered and compiled to reconfigurable hardware with different levels of parallelism. The hardware implementation using Handel-C is shown stressing the correspondence to the CSP refined networks. To complete the synthesis, these designs were compiled to EDIF format and then placed and routed. Accordingly, a performance study has been included for the realised designs. The first design required the largest area with respect to the number of items used, that is 564 Slices per Item for a speed of 257.14 Kips. The modification to the first design which lead to the second design helped in reducing the area to 323 Slices per Item for a speed of 280 Kips. The third and fourth pipelined designs occupied areas of 227 and 237 Slices per Item for a speed of 2.25 Mips. The lastly realised 2D pipelines design has the best area to items ratio of 158.7 Slices per Item and a speed of 3.1 Mips. Future work includes extending the theoretical pool of rules for refinement, the investigation of automating the development processes, and the optimisation of the realisation for more economical implementations with higher throughput.

11 Biography for Issam Damaj

Issam W. Damaj (Ph.D. M.Eng. B.Eng. MIEEE MIEE) received his B.Eng. in Computer Engineering from Beirut Arab University in 1999 (with high distinction), and his M.Eng. in Computer and Communications Engineering from the American University of Beirut in 2001 (with high distinction). He was awarded his Ph.D. degree in Computer Science from London South Bank University, London, United Kingdom in 2004. Currently, he is with the Electrical and Computer Engineering Department at Hariri Canadian Academy for Sciences and Technology, Lebanon. His research interests include reconfigurable computing, parallel processing, h.w./s.w. co-design, computer interfacing and applications, fuzzy logic, and computer security. He has more than 25 international and regional research publications and projects. He is a Member of the IEEE and IEE professional organizations, and the order of Engineers in Beirut.