Ecaterina Madalina
email: madalina.andreica@gmail.com

Andreica

Ionut Mugurel

Nicolae Cataniciu
email: ncataniciu@incsmps.ro

Multidimensional Data Structures and Techniques for Efficient Decision Making

Keywords: Decision making, Multidimensional data structures, Risk management, Range aggregation query, Range selection query, Range minimum query, Weighted median

In this paper we present several novel efficient techniques and multidimensional data structures which can improve the decision making process in many domains. We consider online range aggregation, range selection and range weighted median queries; for most of them, the presented data structures and techniques can provide answers in polylogarithmic time. The presented results have applications in many business and economic scenarios, some of which are described in detail in the paper.

Introduction

The process of decision making is both a permanent necessity and a challenge in many economic and industrial fields, like risk management, banking, selection of financing means, such as leasing or credit, bonity analysis of clients, e-commerce, operational research, and many others. The importance of the decision making process has been confirmed by the large number of publications which develop and propose efficient decision making techniques. These techniques can be classified into several broad categories, such as those handling certain and complete information, those using uncertain data and objectives, and those based on risk assessment. Some of the best known decision making optimization methods are multi-attribute and multi-objective decision making [START_REF] Resteanu | Distributed and Parallel Computing in MADM Domain using the OPTCHOICE Software[END_REF], fuzzy decision rules [START_REF] Stoica | Fuzzy Sets and Their Applications[END_REF] and dynamic programming [START_REF] Ben-Ameur | Dynamic Programming Approach for Valuing Options in the GARCH Model[END_REF].

In order for the decision making process to obtain significant results, we need two factors: a good optimization technique and accurate input values. Although a lot of effort has been directed towards developing highly efficient decision making optimization methods and models, the process of obtaining accurate input information as quickly as possible seems to have been mostly overlooked. This paper presents several techniques and data structures for obtaining aggregate information efficiently from a large database, which can later be used as input data in a decision process.

The main scenario in which we consider the implementation of the presented techniques is the following: A large database, modeled as a multidimensional data cube, is available on a central server, which can be accessed by multiple clients (e.g. economic agents). Every client can ask at any time for the computation of some aggregate information (e.g. min, max, sum, average, median) over several (large) portions of the database. Different clients may be interested in different parts of the database. Using standard techniques in such cases would lead to unacceptably high response times, especially in the case of many clients accessing the database simultaneously. On the other hand, the techniques we propose in this paper lead to response times which are several orders of magnitude better than the common methods. The considered scenario occurs, for instance, in the case of large companies or banks which have all of their data stored in a data center which is accessed by their geographically spread headquarters and subsidiaries.

The rest of this paper is structured as follows. Section 2 presents several use cases for the implementation of the techniques and data structures which are introduced in Section 3. In Section 4 we discuss related work and we conclude.

Use Cases

In this section we present several situations where our techniques can be used successfully. Let's consider the case of a large retail company which handles thousands of transactions per day (both in department stores and over the Internet). The information associated to each transaction (e.g. date and time, purchase price, quantity, personal information about the customer, customer's answers to relevant surveys) is carefully stored in a central database. Periodically, the company's managers have to decide upon the price, quantity and (types of) products which should keep being sold or which should be taken off the market or improved, or if any new products should be launched. In this decision making they need quick access to aggregate information like: the total amount of sales of certain types of products to customers of a certain age and income range; the degree of satisfaction of each consumer target group which acquired products with certain characteristics; the percentage of the sales of each product from the total amount of sales; the most likely characteristics (age and income) of the highest paying consumers from a given age and income range.

Let's consider now a financial consulting company (e.g. specialized on risk assessment and management) which owns a private database with information regarding the outcomes of several strategic decisions made by multiple companies over a large period of time. A company is characterized by multiple parameters, such as total income, profit, investments, financial sources, market share, supply characteristics on the market, and so on. A strategic decision is also characterized by several specific parameters. When advising a client, the consulting company may compute an aggregate of the outcomes of the strategic decisions made by the competing companies with similar characteristics to the customer's company (i.e. using intervals containing the parameter values of the customer's company).

Another example consists of a bank which needs to perform a bonity analysis of its clients soliciting different credits. This is a crucial decision making process, especially in the context of the present global financial crisis. This is why the banks should use as much available data as possible and obtain the results in the most efficient manner.

On the other hand, there are several situations when the access to available data should be significantly improved. For instance, let's consider the case of a leasing company that wants to access some aggregate information regarding the financial status of several competing leasing companies on the market during a specified period of time. At the moment, the company can only obtain the financial balance sheets from the Ministry of Finance, from which it can compute the aggregate information by itself. It would be more efficient if the Ministry provided an e-government service for retrieving the required aggregate information automatically. All the presented use cases can benefit from using the techniques and data structures presented in this paper, as most of the desired information can be easily expressed as range aggregate queries.

Multidimensional Data Structures

Multidimensional Range Minimum Queries

The RMQ technique [START_REF] Bender | The LCA Problem revisited[END_REF] is very versatile and can be extended to multi-dimensional arrays. Each query asks for the minimum (maximum) value in a range [r [START_REF] Resteanu | Distributed and Parallel Computing in MADM Domain using the OPTCHOICE Software[END_REF][START_REF] Resteanu | Distributed and Parallel Computing in MADM Domain using the OPTCHOICE Software[END_REF],r(1,2)] x … x [r(d,1),r(d,2)] of a d-dimensional array a. W.l.o.g., we assume that the array contains n cells in every dimension (and, thus, n d cells overall). We will consider the most general case, in which the dimensions can be split into e≤d groups (group(j)=the group of dimension j). Each group g (1≤g≤e) contains a base dimension bd(g). The length of the interval of any query for every dimension j is always equal to (f(j) x the length of the query interval in dimension bd(group(j))), with f(j)≥1 (f(j)=1 for j=bd(group(j))). Since all the dimensions except the base dimensions are "linked" to the base dimension in their group, we will compute the following values: m(c [START_REF] Resteanu | Distributed and Parallel Computing in MADM Domain using the OPTCHOICE Software[END_REF] ,…,c(d),k(1),…,k(e)) = the minimum value in a range where r(j,1)=c(j) and the length of the query interval in dimension bd(g) is 2 k(g) . We will consider the sequences (k(1), …, k(e)) in increasing lexicographic order. When at least one k(*) value is larger than 0, we can use eq. ( 1), where q(u)=1, if (k(group(u))>0), and 0, otherwise.

)} 1 ) ( ( )) ( ( ),..., 1 ) 1 ( ( )) 1 ( ( )..., 1 ,...( 2 ) ( ) ( ) ( ) ( ( { min )) ( ),..., 1 ( ), ( ),..., 1 ( ( 1 )) ( ( 1 ), ( ) ( } 1 , 0 { )) ( ),..., 1 ( ( - ⋅ - ⋅ ≤ ≤ ⋅ ⋅ ⋅ + = - ≤ ≤ ≤ ∈ e k e bd q k bd q d u u f u s u q u c m e k k d c c m u group k d u u q u s d s s d (1)

A better way to compute m(c(1), …, c(d), k(1), …, k(e))

is to select just one value k(j)>0, set q(j)=1, replace q(bd(u))•(k(u)-1) by k(u) in eq. ( 1) (for u≠j), and instead of considering all the tuples (s(1), …, s(d)) in eq. ( 1), we can consider only the tuples (s(1), …, s(d)) with s(p)=0 or 1, only if group(p)=j (and 0, if group(p)≠j). When all the k(*) values are 0, we distinguish between two cases. If every group contains only one dimension or all the f(j) values of the non-base dimensions j are 1, then m(c(1), …, c(d))=a(c(1), …, c(d)). Otherwise, the problem is reduced to computing the minimum value in a range

[r(1,1),r(1,2)] x … x [r(d,1),r(d,2)],
where r(j,1)=r(j,2) if j is a base dimension in its group, or r(j,2)=r(j,1)+f(j)-1, if j is not a base dimension. We can handle this case as O(n d ) (d-e)-dimensional RMQ queries (one for every tuple of coordinates). The (d-e) dimensions are the remaining (d-e) non-base dimensions, grouped as follows. We remove the base dimension bd(g) from each group g and set as the new base dimension of g the nonbase dimension j' in g with the minimum value f(j') (if it exists). The new f(*) values of the non-base dimensions j'' from a group g will be f(j'')/f(j') (if this is not an integer, we place j'' in a separate group). In order to compute the minimum value in a given range

[r(1,1),r(1,2)]x … x[r(d,1),r(d,2)], we compute the values k'(j)=floor(log 2 (r(j,2)-r(j,1)+1))
and set k(g)=k'(bd(g)) for every group g. The answer is:

)} ( ),..., 1 ( )..., 1 )),...( 1 , ( 1 2 ) ( ) 2 , ( ( ) ( ) 1 , ( ( { min )
) ( ( } 1 , 0 { )) ( ),..., 1 ( ( e k k d u u r u f u r u s u r m u group k d s s d ≤ ≤ - + ⋅ - ⋅ + ∈ (2)
The preprocessing time complexity is at most

O(2 d •d•n d •log d (n)). A query is answered in O(2 d ) time.
It is obvious that the min function can immediately be replaced by max, obtaining identical results for range maximum queries.

Multidimensional Range Aggregate Queries

In this section we consider several alternatives for multidimensional range aggregate queries over a (dynamic) data cube, in which the aggregate function is invertible. The data cube has a fixed number d of dimensions and has m(j) entries in every dimension j (1≤j≤d); the entries are numbered from 1 to m(j). A cell of the data cube has coordinates (c [START_REF] Resteanu | Distributed and Parallel Computing in MADM Domain using the OPTCHOICE Software[END_REF] ,…,c(d 1) prefix subcubes [START_REF] Ho | Range Queries in OLAP Data Cubes[END_REF]. Thus, considering only prefix subcube queries is enough. In fact, we can compute a prefix "sum" cube, with which we can compute the "sum" (aggregate) of a prefix subcube in O(1) time. However, when an update occurs, we need to recompute the prefix "sum" cube (which takes O(n d ) time). Another solution is based on constructing a multidimensional binary indexed tree [START_REF] Fenwick | A New Data Structure for Cumulative Frequency Tables[END_REF], which supports updates and prefix subcube queries in O(log d (n)) time each (with only O(n d ) memory consumption). Yet another solution is based on using a multidimensional block partitioning [START_REF] Andreica | Efficient Data Structures for Online QoS-Constrained Data Transfer Scheduling[END_REF], which supports updates in O(n d/2 ) time and queries in O(1) time. In this section we propose a novel data structure, which supports both updates and queries in O(n d/4 ) time. The main idea is to divide the dimensions into two disjoint sets: the first set consists of q dimensions and the second set consists of the remaining d-q dimensions. The O(n) entries in every dimension j are split into O(n/k) blocks (where k is a function of n) of (approximately) k consecutive entries each. The blocks are numbered starting with 1, as they appear in increasing order (in increasing order of the entries they contain). Moreover, for every entry p in the j th dimension we know the block number of the block which contains the entry, blk(j,p). For every combination (x(1), x(2), …, x(q)), where x(j) is either an entry in dimension j or a block in dimension j (thus, there are O((n+n/k) q ) tuples overall), we maintain a multidimensional block partition BP(x(1), …, x(d)) of the remaining d-q dimensions.

When an update U(u, c(1), …, c(d)) occurs, we proceed as follows. We consider all the tuples (x(1), …, x(q)), where x(j) is either an entry x(j)=c(j), or a block x(j)=blk(j, c(j)) (1≤j≤q; we consider O(2 q )=O(1) tuples overall). For each such tuple, we will update its multidimensional block partition BP(x [START_REF] Resteanu | Distributed and Parallel Computing in MADM Domain using the OPTCHOICE Software[END_REF] ,…,x(d)). We consider all the tuples (y(q+1), …, y(d)), where y(j) is either an entry y(j)≥c(j) with blk(j, y(j))=blk(j, c(j)), or a block y(j)>blk(j, c(j)) (q+1≤j≤d; we consider O((k+n/k) d- q ) such tuples). We update the cell (y(q+1), …, y(d)) of BP(x(1), …, x(q)) by u (e.g. add u to it, multiply it by u, xor it with u). Thus, an update takes O((k+n/k) d-q ) time.

In order to find the answer to a query Q(b(1), …, b(d)), we consider all the tuples (x(1), …, x(q)), where x(j) is either an entry x(j)≤b(j) with blk(j, x(j))=blk(j, b(j)), or a block x(j)<blk(j, b(j)) (we consider O((k+n/k) q ) such tuples). For each such tuple (x(1), …, x(q)), we will compute Qagg(x(1), …, x(q), b(q+1), …, b(d))=the aggregate of the values of the following entries (y(q+1), …, y(d)) of BP(x(1), …, x(q)): y(j) is either an entry y(j)=b(j) or a block y(j)=blk(j, b(j)) (thus, we consider O(2 d-q )=O( 1

) tuples (y(q+1), …, y(d))

). The answer to the query is the aggregate of the values Qagg(x(1), …, x(q), b(q+1), …, b(d)), over all the considered tuples (x(1), …, x(q)). Thus, a query takes O((k+n/k) q ) time. We can choose the parameters k and q as we see fit, possibly according to the expected ratio between queries and updates. Note that when k=n 1/2 and q=d, we obtain O(1) update time and O(n d/2 ) query time; if q=0 we move to the other end of the spectrum, where an update takes O(n d/2 ) time and a query takes O(1) time. If both queries and updates are equally probable, we can choose q=d/2 (and k=n 1/2 ), which, for even d, leads to O(n d/4 ) query and update times.

Multidimensional Medians

We consider n points on the real line, each having an xcoordinate x(i) and a weight w(i) (1≤i≤n). We want to place (at most) K intervals of fixed length L, such that the total sum of distances from the given points to the intervals is minimum. The distance from a point x(i) to an interval [a,b] is defined as: 0, if (a≤x(i)≤b); min{|x(i)-a|, |x(i)-b|}, otherwise. The problem of determining point K-medians (with fixed length L=0) has been studied before [START_REF] Fleischer | Online Maintenance of k-Medians and k-Covers on a Line[END_REF][START_REF] Andreica | Minimum Dissatisfaction Personnel Scheduling[END_REF] and optimal O(n•K) algorithms have been given. We will extend those algorithms to the L>0 case. In order to facilitate the solution to this problem, we will add a new point at coordinate x(i)+L (with weight 0), for every initial point x(i) (1≤i≤n), thus obtaining n'=2•n points. We consider the points sorted in ascending order of their coordinates: x(1)≤x( 2

)≤…≤x(n').

For each point i we compute the value pleft(i)=the smallest index such that x(i)-x(pleft(i))≤L. These values can be computed in O(n') time (if the points are sorted), with a sliding window-type algorithm. In an optimal solution, the medians have the right endpoint positioned at the x-coordinate of some point. We will compute the following values: D min (i,j,0)=the minimum total sum of distances of the first i points after placing j intervals and the rightmost interval has its right endpoint at x(i); D min (i,j,1)=the minimum total sum of distances of the first i points after placing j intervals and the rightmost interval has its left endpoint at an x-coordinate less than or equal to x(i). We have the following equations: D min (0,j,0)= D min (0,j,1)=0 and D min (i,0,0)=D min (i,0,1)= +∞, for i>0. For i>0 and j>0, we have:

{ } ∑ - + = < ≤ - - ⋅ + - = 1 ) ( 1 min ) ( 0 min )) ( ) ( ( ) ( ) 1 , 1 , ( min ) 0 , , ( i pleft p q i pleft p q x L i x q w j p D j i D (3) ∑ + = ≤ ≤ - ⋅ + = i p q i p x(p)) (x(q) w(q) ) (p,j, {D ) (i,j, D 1 min 0 min 0 min 1 (4)
These equations are very similar to the equations used in [START_REF] Andreica | Minimum Dissatisfaction Personnel Scheduling[END_REF] and, thus, the O(n•K) solution presented there can be easily adapted. The main difference consists of the fact that, when computing D min (i,j,0), we only have in the appropriate sorted double-ended queue (deque) values corresponding to candidate positions p (0≤p<pleft(i)), instead of the whole range [0,i-1]; the value corresponding to a position i' is inserted into the deque only when we reach a position i, such that i'<pleft(i).

A well-known case of the interval K-median problem is where K=1. Like in the general case, we add the n extra points x(i)+L with zero weights, obtaining n'=2•n points. We will now slide the interval from left to right, placing its right endpoint at every point. While sliding the interval, we will maintain four values: wd left , wd right and w left , w right , representing the total weighted distance of the points on the left (right) side of the interval and the total weight of these points. We assign a type type(i) to every point i: type(i)=0 if it is one of the original points and type(i)=-j if point i was added as the extra point corresponding to point j (thus, x(i)=x(j)+L). Afterwards, we sort the points in increasing order of their coordinates and renumber them in this order (including the negative type(i) values). The pseudocode below presents the algorithm. We consider the points sorted in increasing order of their coordinates:

Interval-1-Median(): right=1; wd left =w left =0; w right =sum of the values w(i) (2≤i≤n) wdm=wdc=wd right =sum of the values (x(i)-x(1))•w(i) (2≤i≤n) for i=2 to n' do { wd right =wd right -w right •(x(i)-x(i-1)); w right =w right -w(i) wd left =wd left +w left •(x(i)-x(i-1)) if (type(i)<0) then { w left =w left +w(-type(i)) }

wdm=min{wdm, wd right +wd left } wdc=min{wdc, max{wd right , wd left }} // only for L=0 } wdm is the minimum sum of weighted distances, corresponding to the interval 1-median. In the multidimensional case we are given n d-dimensional weighted points, at coordinates (x(i, [START_REF] Resteanu | Distributed and Parallel Computing in MADM Domain using the OPTCHOICE Software[END_REF] ,…,x(i,d)) and with weights w(i) (1≤i≤n). We want to place 1 hyperrectangle with side lengths L(j), 1≤j≤d (L(j) is the side length in dimension j) such that the sum of the distances (L 1 ) from the points to the hyper-rectangle is minimized. The distance from a point i to a hyper-rectangle whose lower corner is at (xr [START_REF] Resteanu | Distributed and Parallel Computing in MADM Domain using the OPTCHOICE Software[END_REF] ,…,xr(d)) is:

∑ =           + > + - ⋅ < - ⋅ + ≤ ≤ d 1 j L(j) xr(j) j) x(i, if L(j))), (xr(j) j) (x(i, w(j) xr(j) j) x(i, if j)), x(i, (xr(j) w(j) L(j) xr(j) j) x(i, xr(j) if 0, ( 5
)
We can decompose the problem into d onedimensional problems. For each dimension j, we can find the coordinate xr(j) independently of the other dimensions, by solving a 1D interval 1-median problem, considering points i (1≤i≤n) with x-coordinates equal to x(i,j) and having weights w(i), and the length of the interval median is L(j).

A situation which arises often in multidimensional data analysis is that of computing the median of a subrange of the given multidimensional data (using the L 1 metric). Let's consider the 1D case first and see how we can support efficient range weighted median queries, i.e. Q(i,j)=find the sum of weighted distances from the median of the points x(r), with i≤r≤j, to all the points in the range. We will show that, with appropriate preprocessing (O(n)), we will be able to answer every such query in O(log(n)) time. We define wsum(i,p) (i≤p), the sum of the weights of all the points between i and p (inclusive), wdsumLR(i,p)=the sum of the weighted distances from every point j (i≤j≤p) to point p, and wdsumRL(i,p)=the sum of the weighted distances from every point j (i≤j≤p) to point i. With O(n) preprocessing, we can compute each of these values in O(1) time. We compute wpsum(i)=the sum of the weights from point 1 to point i (wpsum(0)=0 and wpsum(i≥1)=w(i)+wpsum(i-1)), wdpsum(i)=the sum of the values w(j)•x(j) (1≤j≤i):

wdpsum(0)=0 and wdpsum(i≥1)=w(i)•x(i)+wdpsum(i-1) (obviously,

all of these values can be computed in O(n) time overall). Then, wsum(i,p)=wpsum(p)-wpsum(i-1), wdsumLR(i,p)=wsum(i,p)•x(p) -(wdpsum(p)-wdpsum(i-1)), and wdsumRL(i,p) = (wdpsum(p)-wdpsum(i-1))wsum(i,p)•x(i).

With these values, we can binary search the largest value r between i and j, such that wsum(r+1, j)-wsum(i,r-1)>0 (or r=i if the condition is never met). The optimal location r opt of the median is either r or r+1 (if r+1 is inside the interval [i,j]). The cost of placing the median at position r opt is wdsumLR(i,ropt)+ wdsumRL(ropt,j).

This method can be extended to multiple dimensions, as follows. Let's assume that we have a data cube with d dimensions, having m(j) distinct coordinates in each dimension j (1≤j≤d) (thus, the cube contains m(1)•m(2)•…• m(d) data points). The p th value in the j th dimension has an assigned coordinate, x(j,p) (1≤j≤d; 1≤p≤m(j)). Each value Cube(c(1), …, c(d)) is the weight of a point with coordinates (x( 1,c(1)),…,x(d,c(d))). We want to answer efficiently queries of the following type: find the location of the median (under the L 1 metric) of all the points with the coordinates in the range [a [START_REF] Resteanu | Distributed and Parallel Computing in MADM Domain using the OPTCHOICE Software[END_REF]

, b(1)]x…x[a(d), b(d)] (a point (c(1), …, c(d)) is within the range if x(j,a(j))≤x(j,c(j))≤x(j,b(j)), for every 1≤j≤d).

This multidimensional range weighted median problem can be reduced to d 1D median problems. The j th such problem considers the range [a(j),b(j)] among m(j) distinct points with appropriately chosen weights w(j,i) (1≤i≤m(j)). Let's assume that the median of the j th problem is located at x(j,r(j)). The median of all the points in the d-dimensional range is located at (x(1,r(1)), …, x(d,r(d))). In order to answer range weighted median queries, we will first construct d data cubes with the same size as the initial Cube: DCube j (c(1), …, c(d))=x(j,c(j))•Cube(c(1), …, c(d)). For each such data cube j (1≤j≤d), we will construct a prefix sum data cube: With these data cubes, we can solve the range weighted median problem for every dimension j (1≤j≤d) as follows. For the j th dimension, the interval [a(j), b(j)] is the 1D interval of points on which we will focus; this interval has b(j)-a(j)+1 points, each point p (a(j)≤p≤b(j)) having the coordinate x(j,p) and a weight w(p)= RangeSum (Cube, [u(k), v(k)] (1≤k≤d)), where u(j)=v(j)=p and u(j')=a(j') and v(j')=b(j') (for 1≤j'≤d, j'≠j). We now have to show how to compute the values wsum(i,p), wdsumLR(i,p) and wdsumRL(i,p). wsum(i,p) is RangeSum(Cube, [u(k), v(k)] (1≤k≤d)), where u(j)=i, v(j)=p, and u(j')=a(j') and v(j')=b(j') (for 1≤j'≤d, j'≠j). Let's denote by wdsum(i,p)=RangeSum(DCube j , [u(k), v(k)] (1≤k≤d)), where u(j)=i, v(j)=p, and u(j')=a(j') and v(j')=b(j') (for 1≤j'≤d, j'≠j). We have: wdsumLR(i,p)= wsum(i,p)•x(j,p)-wdsum(i,p) and wdsumRL(i,p)= wdsum(i,p)-wsum(i,p)•x(j,i). We can compute every value wsum(i,p), wdsumLR(i,p) and wdsumRL(i,p) in O(2 d )=O [START_REF] Resteanu | Distributed and Parallel Computing in MADM Domain using the OPTCHOICE Software[END_REF] time. Thus, we can find r(j) in O(log(n)) time. The overall time complexity is

O(d•log(n))= O(log(n)) (since d is a constant).

Multidimensional Range Selection

We consider d sorted arrays w(1), …, w(d), of size n (w(i,j)≤w(i,j+1), 1≤i≤d, 1≤j≤n-1). We want to select the k th smallest weight in the set of points with coordinates (c( 1 d,n))]. The feasibility test for a candidate weight wt consists of computing the number of points p whose weight is at most wt. if p≥k, then wt≥wk; otherwise, wt<wk. Computing p is easy when op=max. For each i (1≤i≤d) we compute limit(i)=the largest index j (0≤j≤n) such that w(i,j)≤wt (using binary search). We consider w(i,0)=-∞. Then, p=limit(1 op -1 denotes the inverse operation of op (i.e. op -1 =for op=+, and op -1 =/ for op=*). The time complexity of the feasibility test for op=+ or * is O(n d-1 •log(n)). The overall complexity of the algorithm is obtained by multiplying the complexity of the feasibility test by log(WMAX). When the weights are integers, the algorithm finds the exact solution. In case of real numbers, it finds the k th smallest weight with any fixed arbitrary precision ε>0 (the binary search ends when the length of the search interval is smaller than ε). If O(n q ) storage is available (1≤q≤d/2), then we could use the following technique. We consider the first q dimensions and generate and store the (multi)set S 1 of all the O(n q ) "sums": (w(1,c( 1)) op w(2, c( 2)) op … op w(q, c(q))), with 1≤c(i)≤n, 1≤i≤q. Then, we sort all these "sums" and denote by S 1 (j) the j th smallest sum in S 1 . Afterwards, we generate sequentially each of the O(n d-q ) "sums" considering the dimensions q+1, …, d. For each such "sum" S=(w(q+1,c(q+1)) op w(q+2, c(q+2)) op … op w(d,c(d))), with 1≤c(i)≤n, q+1≤i≤d, we binary search the largest index j, such that S 1 op -1 S) (j may be 0), where wt is the candidate weight. The sum of all these indices j is the number p used by the feasibility test. The time complexity of this approach is O(n max{q,d-q} •log(n q )). ComputeP is a particular case of this more general solution, with q=1.

Let's denote by ww(k) the k th smallest weight among all the points (1≤k≤n d ). We now want to be able to compute efficiently an aggregate agg of the k smallest weights, where agg=op=+, * or max. When agg=max, the answer is ww(k), since ww(i)≤ww(i+1) (1≤i≤n d -1). For the other cases, we could compute all the weights ww(1), …, ww(k), but this would be too inefficient. We will use the same binary search algorithm as before, but we also compute the prefix "aggregate" array PS 1 , such that PS 1 (j)=(S 1 [START_REF] Resteanu | Distributed and Parallel Computing in MADM Domain using the OPTCHOICE Software[END_REF] agg S 1 [START_REF] Stoica | Fuzzy Sets and Their Applications[END_REF] agg … agg S 1 (j)). We have PS 1 (j)=PS 1 (j-1) agg S 1 (j) (for j≥1) and PS 1 (0)=0 for agg=+ and 1 for agg=*. We will compute a value pagg, initialized to 0 (for agg=+) or 1 (for agg=*). Within the feasibility test, when computing for each "sum" S of the remaining O(n d-q ) "sums" the largest index j such that (S 1 (j) agg S)≤wt, we increment p by j and we set pagg to (pagg agg PS 1 (j) agg mop(S, j)), where mop(u,v)=u•v for agg=+, and mop(u,v)=u v for agg=*. At the end of the binary search, we obtain the weight ww(k), together with the corresponding values p and pagg. It is possible that p≠k, when the weights are not distinct. The aggregate of the k smallest weights will be (pagg agg mop(ww(k), kp)).

Conclusions

In this paper we presented several novel multidimensional data structures and techniques for computing range aggregate queries efficiently in multidimensional databases (modeled as data cubes). The techniques reduce the time complexity dramatically, compared to the naive solutions. The proposed methods were thoroughly analyzed, mostly from a theoretical perspective. They have many applications in a wide range of domains, as was shown in Section 2.

Efficient range query techniques and data structures were developed in the context of OLAP data cubes [START_REF] Ho | Range Queries in OLAP Data Cubes[END_REF], multidimensional databases, data compression [START_REF] Fenwick | A New Data Structure for Cumulative Frequency Tables[END_REF], and even job and communication scheduling [START_REF] Andreica | Efficient Data Structures for Online QoS-Constrained Data Transfer Scheduling[END_REF]. Computing the median of a set of points has applications in personnel scheduling problems [START_REF] Andreica | Minimum Dissatisfaction Personnel Scheduling[END_REF]. As mentioned in the introduction, the presented techniques can be used in order to efficiently provide accurate input data to decision making optimization techniques [START_REF] Resteanu | Distributed and Parallel Computing in MADM Domain using the OPTCHOICE Software[END_REF][START_REF] Stoica | Fuzzy Sets and Their Applications[END_REF][START_REF] Ben-Ameur | Dynamic Programming Approach for Valuing Options in the GARCH Model[END_REF].

  )) (1≤c(j)≤m(j) ; 1≤j≤d) and has a value Cube(c[START_REF] Resteanu | Distributed and Parallel Computing in MADM Domain using the OPTCHOICE Software[END_REF], …, c(d)). A query Q(a(1), b(1), …, a(d), b(d)) consists of computing an aggregate (e.g. sum, product, xor) over the values of the cells (c(1), …, c(d)) whose coordinates are in the range: a(j)≤c(j)≤b(j) (1≤j≤d). An update U(u, c(1), …, c(d)) modifies the value of the cell (c(1), …, c(d)) by the value u (e.g. it increases/multiplies it by or sets it to u). We are interested in supporting both types of operations efficiently. The straight-forward solution is to update every cell in O(1) time (just change the value of Cube(c(1), …, c(d))) and compute the aggregate in O(n d ) time (n=max{m(1), …, m(d)}), by traversing every cell in the given range. For the approaches we present next, we will assume that the query has the form Q(b(1), …, b(d)) and asks for an aggregate of the values of all the cells in the range [1,b(1)]•…•[1,b(d)] (i.e. a prefix subcube of the data cube). The aggregate of any arbitrary range can be computed as a "sum" of the aggregates of O(2 d )=O(

  PSDCube j (b(1), …, b(d))=the sum of the values DCube j (c(1), …, c(d)), with 1≤c(j)≤b(j) (for every 1≤j≤d). As was shown in[START_REF] Ho | Range Queries in OLAP Data Cubes[END_REF], each prefix sum data cube can be computed in O(n d •d) time.Let's denote by RangeSum(X, [u(1), v(1)], …, [u(d), v(d)]) the sum of the values X(c(1), …, c(d)), with u(j)≤c(j)≤v(j) (1≤j≤d); if some v(j)=0, then the sum is 0. Using PSDCube j , we can compute RangeSum(Dcube j , [u[START_REF] Resteanu | Distributed and Parallel Computing in MADM Domain using the OPTCHOICE Software[END_REF], v(1)], …, [u(d), v(d)]) in O(2 d )=O(1) time. We will also compute a prefix sum data cube for the initial cube: PSCube(b(1), …, b(d)) is the sum of all the values Cube(c(1), …, c(d)), with 1≤c(j)≤b(j) (for every 1≤j≤d).

  ), …, c(d)) (1≤c(i)≤n, 1≤i≤d), where the weight of a point (c(1), …, c(d)) is (w(1,c(1)) op w(2, c(2)) op … w(d,c(d))), where op is + (addition), * (multiplication) or max. If O(n d ) storage is available, we could store all the point weights and then select the k th smallest weight in O(n d •(d+log(n d ))), by sorting the weights, or in O(n d •d) time, by using QuickSelect. Instead, we will binary search the k th smallest weight wk in the range [0,WMAX=(w(1,n) op … op w(

  )•…•limit(d). The time complexity of the feasibility test is O(d•log(n)).The other cases can be solved by the following recursive function:

	ComputeP(di, wt):
	if (di=1) then {
	binary search the largest index j (0≤j≤n), s.t. w(1,j)≤wt
	return j } else { // di≥2
	p=0; for j=1 to n do p=p+ComputeP(di-1, wt op -1 w(di, j))
	return p }


Cédric Augonnet
email: cedric.augonnet@inria.fr

Samuel Thibault
email: samuel.thibault@labri.fr

Raymond Namyst
email: ray-mond.namyst@labri.fr

Keywords: GPGPU, Multicore, Scheduling, Performance, StarPU

most famous example of a parallel computer mixing IBM Cell Broadband Engines and AMD opteron processors. Other architectures, featuring GPU accelerators, are expected to appear in the near future. To fully tap into the potential of these hybrid machines, pure offloading approaches, in which the main core of the application runs on regular processors and offloads specific parts on accelerators, are not sufficient. The real challenge is to build systems where the application would permanently spread across the entire machine, that is, where parallel tasks would be dynamically scheduled over the full set of available processing units.

To face this challenge, we propose a new runtime system capable of scheduling tasks over heterogeneous, accelerator-based machines. Our system features a software virtual shared memory that provides a weak consistency model. The system keeps track of data copies within accelerator embedded-memories and features a data-prefetching engine. Such facilities, together with a database of self-tuned per-task performance models, can be used to greatly improve the quality of scheduling policies in this context. We demonstrate the relevance of our approach by benchmarking various parallel numerical kernel implementations over our runtime system. We obtain significant speedups and a very high efficiency on various typical workloads over multicore machines equipped with multiple accelerators.

StarPU: un support exécutif capable d'ordonnancer des tâches sur des machines multicoeurs équippées d'accélérateurs

Résumé : Les machines multicoeurs équippées d'accélérateurs sont de plus en plus répandues. Le RoadRunner, une des machines parallèles actuellement les plus puissantes au monde est d'ailleurs constitué d'une combinaison de processeurs IBM Cell Broadband Engines et d'AMD Opteron. D'autres architectures à base d'accélérateurs (notamment de GPUs) devraient également faire leur apparition dans un futur proche. Afin de réellement tirer la quintessence de ces machines hybrides, il n'est plus suffisant de se contenter d'approches uniquement fondées sur le déport de calcul, où le coeur de l'application tourne sur des processeurs classiques alors que les parties véritablement coûteuses sont déportées sur des accélérateurs. Le véritable défi consiste à concevoir des systèmes où les applications sont réparties sur la totalité de la machine, c'est à dire des qu'il faut être capable d'ordonnancer des tâches parallèles sur la totalité des ressources de calcul disponibles.

Afin de faire face à ce défi, nous présentons un nouveau support exécutif capable d'ordonnancer des tâches sur des machines hétérogènes à base d'accélérateurs. Notre système comporte une mémoire virtuellement partagée qui assure la cohérence des données. Il garde la trace des copies de chacune des données dans les différentes mémoires embarquées sur les accélérateurs, et fourni des mécanismes tels que le pré-chargement de données. Conjointement à l'utilisation de modèles de coût auto-calibrés, ces fonctionnalités permettent d'améliorer substantiellement la qualité de l'ordonancement dans le contexte des architectures hybrides. Nous démontrons l'intérêt de notre approche en analysant les performances de notre support exécutif sur différents algorithmes numériques parallèles. Nous observons une réduction significative des temps de calcul, ainsi qu'une grande efficacité dans l'utilisation des différentes ressources de calcul pour différentes charges de travail typiques, notamment dans le cas de machines multicoeurs équippées de plusieurs accélérateurs.

Mots-clés : GPGPU, Multicoeur, Ordonnancement, Performance, StarPU.

Introduction

The High Performance Computing community has recently witnessed a major evolution of parallel architectures. The invasion of multicore chips has impacted almost all computer architectures, going from laptops to high-end parallel computers. While researchers are still having a hard time trying to bridge the gap between the theoretical performance of multicore machines and the sustained performance achieved by current software, they now have to face yet another architecture trend: the use of specific purpose processing units as side accelerators to speed up computation intensive applications.

The Cell Broadband Engine processor, for instance, features a general purpose core (PowerPC Processing Unit) surrounded by eight co-processors (Synergistic Processing Units), each featuring a specific vector instruction set and equipped with a private local storage. This processor is currently used in several multimedia hardware solutions -it most notably powers the Sony PlayStation 3 -but it is also a key component of the currently most powerful parallel machine: IBM RoadRunner. Although the future of the Cell processor in HPC is not clear, many computer architects predict that future processors will follow a similar design mixing general purpose cores and specialized co-processors.

Another popular practice is to use General Purpose Graphical Processing Units as co-processors to speed up data parallel computations. GPUs are massively parallel devices that can run SIMD programs very efficiently. Many research groups have already reported impressive performance boost thanks to the offloading of part of their application to a GPU device. Although the generalization to multi-GPU computation is currently limited by the bandwidth of the I/O buses connecting these devices to the host machine, some computer manufacturers actually plan to build HPC architectures featuring multiple GPUs racks. It is however expected that GPU-like co-processors and general purpose cores will be part of the same chip in the near future. The future Intel Larrabee processor, for instance, clearly belongs to that category of hybrid processors.

As a result, for an increasing part of the HPC community, the challenge has somehow moved from exploiting hierarchical multicore machines to exploiting heterogeneous multicore architectures. In the former case, many research efforts are merely dedicated to enhancing existing compilers, runtime systems or parallel libraries to better meet the requirements of new multicore architectures. Regarding heterogeneous machines, the gap with traditional multiprocessor architectures is so huge that simply enhancing existing solutions is not an option. Accelerator-based machines indeed introduce many hardware constraints -related to several memory limitations and to the execution model of accelerators -that make traditional multithreaded approaches impractical. Indeed, accelerators such as GPUs or Cell's SPUs have their own local storage which is not implicitly kept consistent with the main host memory. Thus, all data transfers between the main memory and the accelerators must be explicitly done by the software. Furthermore, because accelerators feature their own instruction set and a specific execution model, the code offloaded to these co-processors must be generated by a specific compiler and is completely different from the code compiled for the regular cores, which makes it impossible to migrate tasks between processors and accelerators during their execution.

As a result, exploiting this new generation of parallel computers requires to revise the complete software stack, from the languages down to the runtime systems. Our contribution is the design of a new runtime system, called StarPU, that meets the needs of heterogeneous, accelerator-based architectures. We explain how StarPU is capable of scheduling tasks efficiently over the whole range of processing units, thanks to the use of auto-tuned performance models and optimized data transfers between processing units. Our runtime system is typically intended to serve as a target for numerical kernel libraries and parallel language compilers. We use the StarPU implementation over NVIDIA GPU cards and Intel Nehalem processors to perform detailed performance studies of our scheduling techniques on several numerical kernels. The results of these studies provide insights into the impact of the accuracy of performance models and the granularity of the data partitioning.

Background

We first review the main types of approaches to the problem of programming accelerator-based machines.

New programming tools

To deal with accelerator programming, a wide spectrum of languages, techniques and tools that have been specifically designed, ranging from low-level development kits to new parallel languages.

Low-level development kits

Constructors usually provide low-level software development toolkits (SDK) to program their particular type of accelerators. Since such accelerator-specific APIs are usually the first available for the corresponding acceleration technologies, they are usually also the most widely used. CUDA, for instance, is the official SDK provided by NVIDIA for their GPU cards, while the AMD company provides Stream, a totally different SDK for their ATI cards. Concerning programming the Cell, LibSPE provided by IBM still remains the main SDK.

All these SDKs look similar at first sight, but they also export many hardware capabilities to the application, allowing to develop very efficient, carefully tuned code. However, their low-level API is hard to use and leads to nonportable programs. They also do not provide support for distributing work between regular cores and accelerators, they merely provide tools to drive computations on the latter.

The recently announced OpenCL standard provides a portable low-level interface, but it does not really help people in handling the complexity of such challenging hardware: such kind of low-level library does not really take any decision, and just obeys to the programmer's orders. In order to achieve the best performance a programmer would still have to know the intricacies of the underlying hardware, to give judicious orders.

Optimized Kernel Libraries

Because many HPC applications are making extensive use of common numerical kernels (BLAS for linear algebra [START_REF] Volkov | Benchmarking gpus to tune dense linear algebra[END_REF], FFT for Fourier transforms [START_REF] Naga | High performance discrete fourier transforms on graphics processors[END_REF]), it is no surprise that numerous research groups have designed implementations of numerical kernel libraries specific to GPU or Cell accelerators [START_REF] Owens | A survey of general-purpose computation on graphics hardware[END_REF]. These highly optimized libraries are typically built using the aforementioned software development toolkits. For instance, the CUBLAS library provides optimized BLAS routines for NVIDIA GPU cards, and is implemented in CUDA.

Most existing libraries are only capable of exploiting a single accelerator at a time, though. Some libraries such as MAGMA consider not only multiaccelerator setups, but they are also currently getting support for using CPUs and accelerators simultaneously [START_REF] Tomov | Towards Dense Linear Algebra for Hybrid GPU Accelerated Manycore Systems[END_REF]. However, their implementation is based on a static distribution of tasks and data, which would not be suitable for complex systems composed of multiple heterogeneous accelerators. Furthermore, because existing implementation do not feature a true task scheduler, a parallel numerical kernel can not be run in parallel with another instance of a kernel without resulting to interferences. This problem prevents efficient parallel code composability.

Programming Languages

Several new specific languages have been designed to write portable programs in a high-level fashion that should take advantage of accelerator-based machines.

Many new languages are based on a streaming model. RapidMind and BrookGPU provide convenient ways to explicitely manipulate data streams. Other languages such as Sequoia or the MapReduce framework let the programmer express computations following some structure that exposes a sufficient amount of parallelism.

There are also numerous propositions to extend existing languages with some OpenMP-like constructs. StarSs and CellGen both propose a source-to-source approach to perform target-specific code optimization (eg. vectorization or tiling) and propose pragmas to delimit tasks and describe data accesses. Similarly, PGI accelerator and HMPP also take care of transforming annotated C/Fortran codes into efficient kernels for accelerators.

While those efforts attempt to save programmers from low-level architecturespecific interfaces thanks to high-level and portable constructs, most of them concentrate on generating efficient code. For instance, HMPP [START_REF] Dolbeau | HMPP: A hybrid multi-core parallel programming environment[END_REF] relies on a minimalist runtime system that simply uses a coarse-grain lock to protect each processing units. Some languages such as Sequoia even map computations in a static fashion [START_REF] Fatahalian | Sequoia: Programming the memory hierarchy[END_REF]; this may not be realistic when it comes to irregular applications and/or complex heterogeneous platforms. Runtime systems permit to do dynamically what cannot be done statically anymore.

Towards new runtime systems

The numerous research efforts presented in the previous sections cover a wide range of functionalities, but still suffer from relying on basic runtime systems. We claim that implementations of numerical kernel libraries or parallel languages could achieve better performance over complex heterogeneous architectures if they would use a runtime system capable of dynamically scheduling tasks on the most appropriate processing unit.

Designing such a runtime systems raises a number of challenges related both to the scheduling itself and to the management of data movements between the host memory and the accelerators. Scheduling. To allow the scheduler to dynamically assign tasks to any processing unit, tasks should have several implementations, one per type of processing unit that is capable of executing it. These implementations will obviously exhibit very different performance, but the gap between the implementation for an Intel Nehalem processor and the implementation for an NVIDIA Quadro GPU can also vary a lot depending on the kind of computation itself. Thus, the scheduler should be able to use performance prediction information as an input for its scheduling policy, so as to compute the best possible task assignment. Moreover, we think that the ability to easily tune and optimize scheduling policies incrementally is fundamental on this new generation of architecture, as we will show in the following sections. Data Management. Data management was already a delicate issue on multicore machines, due to the Non-Uniform Memory Accesses introduced by the topology formed by the memory banks. On archictectures equipped with accelerators, the situation is even worse, especially if the accelerators are connected via a slow I/O bus. Numerous data movements between the main memory and the embedded memory of an accelerator can really hurt performance. Thus, it is crucial that the runtime system keeps track of data copies linked to each processing unit so as to avoid unnecessary data movement whenever possible. Moreover, the runtime can go much further and estimate the data transfer times to better compare different scheduling alternatives, or even provide a powerful data management system able to optimize transfers of complex data structures (cyclic column distributions, etc.) and perform data prefetching to overlap part of the transfer cost. The next section presents our answers to these challenges.

The StarPU runtime system

StarPU is a runtime system for scheduling a graph of tasks onto a heterogeneous set of processing units, intended to provide the powerful features detailed in Section 2.2 in a portable but still efficient way. It is meant to be used as a backend for e.g. parallel language compilation environments and High-Performance libraries, here called the application. The two basic principles of StarPU is firstly that tasks can have several implementations, for some or each of the various heterogeneous processing units available in the machine, and secondly that necessary data piece transfers to these processing units are handled transparently by StarPU.

Figure 1 gives an overview of the journey of tasks within StarPU, thus providing a general idea of how StarPU modules are related. The application first submits tasks [START_REF] Augonnet | Automatic Calibration of Performance Models on Heterogeneous Multicore Architectures[END_REF]. When a task gets ready [START_REF] Augonnet | Exploiting the Cell/BE architecture with the StarPU unified runtime system[END_REF], it is dispatched to one of the device drivers by the scheduler (3). The DSM ensures the availability of all pieces of data (4). The driver then offloads the proper implementation for the task [START_REF] Ayguadé | An Extension of the StarSs Programming Model for Platforms with Multiple GPUs[END_REF]. When the task completes, tasks depending on it are released and an application callback for the task is executed [START_REF] Bellens | Exploiting Locality on the Cell/B.E. through Bypassing[END_REF]. Before further sections provide more details, the following subsections give the rationale for the DSM and the scheduler modules, as well as a complete StarPU programming example.

A virtually shared memory system for heterogeneous chips

Because StarPU schedules tasks at runtime, data transfers have to be done automatically and "just-in-time" between processing units, relieving the application programmer from explicit data transfers. Moreover, to avoid unnecessary transfers, StarPU keeps data where it was last needed, even if was modified there, and it allows multiple copies of the same data to reside at the same time on several processing units as long as it is not modified. To summarize, the data management part of StarPU, detailed in Section 4, actually implements a Software DSM (Distributed Shared Memory), with relaxed consistency and data replication capability thanks to an MSI (Modified/Shared/Invalid) protocol. The application just has to register the different pieces of data by giving their addresses and sizes in the main memory. Data can thus be dynamically as well as statically allocated.

StarPU also provides an additional high-level abstraction to easily handle partitioned data for block-and tile-based parallelism: filters. These can be used to logically divide data into smaller blocks, according to the application algorithm needs (e.g. MAGMA [START_REF] Tomov | Towards Dense Linear Algebra for Hybrid GPU Accelerated Manycore Systems[END_REF]). For instance, Figure 2 illustrates how a matrix can be first divided into 3 * 3 blocks, and then one of these blocks can be further divided into 2 * 2 sub-blocks, resulting in a tree of sub-data. Multiple filters are provided and can be extended to match various splitting strategies (bloc, cyclic, etc.) and various native memory access patterns (dense BLAS matrices, sparse CSR matrices, etc.). Such filters can be applied dynamically at runtime, to stick to the algorithm behavior, in particular application featuring adaptive refinement methods, and allow the DSM to transfer to the processing unit only the exact needed data, or to find it inside an already-transferred data containing it.

A tasking model enabling heterogeneous scheduling

StarPU tasks can be executed by as many processing units as possible, those which the programmer has provided an implementation for, including accelerators which do not have access to the main memory. Therefore, the programmer has to also explicit all the input and output data of tasks by using the handles returned by the DSM, so that the latter can automatically fetch data as needed right before execution.

The submission of tasks is asynchronous and termination is signaled through a callback. This lets the application submit several tasks, including tasks which depends on others. To express dependencies and thus actual graphs of tasks, integers called tags are used to characterize the termination of some task. Figure 3 illustrates how the dependencies between 4 tasks can be expressed. The advantage of decoupling tags from the tasks themselves is that a tag can depend on another tag which has not been associated with a task yet, or on another tag whose task is already finished. This lets the application avoid having to submit all the tasks at once (because it would completely fill the memory or just does not make sense, in the case of flow applications). A tag can also be explicitly waited for by the main application to synchronize with the execution of StarPU tasks.

Developing on top of StarPU: a case study

Figure 6 gives a global idea of how using StarPU looks like. It shows the implementation of a blocked matrix multiplication taking advantage of all CPUs and GPUs by leveraging both the existing CUBLAS and BLAS SGEMM kernels at the same time.

On the left part, the CUBLAS and BLAS SGEMM operations are first wrapped into generic functions to harmonize the calling convention, and the cl structure groups them into just one codelet which will take 3 matrix parameters and which is able to run on either a CPU or a GPU.

The right part shows how to submit tasks executing that codelet to perform the computation block per block. The A, B, and C matrices are first registered, then partitioned into blocks as represented on the bottom left figure: N rows and M columns. A 2-dimension loop then creates N * M tasks executing the SGEMM codelet and assigns to them the proper sub-parts of the matrices. After waiting i,j i get_sub_data(C_handle, 2, i, j) for the termination of all tasks, the C matrix can be unpartitioned and synchronized, to make sure parts computed on GPUs are written back to the main memory.

j filter_col (M=4) filter_row (N = 3)
It is worth noting that this fast_SGEMM function, instead of waiting for task completion itself, could also just return a tag depending on the termination of all tasks, thus allowing the application to start other kinds of computations concurrently, and StarPU to schedule them as efficiently as possible along the SGEMM kernels. Also, such kind of code could be automatically generated for a range of BLAS operations.

Managing Data

We now detail the internal principles of the data management module of StarPU, in particular the DSM mechanisms being used.

Data replication and memory consistency

When the application registers some data to the DSM, StarPU allocates in the main memory an array to record the MSI states of that data on the different memory nodes, as illustrated by Figure 7(b). This permits to replicate the data on multiple memory nodes (multiple-reader Shared access) or to keep a modified version on one memory node (single writer Modified access) while keeping readwrite consistency thanks to the MSI protocol. Here, for instance, the memory node number 2, a GPU, requests both Read and Write access, and it does not have a copy of the data (Invalid state). The data is transferred for the Read access, and other copies are invalidated due to the requested Write access.

By default, modified data is not immediately written back to the main memory. This lazy approach permits StarPU to save bus bandwidth when data is re-used by several tasks. If however the programmer knows that the data will not be re-used he can request for a write-through strategy. When enough memory is not available on some processing unit for the data required by a task, StarPU performs a memory reclaiming pass to flush some of the useless data out.

As a result, even for a complex application with a graph of tasks dealing with various data whose total size can not fit into the local memory of processing units and for which having to issue data transfers by hand would be very difficult, StarPU lets the programmer just express the blocked or tiled computation tasks and data dependencies, and StarPU will transparently optimize the usage of the limited local memory.

Asynchronous data management and data requests

Data transfers can be very long due to the main bus typically being a bottleneck. Nowadays acceleration cards support asynchronous data transfer, i.e. these transfers can be overlapped with computations. StarPU's data management is also completely asynchronous by associating each memory node with a queue of pending data requests, possibly queued by various parts of StarPU, another memory node requesting data for some task for instance. It also makes dynamic prefetching of data quite natural: as soon as a task is scheduled to be run by some processing unit, the data transfer order can be queued so that the execution of the task can hopefully happen as soon as the processing unit has finished previous tasks, thanks to data being transferred in parallel. More generally the scheduling and the data management library of StarPU can collaborate to optimize the whole execution.

Scheduling tasks and data transfers

Because a generic scheduler can not achieve the best performance for all kinds of application with all hardware, StarPU aims at hosting a testbed to implement and choose between various powerful scheduling strategies in a portable way. This section details how schedulers can be implemented and describes a few scheduler examples already provided along StarPU.

Application

A portable scheduling framework

The role of the StarPU scheduler is to dispatch tasks onto the different processing units. From the point of view of the application, tasks are just pushed to the scheduler, while from the point of view of the drivers of the processing units, tasks are just poped from the scheduler, i.e. the simplest implementation is a single FIFO queue as is shown on Figure 8(a). A more complex scheduler, like the one illustrated in Figure 8(b) thus can boil down to organizing a set of queues (which can be FIFOs, stacks, deques, etc.) and implement the above-mentioned push and pop operations. All distribution, load balancing, etc. decisions are matters internal to the scheduler. This approach permits to design various schedulers independently of the application and the drivers being used, and to choose which one should be used at runtime. Table 1 gives a list of the different portable policies that have been implemented in StarPU and which will be described in the following subsections. More details can be found in a previous paper [START_REF] Augonnet | StarPU: A Unified Platform for Task Scheduling on Heterogeneous Multicore Architectures[END_REF].

Greedy Strategies

The simplest strategy consists in having all processing units share a single FIFO as illustrated in Figure 8(a). A refinement of this is to make the push method of the FIFO take programmer-provided priorities of tasks into account: more prioritized tasks will be queue closer to the output of the FIFO, resulting to the greedy StarPU scheduler, as opposed to the no-prio scheduler. Such strategy however offers poor scalability w.r.t. the number of accelerators due to locking contention, especially in case the tasks are fine grained. This is commonly solved by decentralizing task queues as shown in Figure 8(b). A Cilk-like [START_REF] Frigo | The implementation of the cilk-5 multithreaded language[END_REF] work-stealing strategy was implemented in the pop method of the scheduler to steal tasks from other queues if the local queue is empty, resulting to the ws scheduler.

These schedulers are very portable in that they work with generic tasks. However, as discussed in details later in Section 6, the completion time of tasks can very often be predicted by automatic cost models, at least with rough precision. Greedy strategies can get huge benefit from this kind of information.

Cost model-guided scheduling strategies

Figure 9 shows an typical example of the distribution of tasks by the greedy or ws scheduler onto a dualcore machine equipped with two GPUs with differing performance. The heterogeneity of performance between different processing units leads to a load imbalance, which could have been avoided if the scheduler knew that task completion takes so much time on CPUs.

A common approach to tackle it is to consider the computational power of each unit (e.g. sustained speed in GFlop/s). The w-rand StarPU scheduler takes it into account by randomly dispatching tasks onto the processing units with a probability proportional to the respective computational power. This avoids the situation of Figure 9 and brings much better load balancing at the end of the application completion. It however does not take into account that some tasks may e.g. perform really well on a GPU and actually tend to reach the theoretical computational power, while others do not perform so well, and should rather be scheduled on the CPU, leaving the GPUs to the former. More refined StarPU schedulers thus rather use models of the execution time of the different tasks instead of just the processing units.

Taking task duration into account

The HEFT scheduling algorithm (Heterogeneous Earliest Finish Time [START_REF] Topcuoglu | Performance-effective and lowcomplexity task scheduling for heterogeneous computing[END_REF]), implemented in the heft-tm StarPU policy, is illustrated in Figure 10. It makes use of performance prediction to keep track of the expected dates Avail (P i ) at which each processing unit will become available (after all the tasks already assigned to it complete). A new task T is then assigned to the processing unit P i that minimizes the new termination time with respect to the expected duration Est Pi (T ) of the task on the unit i.e.

min Pi Avail (P i ) + Est Pi (T )

Taking data transfers into account

The time to transfer data to e.g. accelerators is actually far from negligible compared to task execution times, and can sometimes even become a bottleneck. We therefore extended the heft-tm policy into the heft-tmdp policy which takes data locality into account thanks to the tight collaboration between the scheduler and the data management library presented in Section 4. In addition to the computation time, the scheduler computes a penalty based on the times T j→i (d) required to move each data d from P j (where a valid copy of d resides) to P i . Such penalty of course reduces to 0 if the target unit already holds the data, i.e. j = i. The resulting minimization is min Pi Avail (P i ) + Est

Prefetching data

As explained in Section 4.2, hardware often allows data transfers to overlap with the actual computations, thus maximizing the time spent computing. This is exploited by StarPU by having the scheduler request from the DSM data transfers for a task as soon as the placement of that task gets decided. The DSM then queues the required data transfer orders to the various device drivers and the scheduler can reduce its estimation of the transfer time. This can be particularly important when using multiple accelerators which do not support direct transfer, in which case the data has to be written back to memory first, hence at least doubling the transfer time. The heft-tm and heft-tmdp policies have thus been extended into the heft-tm-pr and heft-tmdp-pr policies.

Predicting task completion time

We described how task scheduling can benefit from models that predict the performance of task execution and data transfer. Diamos et al. mention [START_REF] Diamos | Harmony: an execution model and runtime for heterogeneous many core systems[END_REF] that their Harmony runtime could also benefit from such models. It is worth noting that in the case of StarPU, such models do not necessarily need to be extremely accurate, since scheduling is done at runtime and can compensate not-so-optimal decisions according to what actually happened. We now detail how StarPU automatically builds these models.

Task performance models

More or less automatic ways exist to build task performance models. An extreme approach is to have the programmer explicitly construct them. This however requires extra experimentation, as well as a lot of knowledge of both the algorithm and the underlying machine. Such work also needs to be repeated for each platform, due to details as simple as cache size, so that this is not a realistic approach in general.

For very common kernels and libraries (BLAS, typically), an important literature includes performance models for the main kernels. Using them is possible with some extra programming efforts to tune their parameters according to the architecture. These models are however sometimes pretty simplistic (e.g. limitations in the cache size), and offline precalibration often does not capture real-world use, where contention and cache sharing can severely impact performance.

In order to get performance models which actually match what can be observed with real applications, we use history-based prediction models. When a task is executed, StarPU measures its duration, and to schedule a task, it can consult the history to take into account the performance obtained during previous executions.

This assumes some regularity hypothesis: execution time should be independent from the content of data itself, thanks to the flow control being static for instance. This is accurate for the kernels mentioned above. When provided with their literature performance models, StarPU can thus tune the latter during the execution by regression analysis, according to the underlying architecture, but also according to the general contention impact of the application.

In addition to that, StarPU can use history-based prediction for kernels without literature performance model, by assuming that for a given application they are mostly always called with the same parameters, the typical task granularity for instance. Given task parameters characteristics (the size of the data, typically), StarPU computes a hash characterizing this task's complexity, and that can thus be used as an index in the history tables from which an average of the previous execution times can be obtained. Figure 11 illustrates the resulting feedback loop between these tables, the Scheduler, and the actual execution measurement The history is saved to a file to be re-used for further executions. More details can be found in a previous paper [START_REF] Augonnet | Automatic Calibration of Performance Models on Heterogeneous Multicore Architectures[END_REF].

The key advantages of this approach is not only being efficient, but also being transparent for the programmer (he already has to specify the data size for transfers anyway), and it can thus be applied to any kernel which has regular execution times. When unsure about the latter property, StarPU also computes the standard deviation of the measurements and can even provide a histogram of the distribution of execution times, so that the variability of performance can be checked. StarPU can also generate a detailed Gantt chart of the actual execution, which helps a lot in understanding performance, particularly on complex heterogeneous architectures for which interference can happen between several different parts of the application for instance..

Data transfer models

Compared to the internal memory bandwidth of accelerators (several hundreds of GB/s), the speed of the PCI bus (half a dozen GB/s) is a real concern, and data transfer time can become as important as the computation time itself and thus be critical for scheduling decisions. To estimate the time to transfer some data, StarPU simply uses an estimation of the latency and the asymptotic bandwidth between these two memory nodes, which can be either measured through offline experiments, or, better, inferred from online measurements through regression analysis.

This approach catches various architecture artifacts, like heterogeneity and topology of PCI buses (8x vs. 16x and NUMA effects), as well as asymmetry of transfer paradigms (upload and download speeds are usually different on CUDA devices).

Porting StarPU to various hardware

As a reminder of Figure 1, when a task is submitted by the application, it first has to wait for its dependencies, then gets dispatched by the scheduler to some driver, the DSM ensures that all the data are available and the driver eventually offloads the codelet. When the latter completes, tasks depending on it are released and a callback is executed.

Porting StarPU to a new kind of device actually boils down to implementing functions transferring data between the device and the main memory, and executing tasks on the accelerator. To achieve the latter, drivers run an execution loop. For instance, Figure 12 summarizes the loop of the driver for a CPU core. It continuously fetches tasks from the scheduler (pop_task), requests the DSM to make sure all data are available (fetch_task_input, which possibly calls transfer functions of other drivers to fetch data back to main memory), and actually calls the codelet. After codelet termination, data can be released (push_task_output) and tasks depending on this task can be released and the callback called (handle_task_termination). The driver work is thus completely decoupled from the scheduling decision and DSM operation.

Heterogeneous multi-GPU

The implementation of the driver for an NVIDIA GPU is relatively straightforward: data transfer functions use the API from NVIDIA and the driver can directly call the codelet's cuda_func() function1 which makes standard CUDA or CUBLAS kernel calls. Some shortcomings of CUDA raise issues when handling multiple GPUs, but the design of StarPU overcomes them quite naturally.

Due to the CUDA specification and implementation, switching between the CUDA context of different GPUs is very costly. StarPU's CUDA driver thus runs one execution loop in a separate kernel thread for each GPU device (this is hence where the codelet's cuda_func() function is called). Because the CPU consumption of these threads is far from negligible, StarPU devotes CPU contexts (cores or hyper-threaded contexts) for them.

The CUDA API also does not permit to issue direct GPU-GPU memory transfers, they have to be achieved in two steps through the main memory. Thanks to its chained asynchronous data transfer described in Section 4.2, StarPU overcomes this by issuing the two transfer orders to the two driver instances.

The Cell B/E

The Cell B/E processor is a heterogeneous multicore chip composed of a main hyper-threaded core, named PPU (Power Processing Unit), and 8 coprocessors called SPUs (Synergistic Processing Units).

The StarPU port to the Cell [START_REF] Augonnet | Exploiting the Cell/BE architecture with the StarPU unified runtime system[END_REF] actually relies on the Cell-RTL (Cell Run Time Library) [START_REF] Nijhuis | Mapping and synchronizing streaming applications on cell processors[END_REF] to perform task offloading and manage data transfers between the main memory and local stores (LS) on the SPUs. Figure 13 depicts how the Cell-RTL works: when a task for some SPU is submitted on the PPU, the latter sends a message to the former (1); an automaton running on each SPU reads the message, fetches data into the LS (2), executes the corresponding task (3), commits the output data back to the main memory (4), and sends a signal to the PPU (5); when the latter detects this signal (6), a termination callback is executed [START_REF] Bellens | Cellss: Scheduling techniques to better exploit memory hierarchy[END_REF]. To greatly reduce the synchronization overhead, the Cell-RTL can submit chains of tasks to SPUs.

To fully exploit the possibilities of the Cell-RTL, the StarPU Cell B/E driver loop automatically builds chains of tasks before submitting them. Thanks to SPUs being almost full-fledged cores, the flexible Cell-RTL API permits to devote only one thread to submitting tasks to the different SPUs.

In a way, StarPU leverages the Cell-RTL by adding scheduling facilities and providing with the high-level data management and task dependencies enforcement, permitting efficient task chaining. StarPU could leverage other backends like IBM's ALF [START_REF] Crawford | Accelerating computing with the cell broadband engine processor[END_REF] or CellSs' runtime [START_REF] Bellens | Exploiting Locality on the Cell/B.E. through Bypassing[END_REF].

Evaluation

To validate our approach, we show that StarPU efficiently uses all the available processing units with a low overhead, and we analyze the performance of some numerical algorithms under different scheduling strategies. Eventually, we show that selecting the best granularity is a challenging issue on accelerator-based platforms.

Experimental testbed

To evaluate the performance of StarPU, we have used two platforms both running Linux 2.6 and CUDA 2.3:

A homogeneous multi-GPU platform composed of two X5550 (Nehalem) hyper-threaded quad-core processors running at 2.67 GHz with 48 GB of memory divided in two NUMA nodes. It is equipped with 3 NVIDIA Quadro FX5800 with 4 GB of memory.

A heterogeneous multi-GPU platform composed of a E5410 Xeon quad-core processor running at 2.33 GHz with 4 GB of memory. It is equipped with an NVIDIA Quadro FX5800 with 4 GB of memory and a: NVIDIA Quadro FX4600 with 768 MB of memory; the second GPU does not support asynchronous data transfers.

Distribution of the computation among heterogeneous platforms

Even with its simplest greedy, scheduler, StarPU simultaneously takes advantage of all available resources according to their respective computational power: Figure 14 (resp. 15) shows the proportion of tasks attributed to the different processing units in a multi-GPU system (resp. heterogeneous multi-GPU system). On the left-hand side, we have a blocked matrix-product, and on the right-hand side, we have a band-pass filter implemented using FFTW and CUFFT. In both benchmarks, the GPUs become relatively more efficient than the CPUs and thus get attributed more tasks when the granularity increases. Figure 15 illustrates that StarPU is able to distribute task onto different models of GPUs with respect to their respective speed: a Quadro FX5800 is given more tasks than a Quadro FX4600 which is much less powerful.

Overhead evaluation

It is important to make sure that StarPU only has a minimal overhead compared to its potential gains. Table 2 gives the typical overhead measured during the execution of chains of one million empty tasks, and by repeating this experiment 128 times. We distinguish the overhead of task submission, data management (for a single piece of data), and the overhead of the task execution itself. The second line gives the performance when this experiment is performed on 8 CPUs which are all associated with only one single queue. Even in this extreme case the overhead remains low, which shows that StarPU properly handles highly contended situations with a typical overhead smaller than a few micro-seconds.

In comparison, Volkov and Demmel report that CUDA functions typically incur a 3 to 7 µs overhead for asynchronous launch, and 11 µs for synchronous calls [START_REF] Volkov | Benchmarking gpus to tune dense linear algebra[END_REF]. Since data transfers are also taking tens of micro-seconds, we conclude that the overhead introduced by StarPU is acceptable with respect to the performance improvement resulting from scheduling policies.

A case study: LU and Cholesky decompositions

Implementing LU and Cholesky blocked decompositions for multicore machines equipped with multiple GPUs is straightforward with StarPU, even for problems that are larger than the memory available on the accelerators.

Prioritizing tasks

It is however crucial to maintain enough parallelism to feed such a massively parallel machine during the entire application. The LU and Cholesky decomposition algorithms are known to suffer from a lack of parallelism, that is usually avoided by performing critical tasks as soon as possible by the means of task priorities.

Figure 16 shows the evolution of the number of ready tasks during the execution of a Cholesky decomposition running on the homogeneous multi-GPU machine with our best scheduling strategy, heft-tmdp-pr. While every tasks has the same priority on the top curve, we put a maximum priority for the critical tasks in the second case. Priority-aware scheduling here prevents substantial loss of parallelism, so that we observe a typical 15 % speed improvement.

Modeling task performance

Figure 17 shows the speed of LU decomposition according to the problem size, and for various scheduling policies. The second bottom-most line gives the performance obtained by the heft-tm strategy with the support of the automatically tuned history-based performance models presented in Section 6.1. While no application code modification was required compared to the greedy policy, the heft-tm strategy significantly outperforms the simple greedy scheduling policy, especially for small problems which tend to have little parallelism. In these cases, it is indeed critical that the most efficient processing units (ie. GPUs) get the most critical tasks.

In a previous paper [START_REF] Augonnet | StarPU: A Unified Platform for Task Scheduling on Heterogeneous Multicore Architectures[END_REF], we have shown that in the case of a quad-core machine equipped with a single GPU, this strategy obtains a super-linear efficiency: the processing power of the hybrid system is slightly higher than the sum of the powers of a GPU and of the CPUs, taken individually. This is possible because the LU decomposition is composed of different types of tasks, with different relative speedups between their CPU and their GPU implementations. A matrix product (gemm) could be accelerated 20 times by a GPU, while the resolution of a triangular system (trsm) may only be accelerated 10 times by a GPU. With the heft-tm policy, the processing units therefore tend to receive the types of tasks that they most efficiently handle, and they process less tasks for which they are relatively inefficient. StarPU thus actually takes advantage of the heterogeneous nature of accelerator-based platforms.

Scheduling data transfers

Taking data transfers into account while scheduling tasks is critical. This is especially important for multi-accelerator platforms where memory buses are major bottlenecks.

The heft-tm-pr policy attempts to mask the cost of memory transfers by overlapping them with computations thanks to the prefetch mechanism described in Section 5.3.3. Figure 17 shows a typical 20 % speed improvement over the heft-tm strategy which does not prefetch data.

However, data prefetching does not reduce the total amount of memory transfers. In Section 5.3, we described how the heft-tmdp-pr strategy extends heft-tm-pr by penalizing non-local data accesses. The average activity on the bus is shown at the top of Figure 18: penalizing remote data accesses, and therefore favoring data locality has a direct impact on the total amount of data transfers which drops almost by half. Likewise, the bottom of Figure 18 shows that the processing units tend to work more locally. On Figure 17, this translates into another 15 % reduction of the execution time.

How about granularity?

In this section, we show that selecting the best granularity is not trivial, especially in the context of heterogeneous platforms, and we give some insights about how StarPU could handle this challenge. Figure 19 shows that depending on the problem size, the granularity for which the best result is obtained can vary. Separately precalibrating the various computational kernels individually for the different architectures is thus not sufficient. Libraries relying on precalibration to detect optimal block sizes [START_REF] Whaley | Automatically Tuned Linear Algebra Software[END_REF][START_REF] Li | A Note on Auto-tuning GEMM for GPUs[END_REF] should therefore take into account the factors that influence the value of the optimal granularity.

Different factors may indeed affect the choice of the grain size in the context of heterogeneous machines equipped with one of more accelerators. The relative efficiency of a task on various architectures depends on the data size: GPUs are for instance large vector systems that are only efficiently used when all the vector processors are used simultaneously, which is not possible with tiny data input. On the one hand, a small grain size means that there is a lot of small tasks, which potentially leads to a high overhead. On the other hand, too big a grain size causes processing units to become idle due to lack of parallelism.

Selecting the most appropriate granularity is therefore a matter of trade-offs when dealing with heterogeneous platforms equipped with accelerators.

Dealing with multiple grain sizes

As a result, while there can be an optimal grain size for a specific problem on a homogeneous machine, it may be that the best solution in a heterogeneous context is a mixture of multiple granularities.

A first simple approach is to explicitely decompose an application into multiple phases with different granularity, according to the evolution of the behavior of the computation. In the case of LU decomposition, we observe that there is a lack of parallelism at the end of the algorithm. Intuitively, we would like to offload tasks with a big granularity at the beginning of the algorithm to get goot performance on GPUs, and we should reduce the grain size when the amount of parallelism becomes too low to avoid load imbalance. Table 3 gives the speed of an LU decomposition on a 1 GB problem with the heft-tmdp-pr strategy. The first two columns gives the speed obtained with either (512×512) or (1024×1024) block sizes. The last column gives the performance obtained by starting with blocks of (1024×1024), and by reducing the grain size to (512×512) at the middle of the algorithm execution. The performance improvements by this approach which dynamically adapts the granularity are promising: the hybrid solution gets better performance than the best performance obtained with a single grain size (see the 1 GB curve of Figure 19 which does not reach 500 GFlop/s).

Besides, it is possible to use multiple granularities at the same time, either explicitly or implicitly. Programmers could submit tasks with multiple grain size by hand, but this requires significant programming efforts. This approach is applicable for applications which explicitly select the tasks that should be executed on the accelerators, such as MAGMA [START_REF] Tomov | Towards Dense Linear Algebra for Hybrid GPU Accelerated Manycore Systems[END_REF], but such static schedules are hardly conceivable in a portable fashion. Another solution is to have dividable tasks: when StarPU detects that the granularity is not adapted, it could dynamically perform some local transformations on the graph of tasks, to either divide or merge tasks (similarly to the task chaining mechanism that we use on the Cell).

The difficulty of this implicit method is that StarPU would have to automatically make appropriate decisions about when to change the grain size, and how to select a better granularity. This could be achieved with performance feedback mechanisms, either offline or online. A post-mortem analysis would help in choosing the optimal block size or in detecting the different phases of a regular algorithm. A dynamic approach could use the metrics available from the scheduling engine, such as load imbalance, in collaboration with hardware performance counters (e.g. to detect a sudden contention).

Related Work

Recent years have seen the democratization of accelerator-based computation. A vast majority of the interests has however been devoted to writing efficient code for the accelerators. Relatively few projects have investigated the problem of tasks scheduling. To our knowledge we are the first to study this specific aspect in depth while offering a portable and high-level platform to experiment with task scheduling on such heterogeneous machines equipped with accelerators.

Various runtime systems have been designed (or adapted) to support heterogeneous multicore systems. Most of them are however offering an interface to offload tasks on only a single type of accelerator. While those interfaces hide most low-level technical issues, they also typically require the programmer to decide where the computation should take place, so that some knowledge of the underlying platform is still required.

Jimenez et al. [START_REF] Víctor | Predictive Runtime Code Scheduling for Heterogeneous Architectures[END_REF] consider the problem of scheduling computations on multicore machines equipped with accelerators. However, contrary to StarPU which schedules tasks at a fine granularity, they only dispatch whole applications with respect to the relative speedups observed on the different processing units during previous executions, regardless of the data input. The burden of data management is also left to the programmer. The Charm++ runtime system was also implemented for both GPUs [START_REF] Wesolowski | An application programming interface for general purpose graphics processing units in an asynchronous runtime system[END_REF] and Cell processors [START_REF] Kunzman | Charm++, Offload API, and the Cell Processor[END_REF] in a cluster environment. Even though the Offload API they proposed in [START_REF] Kunzman | Charm++, Offload API, and the Cell Processor[END_REF] was adopted by most task-based approaches, there is currently no performance evaluation available yet, to the best of our knowledge.

Similarly to the low-level interface of OpenCL, IBM ALF [START_REF] Crawford | Accelerating computing with the cell broadband engine processor[END_REF] provides an interface to offload tasks onto heterogeneous platforms, such as the RoadRunner. ALF therefore supports both x86 and Cell processors. While ALF relieves programmers from numerous low-level concerns, it only offers very basic load balancing mechanisms thanks to the DaCS framework, and it is done at the process level contrary to StarPU which actually schedules tasks.

Likewise, the Cell Run Time Library (or Cell-RTL) [START_REF] Nijhuis | Mapping and synchronizing streaming applications on cell processors[END_REF] provides a convenient interface to offload tasks onto the SPUs of a Cell processor. As we have shown in Section 7.2, StarPU and the Cell-RTL are very complementary on this platform. The Cell-RTL indeed provides low-level Cell-specific optimiza-tions (e.g. very efficient synchronization mechanisms) but it only has basic load balancing mechanisms which StarPU can replace. Cell-RTL also handles data movements within a Cell chip, but StarPU provides the higher-level abstractions required to manage them consistently.

Diamos et al. [START_REF] Diamos | Harmony: an execution model and runtime for heterogeneous many core systems[END_REF] present a "collection of techniques" for the Harmony runtime system which they "would like to implement in a complete system". Some of those techniques are similar to those implemented in StarPU: for instance they consider the problem of task scheduling with the support of performance models. Harmony however does not allow user-provided scheduling strategies contrary to StarPU which makes it possible to select the most appropriate policy at runtime. Besides the regression-based model also proposed by Harmony, the tight integration of StarPU's data management library along with the scheduler allows much simpler, yet more accurate, history-based performance models [START_REF] Augonnet | Exploiting the Cell/BE architecture with the StarPU unified runtime system[END_REF]. The data management facilities offered by StarPU are also much more flexible as it is possible to manipulate data in a high-level fashion that is much more expressive than a mere list of addresses.

StarSs introduces #pragma language annotations similar to those proposed by Ayguadé et al. [START_REF] Ayguade | A proposal to extend the openmp tasking model for heterogeneous architectures[END_REF]. They rely on a source-to-source compiler to generate offloadable tasks. Contrary to StarPU, the implementation of the StarSs model is done by the means of a separate runtime system for each platform (SMPSs, GPUSs [START_REF] Ayguadé | An Extension of the StarSs Programming Model for Platforms with Multiple GPUs[END_REF], CellSs [START_REF] Bellens | Exploiting Locality on the Cell/B.E. through Bypassing[END_REF]). Some efforts are done to combine those different runtime systems: Planas et al. allow programmers to include CellSs tasks within SMPSs tasks [START_REF] Planas | Hierchical task based programming with StarSs[END_REF], but this remains the duty of the programmer to decide which tasks should be offloaded. In contrast, StarPU considers tasks that can indifferently be executed on multiple targets.

GPUSs permits to use multiple GPUs [START_REF] Ayguadé | An Extension of the StarSs Programming Model for Platforms with Multiple GPUs[END_REF], but its load balancing mechanisms and data management are still very simple. CellSs also consider advanced task scheduling mechanisms taking data locality into account [START_REF] Bellens | Cellss: Scheduling techniques to better exploit memory hierarchy[END_REF]. Contrary to the Cell-RTL (which we used to offload StarPU tasks onto SPUs), the runtime system of CellSs greatly benefits from software data caching which permits to directly direct transfers between SPUs. Adding this feature into Cell-RTL or using CellSs within StarPU's Cell driver would significantly improve our support for the Cell processor. Likewise, we plan to use StarPU as a target of the source-to-source compiler from StarSs.

The approach of the Anthill runtime system is very similar to StarPU. Teodoro et al. experimented [START_REF] Teodoro | Coordinating the Use of GPU and CPU for Improving Performance of Compute Intensive Applications[END_REF] the impact of task scheduling for clusters of machine equipped with a single GPU and multiple processors. They implemented two scheduling policies which are equivalent to the simple greedy strategies we have shown in Section 5.2. The tight integration of data management and task scheduling within StarPU makes it possible to tackle the problem of task scheduling onto multiple GPUs while taking data movements overhead into account, while Anthill only considers relative speedups to select the most appropriate processing unit.

Conclusion and Future Work

The StarPU runtime system and the experiments presented in this paper can be downloaded at http://runtime.bordeaux.inria.fr/StarPU/. The significant performance improvements that we exhibit on top of already optimized RR n➦ 7240 computation kernels show that task scheduling and efficient code generation are orthogonal problems: producing the most efficient computational kernels is wasteful if those are not scheduled judiciously. We advocate that such scheduling techniques are required to cope with the growing number of accelerators in a scalable fashion: StarPU massively reduces the amount of data transfers without any application code modification.

The key contribution of StarPU is the tight collaboration between its highlevel data management library and its portable scheduling engine which permits to easily design powerful scheduling policies. StarPU unlocks the portability of performance on complex accelerator-based platforms: it is for instance generic enough to transparently handle heterogeneous multi-GPU setups by hiding both low-level heterogeneity and by dispatching tasks according to the capabilities of the different units.

StarPU is not limited to multicore machines equipped with GPUs and Cell processors. Its asynchronous event-driven design will for instance make it straightforward to implement an OpenCL backend. The Larrabee processor should also benefit from the techniques adopted to exploit the Cell processor. Likewise, the dynamic code loading mechanism that we implemented in the Cell StarPU driver would be applicable to reconfigurable architectures such as FPGAs. StarPU would also be helpful in the case of fixed-functions hardware such as multicore DSPs which also raise heavy task and data scheduling requirements.

In the future, we plan to use StarPU as a backend for compiling environments (such as StarSs) which would generate StarPU tasks. We intend to improve our performance models with regards to the actual hardware: first by taking NUMA effects into account when scheduling tasks, and then by considering the non-uniform performance of the IOs on such NUMA platforms [START_REF] Moreaud | Impact of NUMA Effects on High-Speed Networking with Multi-Opteron Machines[END_REF]. We will experiment with better data prefetching and data reclaiming policies to improve the performance in the case of out-of-core computations and data streaming applications. The importance given to the penalty for non-local data accesses should also be dynamically throttled depending on the load balance conditions. Similarly, the scheduling engine should be assisted by some execution feedback (e.g. from the hardware performance counters). Eventually, having different instances of StarPU collaborating through MPI would allow to exploit clusters of machines equipped with accelerators. Granularity is a key issue that we think has yet been ignored in the context of accelerator-based computation until now: we believe that our proposition of dividable tasks is an interesting basis to attack this challenging problem.

More generally, we envision StarPU as a meeting point to leverage the efforts in various research domains: it stands as a powerful backend for compiling environments (eg. StarSs), it offers a high-level playground to experiment with algorithms from the scheduling literature on actual platforms, and it provides portable performance to HPC libraries such as MAGMA.

4 Figure 1 :

 41 Figure 1: Execution of a Task within StarPU.

1 / 4 / 7 / 8 * 12 /Figure 2 :

 1478122 Figure 2: An example of partitioned data, its tree representation and the corresponding StarPU code

Figure 3 : 5 ( 7 } 8 9Figure 4 :

 35784 Figure 3: Expressing task dependencies with tags

Figure 5 :

 5 Figure 5: Manipulating sub-matrices with filters.

Figure 6 :

 6 Figure 6: Blocked matrix multiplication written with StarPU using the codelet defined in Figure 4 and the filters described on Figure 5.

Figure 7 :

 7 Figure 7: The MSI protocol maintains the state of each data on the different memory node. This state (Modified, Shared or Invalid) is updated accordingly to the access mode (Read or Write).

  An example of strategy.

Figure 8 :

 8 Figure 8: All scheduling strategies implement the same queue-based interface.

Figure 9 :Figure 10 :

 910 Figure 9: A pathological case: Gantt diagram of a blocked matrix multiplication with a greedy strategy. CPU CPU FX5800 FX4600

Figure 11 :

 11 Figure 11: Performance feedback loop.

  cl -> cpu_func ( task -> interface , task -> cl_arg ) ; 7 push_task_output ( task ) ;

8 9h 10 }Figure 12 :

 81012 Figure 12: Driver for one CPU core

Figure 13 :

 13 Figure 13: Offloading tasks with the Cell-RTL.

1 FX5800 # 2 Figure 14 :FX5800Figure 15 :

 121415 Figure 14: Distribution of the computation on a multi-GPU system.

Figure 16 :

 16 Figure 16: Impact of priorities on Cholesky decomposition.

Figure 17 :

 17 Figure 17: Impact of scheduling policies on LU decomposition.

Figure 18 :

 18 Figure 18: Penalizing non local data accesses reduces the total amount of data transfers by lowering the ratio of data accesses that require memory transfers.

Figure 19 :

 19 Figure 19: Impact of granularity on the efficiency of LU decomposition on GPUs+CPUs.

Table 1 :

 1 Scheduling strategies implemented with StarPU

	Name	Policy description
	greedy	Greedy policy with support for priorities
	no-prio	Greedy policy without support for priorities
	ws	Greedy policy based on Work Stealing
	w-rand	Random weighted by processor speeds
	heft-tm	HEFT based on Task duration Models
	heft-tm-pr	heft-tm with data PRefetch
	heft-tmdp	heft-tm with remote Data Penalty
	heft-tmdp-pr	heft-tmdp with data PRefetch

Table 2 :

 2 StarPU critical path overhead Task submission

			Data Management	Execution
	1 CPU	468 ± 7.3 ns	94 ± 8 ns	364 ± 17 ns
	8 CPUs 470 ± 8.2 ns	113 ± 11 ns	733 ± 26 ns

Table 3 :

 3 Dynamically adapted granularity for the LU decomposition of a 1GB

	matrix	Block size	(512x512) (1024x1024)	hybrid
		Time (ms)	5991 ± 29	5974 ± 57	5691 ± 83
		Speed (GFlop/s) 489.4 ± 2.4 490.8 ± 4.7 515.2 ± 7.4
	8.5.1 Making trade-offs to select the a granularity

RR n➦ 7240

in the case of Figure6, the cublas SGEMM function RR n➦ 7240

Acknowledgments

This work has been supported by the ANR through the COSINUS (PROHMPT ANR-08-COSI-013 project) and CONTINT (MEDIAGPU ANR-09-CORD-025) programs. We thank NVIDIA and NVIDIA's Professor Partnership Program for their hardware donations.