Etienne Berthier
email: etienne.berthier@legos.obs-mip.fr

Ted A Scambos
email: teds@nsidc.org

Christopher A Shuman
email: christopher.a.shuman@nasa.gov

Mass loss of Larsen B tributary glaciers (Antarctic Peninsula) unabated since 2002

Ice mass loss continues at a high rate among the large glacier tributaries of the Larsen B Ice Shelf following its disintegration in 2002. We evaluate recent mass loss by mapping elevation changes between 2006 and 2010/11 using differencing of digital elevation models (DEMs). The measurement accuracy of these elevation changes is confirmed by a 'null test', subtracting DEMs acquired within a few weeks. The overall 2006-2010/11 mass loss rate (9.0±2.1 Gt a -1 ) is similar to the 2001/02-2006 rate (8.8±1.6 Gt a -1 ), derived using DEM differencing and laser altimetry. This unchanged overall loss masks a varying pattern of thinning and ice loss for individual glacier basins. On Crane Glacier, the thinning pulse, initially greatest near the calving front, is now broadening and migrating upstream. The largest losses are now observed for the Hektoria/Green glacier basin, having increased by 33% since 2006. Our method has enabled us to resolve large residual uncertainties in the Larsen B sector and confirm its state of ongoing rapid mass loss.

Introduction

It is now well-demonstrated that the larger, deeper tributary glaciers of the Larsen A and B ice shelves have dramatically accelerated, retreated and thinned in response to the disintegration events of 1995[START_REF] Rott | Northern Larsen Ice Shelf, Antarctica: further retreat after collapse[END_REF][START_REF] De Angelis | Glacier Surge After Ice Shelf Collapse[END_REF][START_REF] Pritchard | Extensive dynamic thinning on the margins of the Greenland and Antarctic ice sheets[END_REF][START_REF] Rignot | Accelerated ice discharge from the Antarctic Peninsula following the collapse of Larsen B ice shelf[END_REF][START_REF] Rott | Northern Larsen Ice Shelf, Antarctica: further retreat after collapse[END_REF][START_REF] Rott | The imbalance of glaciers after disintegration of Larsen B ice shelf, Antarctic Peninsula[END_REF][START_REF] Scambos | Glacier acceleration and thinning after ice shelf collapse in the Larsen B embayment, Antarctica[END_REF][START_REF] Shuman | 2001-2009 elevation and mass losses in the Larsen A and B embayments, Antarctic Peninsula[END_REF]. Although the collapse of a floating ice shelf has no direct impact on sea level, the increased ice discharge from grounded tributary glaciers does contribute significant mass to the ocean. However, there are still large and unexplained discrepancies (range: 4-22 Gt a -1 ) between the ice losses inferred for the larger tributary glaciers that fed the Larsen B ice shelf (hereafter referred to as the northern Larsen B tributary glaciers, NLBTG, comprising the Hektoria-Green-Evans glacier system, Jorum and Punchbowl glaciers, Crane Glacier, and Mapple, Melville, and Pequod glaciers; see Figure 1) using different assessment methods [START_REF] Rignot | Accelerated ice discharge from the Antarctic Peninsula following the collapse of Larsen B ice shelf[END_REF][START_REF] Rott | The imbalance of glaciers after disintegration of Larsen B ice shelf, Antarctic Peninsula[END_REF][START_REF] Shuman | 2001-2009 elevation and mass losses in the Larsen A and B embayments, Antarctic Peninsula[END_REF].

Nearly a decade after the ice shelf collapsed, reconciling estimates of the NLBTG mass imbalance contribution is essential.

Three methods are presently available to measure the changes in the mass of an ice sheet, or a portion of it: space gravimetry; the mass budget method (MBM); and the geodetic method (GM).

Space gravimetry from the Gravity Recovery and Climate Experiment (GRACE) is currently not able to resolve mass losses occurring at a length scale of a few tens of kilometers, the typical size of the NLBTG basins, and thus has been applied broadly to the Antarctic Peninsula north of 70°S (Graham Land) where a steady decrease in mass is observed since 2002 (-32±6 Gt a -1 [START_REF] Ivins | On-land ice loss and glacial isostatic adjustment at the Drake Passage: 2003-2009[END_REF]; -28.6 Gt a -1 [START_REF] Chen | Accelerated Antarctic ice loss from satellite gravity measurements[END_REF]). The MBM consists of comparing the input (net accumulation) to the output (ice flux through a cross-sectional gate at, or close to, the grounding line). One determination using MBM for all basins between Hektoria and Crane glaciers (Figure 1), indicated mass losses of 21.9±6.6 Gt a -1 in 2003 [START_REF] Rignot | Accelerated ice discharge from the Antarctic Peninsula following the collapse of Larsen B ice shelf[END_REF] (30% uncertainty from [START_REF] Rignot | Changes in ice dynamics and mass balance of the Antarctic ice sheet[END_REF]). Recently, the same technique was applied to the same glaciers by another group using velocity fields measured in 2008/2009 in comparison with velocities from 1995-1999 [START_REF] Rott | The imbalance of glaciers after disintegration of Larsen B ice shelf, Antarctic Peninsula[END_REF]. Despite small changes in surface velocities since 2003 overall, they reported mass losses of 4.1±1.6 Gt a -1 , a factor of 5 lower than the previous MBM study.

A geodetic estimate of the regional mass loss can be obtained by the differencing of digital elevation models (DEMs) acquired a few years apart. Differencing of DEMs acquired prior to the Larsen B ice shelf disintegration in November 2001 (extended to the upper part of Crane Glacier with a DEM acquired in November 2002) with one acquired in November 2006, yields a distinctly different estimate of the NLBTG mass loss, 8.8±1.6 Gt a -1 [START_REF] Shuman | 2001-2009 elevation and mass losses in the Larsen A and B embayments, Antarctic Peninsula[END_REF].

The causes of the differences between the three published NLBTG loss estimates are not understood yet. Likely factors are (i) different time periods surveyed combined with non-steady mass loss response, (ii) uncertainties in net accumulation for the MBM, (iii) unknown bed topography close to the grounding line for most glaciers, (iv) unaccounted surface drawdown at MBM flux gates, (v) errors in the DEMs for the GM; and (vi) grounding line migration for the GM and the MBM. Both the MBM estimates discussed above and our GM study use 900 kg/m 3 as density for converting volume to mass. Using this density (instead of that for pure ice, 917 kg/m 3 ) in the Larsen B sector is justified by the fact that the elevation changes are nearly entirely dynamically-driven, and that the entire column of ice (pure ice, firn, and snow) is being lost by calving at the ice fronts. By our estimate about 80% of the losses occur on fast flowing glacier trunks below 500 m a.s.l.

In this study, we use new 2010 and 2011 satellite DEMs to infer the mass loss of NLBTG between 2006 and 2010/11. These updated mass losses are then compared to 2001/02-2006 losses to determine how NLBTG losses have evolved over time since the break up and, thus, to partly address issue (i). In addition, satellite DEMs acquired within a few weeks in late 2006 are compared to better constrain the uncertainties associated with the GM (issue v). A basin-by-basin comparison is also performed to identify the glaciers for which the discrepancies between MBM and GM estimates are largest and help to target future data acquisition to reconcile estimates of the NLBTG mass imbalance (issue iii). Our elevation changes time series from DEMs (augmented with airborne and satellite laser altimetry) helps to address issue (iv). Grounding line migration since 2002 has probably been rapid, with the ice fronts quickly retreating past its 1990s position [START_REF] Shuman | 2001-2009 elevation and mass losses in the Larsen A and B embayments, Antarctic Peninsula[END_REF]. Satellite imagery suggests the present ice fronts are partially floating, complicating flux gate assessments [START_REF] Zgur | Geophysical survey of the thick, expanded sedimentary fill of the new-born Crane fjord (former Larsen B Ice Shelf, Antarctica)[END_REF][START_REF] Rott | The imbalance of glaciers after disintegration of Larsen B ice shelf, Antarctic Peninsula[END_REF].

Data Sets and Methods

Our analysis is based on seven DEMs of the NLBTG derived from ASTER [START_REF] Fujisada | ASTER DEM performance[END_REF] and SPOT5 [START_REF] Korona | SPIRIT. SPOT 5 stereoscopic survey of Polar Ice: Reference Images and Topographies during the fourth International Polar Year (2007-2009[END_REF] optical stereo-imagery (Table S1, Figure S1 andS2). The processing steps followed to horizontally/vertically adjust the DEMs have been described in detail previously [e.g., [START_REF] Shuman | 2001-2009 elevation and mass losses in the Larsen A and B embayments, Antarctic Peninsula[END_REF] and are briefly summarized here. Cloudy and unreliable pixels are masked. All DEMs are first horizontally co-registered to the reference DEM (the 25 November 2006 SPOT5 DEM) by minimizing the standard deviation of the elevation differences.

Then, all DEMs are vertically adjusted to the 25 November 2006 DEM using the (assumed) stable regions outside of the fast changing outlet glaciers. Only a constant vertical offset is corrected for each DEM around each major basin (the average of the vertical bias within each basin, typically less than 5 m), neglecting any spatial variations in the vertical bias within each basin.

Null tests: accuracy of the mass loss from sequential DEMs

We analyze three DEMs acquired within a short time span in late 2006 (Table S1) to assess the accuracy of the ASTER/SPOT5 and SPOT5/SPOT5 basin-wide elevation changes. The assumption of insignificant elevation change is valid between the SPOT5 and ASTER DEMs as they were acquired only 16 minutes apart on 25 November 2006. This assumption is not as appropriate between the 25 November 2006 and the 31 December 2006 SPOT5 DEMs, with 36 days elapsed.

The rates of elevation changes measured during 2001/02-2006 [Shuman et al., 2011, their Table 2] are used to estimate 36-day elevation changes within each major basin. Based on this, the maximum change is -0.7 m for Evans Glacier and the smallest is for Jorum Glacier at -0.2 m.

Basin-wide mean DEM differences (Table 1) provide null test errors for the GM applied to the NLBTG.

For individual drainage basins, the errors are larger for ASTER/SPOT5 than for SPOT5/SPOT5. This is expected given the higher image resolution and the better orbit knowledge for the SPOT5 sensor so that SPOT5 DEMs are about a factor of two more precise than ASTER DEMs [Berthier et al., 2010, their Table S5]. The basin-wide elevation difference in the ASTER/SPOT5 comparison reaches 5.0 m for Crane Glacier. This result confirms the ±5 m elevation error used previously [START_REF] Shuman | 2001-2009 elevation and mass losses in the Larsen A and B embayments, Antarctic Peninsula[END_REF]. In contrast, the basin-wide elevation difference is always less than 2 m for SPOT5/SPOT5 differences. When the mean elevation difference is computed for all NLBTG basins, errors are 2.4 m for ASTER/SPOT5 and 1.0 m for SPOT5/SPOT5.

In the following assessments of volume and mass change, a ±2 m (respectively, ±5 m) uncertainty is used when two SPOT5 (respectively, ASTER and SPOT5) DEMs are subtracted.

Our null tests indicate that these uncertainties are reasonable at the individual basin scale, and are conservative for elevation differences averaged over all NLBTG. 2).

Elevation changes and mass loss in

Between the two epochs, the maximum mean annual thinning rate considerably diminished for Crane Glacier from over 35 m a -1 to 18 m a -1 . We note that the peak elevation loss during 2001/02-2006 was likely influenced by a sub-glacial lake drainage event that occurred between September 2004 and September 2005 in the lower Crane Glacier [START_REF] Scambos | The triggering of sub-glacial lake drainage during the rapid glacier drawdown: Crane Glacier, Antarctic Peninsula[END_REF], but a region of >20 m a -1 loss extended over a far larger area than the inferred lake extent in the lower glacier trunk.

During 2006-2010/11, thinning has propagated further upstream, and the 10 m a -1 thinning rate contour has moved from 500 m a.s.l. to 700 m a.s.l between the two assessments. An upstream migration of thinning is also observed for Hektoria and Green glaciers with the upper extent of the 10 m a -1 thinning contour moving from 350 to 500 m a.s.l. between the two epochs. For the latter two glaciers, the peak thinning rates are higher during 2006-2010/11 (35 m a -1 ) than during 2001/02-2006 (23 m a -1 ). For the main trunks of Hektoria, Green and Crane glaciers, the horizontal speed of inland propagation of the 10 m a -1 thinning contour is similar, at about 2 km a -1 .

Total mass losses from Hektoria and Green glaciers have increased by one third since the earlier period (from 4.2 to 5.6 Gt a -1 ). Their ice front positions have also varied considerably since 2002, with both retreats and advances [START_REF] Shuman | 2001-2009 elevation and mass losses in the Larsen A and B embayments, Antarctic Peninsula[END_REF]. This makes these glaciers the main contributors to the regional losses (Table 2), ahead of Crane Glacier where net mass loss and ice front position (since ~2006) have remained nearly unchanged.

South of the major NLBTG glaciers, elevation changes of Mapple, Melville and Pequod glaciers are small despite the fact that they have been less constrained since the Larsen B ice shelf disintegration. In our earlier study, we noted that these glaciers have shallower seabed bathymetric troughs in front of them, and likely had less of their longitudinal resistive stresses derived from the former ice shelf [START_REF] Shuman | 2001-2009 elevation and mass losses in the Larsen A and B embayments, Antarctic Peninsula[END_REF][START_REF] Zgur | Geophysical survey of the thick, expanded sedimentary fill of the new-born Crane fjord (former Larsen B Ice Shelf, Antarctica)[END_REF].

Discussion

The null test in Section 3 demonstrates that for all NLBTG (~3000 km²), area-average elevation changes can be measured with an accuracy of ±1 m (using two SPOT5 DEMs) and ±2.5 m (using one ASTER and one SPOT5 DEM). These null test errors can be compared to the standard errors, often applied in differential DEMs studies [e.g., [START_REF] Nuth | Co-registration and bias corrections of satellite elevation data sets for quantifying glacier thickness change[END_REF]. The standard errors are computed from the standard deviation of the DEMs (typically ±10 m for SPOT5 and ±15 m for ASTER) after accounting for the number of independent samples. For all NLBTG basins, and assuming autocorrelation lengths of 200 m [START_REF] Howat | Rates of southeast Greenland ice volume loss from combined ICESat and ASTER observations[END_REF] / 500 m [START_REF] Berthier | Contribution of Alaskan glaciers to sea-level rise derived from satellite imagery[END_REF] / 1000 m [START_REF] Nuth | Glacier geometry and elevation changes on Svalbard (1936-90): a baseline dataset[END_REF], standard error of, respectively ±0.7 m / ±1 m / ±1.5 m are derived for the ASTER/SPOT5 comparison. Those standard errors are all lower than our null test error (±2.5 m), probably due to spatially-varying vertical biases in the DEMs [START_REF] Nuth | Co-registration and bias corrections of satellite elevation data sets for quantifying glacier thickness change[END_REF] that cannot be accounted for by a single shift measured on (assumed) stable regions. Modeling these complex vertical biases in the DEMs using ICESat laser altimetry data was attempted. This is not discussed here because ICESat data are too scarce in our study region to further refine the error bars calculated in the null test. Such a strategy of DEM adjustment using precise external data would be certainly useful in glaciated areas where stable ground is effectively absent, where the density of altimetry data is higher (e.g., closer to 86° latitude for ICESat) and where the altimetry surveys are performed close-in-time to the DEMs.

Geodetic mass losses upstream of the pre-collapse grounding line [START_REF] Rack | Pattern of retreat and disintegration of the Larsen B ice shelf, Antarctic Peninsula[END_REF] 2). This finding is in agreement with limited changes in velocities for most glaciers between [START_REF] De Angelis | Glacier Surge After Ice Shelf Collapse[END_REF]and 2008/2009[START_REF] Rott | The imbalance of glaciers after disintegration of Larsen B ice shelf, Antarctic Peninsula[END_REF], despite some shorter term flow variability for Crane [START_REF] Rignot | Changes in ice dynamics and mass balance of the Antarctic ice sheet[END_REF][START_REF] Scambos | The triggering of sub-glacial lake drainage during the rapid glacier drawdown: Crane Glacier, Antarctic Peninsula[END_REF] and Hektoria [START_REF] Rignot | Changes in ice dynamics and mass balance of the Antarctic ice sheet[END_REF][START_REF] Rott | The imbalance of glaciers after disintegration of Larsen B ice shelf, Antarctic Peninsula[END_REF] glaciers. The consistency of the mass loss through time between 2002 and 2010/11 is also independently supported by the regional GRACE time series [START_REF] Ivins | On-land ice loss and glacial isostatic adjustment at the Drake Passage: 2003-2009[END_REF] and by a continuous elastic uplift of the solid earth since 2002 at the Palmer GPS station, ~100 km west of NLBTG [START_REF] Thomas | Widespread low rates of Antarctic glacial isostatic adjustment revealed by GPS observations[END_REF]. Our analyses reveal an evolving pattern of elevation changes, with a wave of glacier thinning that has broadened and migrated rapidly upstream over time.

Our GM mass loss estimates for NLBTG lie between the losses inferred by the two earlier studies using the MBM, suggesting that losses are overestimated by 12.9±6.9 Gt a -1 in Rignot et al.

[2004] and underestimated by 4.7±2.6 Gt a -1 in [START_REF] Rott | The imbalance of glaciers after disintegration of Larsen B ice shelf, Antarctic Peninsula[END_REF]. Quantifying and understanding these large discrepancies is important because these methods (MBM and GM) are employed to determine the mass balance of glaciers in the Peninsula, a region which alone contribute about one fourth of the continent-wide mass imbalance [e.g., [START_REF] Rignot | Recent Antarctic ice mass loss from radar interferometry and regional climate modelling[END_REF]. It is also crucial to provide realistic scenarios of mass losses following ice shelf collapse to test the ability of ice flow models to simulate them. We note that the GM method requires far fewer assumptions than the MBM in the Larsen B region, with well-constrained differential DEM errors (section 3). Both bed topography and surface mass balance are poorly known for the NLBTG. For Crane Glacier only, bed topography at the present grounding line can be inferred by extrapolating upstream bathymetry data collected in 2006 [START_REF] Rott | The imbalance of glaciers after disintegration of Larsen B ice shelf, Antarctic Peninsula[END_REF][START_REF] Zgur | Geophysical survey of the thick, expanded sedimentary fill of the new-born Crane fjord (former Larsen B Ice Shelf, Antarctica)[END_REF]. The authors of the preceding MBMbased studies had, by necessity, to assume the glaciers to be in equilibrium before the Larsen B ice shelf collapse and computed the pre-collapse ice discharge from a model estimate of net accumulation over the catchment basins. However, recent assessments of surface mass balance (net accumulation) for Antarctica [START_REF] Lenaerts | A new, high-resolution surface mass balance map of Antarctica (1979-2010) based on regional atmospheric climate modeling[END_REF], and a review of some field measurements for the region reported in [START_REF] Rott | The imbalance of glaciers after disintegration of Larsen B ice shelf, Antarctic Peninsula[END_REF], as well as field evidence from ongoing monitoring by one of us (TAS) indicate that the effective net accumulation (~1900 kg m -2 a -1 ) used by [START_REF] Rignot | Accelerated ice discharge from the Antarctic Peninsula following the collapse of Larsen B ice shelf[END_REF] was too large. A strong accumulation gradient exists across the NLBTG area, with progressively decreasing snow input east of the Antarctic Peninsula divide. Measured accumulation rates as part of an ongoing study (LARISSA: Larsen Ice Shelf System Antarctica) reach 2000 to 3000 kg m -2 a -1 at the divide [Zagorodnov et al., in press], but are < 300 kg m -2 a -1 near 450 m a.s.l. on the nearby Flask and Leppard glacier outlets. Applying the basin-wide mean net accumulation value indicated by [START_REF] Rott | The imbalance of glaciers after disintegration of Larsen B ice shelf, Antarctic Peninsula[END_REF], 1087 kg m a -1 , would adjust the mass imbalance reported by [START_REF] Rignot | Accelerated ice discharge from the Antarctic Peninsula following the collapse of Larsen B ice shelf[END_REF] to 57% of the value reported, or ~12.2 Gt a -1 for the NLBTG region (if the accumulation correction ratio holds for the entire area).

We partially reconcile our net imbalance estimate on Crane Glacier with the value reported by [START_REF] Rott | The imbalance of glaciers after disintegration of Larsen B ice shelf, Antarctic Peninsula[END_REF] by considering the velocity variations observed over the past decade [Scambos et al., 2011]. A sequence of visible and near-infrared satellite images shows more than one acceleration period, and in particular ice flow speed in late 2006 was 1.3 times the level in late 2008 (the time of the Rott et al. measurement). Further slowing occurred on Crane Glacier between 2008 and 2009 [START_REF] Rott | The imbalance of glaciers after disintegration of Larsen B ice shelf, Antarctic Peninsula[END_REF][START_REF] Scambos | The triggering of sub-glacial lake drainage during the rapid glacier drawdown: Crane Glacier, Antarctic Peninsula[END_REF]. Thus, our geodetic method, which integrates over a 5-year period, may be expected to show a higher value than a shorter-term assessment during a single period of slower flow speed.

During the two epochs studied here, Hektoria and Green glaciers have been the major contributors to the regional mass loss (over 60% of total loss during 2006-2010/11). This is also where the largest discrepancies in the MBM estimates are observed in our study area. The data most needed to reconcile estimates of the mass loss in the Larsen B embayment are bed topography profiles from ice-penetrating radars for the Hektoria and Green glaciers.

Conclusion

At 8.9 Gt a -1 , our 2002-2011 NLBTG mass loss estimate represents about one third of the overall loss observed with GRACE in the Graham Land of the Antarctic Peninsula [START_REF] Chen | Accelerated Antarctic ice loss from satellite gravity measurements[END_REF][START_REF] Ivins | On-land ice loss and glacial isostatic adjustment at the Drake Passage: 2003-2009[END_REF]]. An implication is that rapid ice loss and surface lowering are occurring elsewhere in the northern Antarctic Peninsula as indicated by other studies [START_REF] Glasser | From ice-shelf tributary to tidewater glacier: continued rapid recession, acceleration and thinning of Rohss Glacier following the 1995 collapse of the Prince Gustav Ice Shelf, Antarctic Peninsula[END_REF][START_REF] Pritchard | Widespread acceleration of tidewater glaciers on the Antarctic Peninsula[END_REF][START_REF] Pritchard | Extensive dynamic thinning on the margins of the Greenland and Antarctic ice sheets[END_REF]. Our results suggest that differential DEM analysis would provide similar insights on the mass balance for all of Graham Land and similar glaciated regions with mostly unknown ice thickness and spatially varying surface mass balance.

We have re-assessed the glaciers of the Larsen B embayment where variations in published mass balance assessments are largest (up to 300%) and suggest that poorly constrained bedrock topography and net accumulation are the reasons for this range, primarily affecting the MBM method. Our differential DEM analysis shows continuing steady net losses from the Larsen B embayment glaciers overall but with accelerating losses for the northern Hektoria/Green basin. Table 1. Basin-wide mean elevation differences (m) between DEMs acquired the same day (25 Nov. 2006, ASTER/SPOT5) and 36-days apart (SPOT5 25 Nov./SPOT5 31 Dec.). In bold, the maximum absolute error for each method. The last row shows the mean differences for the entire study region. MMP stands for the Mapple, Melville, and Pequod glaciers.

Tables

  2006-2010/11 Recent elevation changes for NLBTG are measured by subtracting two SPOT5 DEMs (31 December 2006 and 14 March 2011) for a 2633 km² area in the west and a SPOT5 (25 November 2006) DEM and an ASTER (15 December 2010) DEM for a 906 km² area in the east (Figure 1). The rates of elevation changes during 2006-2010/11 are compared to those reported in Shuman et al. [2011] during 2001/02-2006 (Figures 2 and 3, Table

  are virtually unchanged since 2002 for the grounded NLBTG (8.8±1.6 Gt a -1 for 2001/02-2006 and 9.0±2.1 Gt a -1 for 2006-2010/11, Table

ASTER*FiguresFigure 1 .

 1 Figures

  Figure 1. Study area and the main drainage basins of the northern Larsen B embayment. All glaciers studied here flowed into the Larsen B ice shelf before its collapse in 2002. Southern tributaries (e.g., Flask and Leppard glaciers) still constrained by a remnant of the Larsen B ice shelf are not considered. The area shaded in light blue is where elevation changes are measured during 2006-2011; the area shaded in yellow has difference measurements during 2006-2010. Background: Mosaic of 25 November 2006 and 31 December 2006 SPOT5 images (© CNES 2006, Spot Image). Inset: location of the study area in the north of the Antarctic Peninsula.

Figure 2 .

 2 Figure 2. Rate of ice elevation changes (m a -1 ) between 2001/02-2006 (left) and 2006-2010/11(right). The 10 m a -1 thinning contour is shown with a gray dotted line, the 20 m a -1 with a dark dotted line.

Figure 3 .

 3 Figure 3. Rate of elevation change averaged for 50-m altitude bands for Crane (a) and Hektoria/Green (b) during 2001/02-2006 (filled symbols) and 2006-2010/11 (unfilled symbols). For clarity, errors bars are only shown in the legends. Insets are 3D surface views showing the glaciers in 25 November 2006 (Copyright CNES 2006, Spot Image).

Acknowledgements

We thank G. Durand, B. Kulessa and four anonymous referees for their comments on earlier versions of the manuscript. EB acknowledges support from CNES (TOSCA and ISIS proposals), PNTS, and ANR-09-SYSC-001. TAS acknowledges support from NSF-OPP, ANT-0732921, and NASA NNX10AR76G grants. CAS was supported by grants from the NASA Cryospheric

Reference [START_REF] Rignot | Accelerated ice discharge from the Antarctic Peninsula following the collapse of Larsen B ice shelf[END_REF][START_REF] Rott | The imbalance of glaciers after disintegration of Larsen B ice shelf, Antarctic Peninsula[END_REF][START_REF] Shuman | 2001-2009 elevation and mass losses in the Larsen A and B embayments, Antarctic Peninsula[END_REF] This


Camille Couprie

Hugues Talbot

Jean-Christophe Pesquet

Laurent Najman

Leo Grady

DUAL CONSTRAINED TV-BASED REGULARIZATION

Keywords: Proximal algorithm, inverse problems, image denoising, convex optimization, image restoration

Algorithms based on the minimization of the Total Variation are prevalent in computer vision. They are used in a variety of applications such as image denoising, compressive sensing and inverse problems in general. In this work, we extend the TV dual framework that includes Chambolle's and Gilboa-Osher's projection algorithms for TV minimization in a flexible graph data representation by generalizing the constraint on the projection variable. We show how this new formulation of the TV problem may be solved by means of a fast parallel proximal algorithm, which performs better than the classical TV approach for denoising, and is also applicable to inverse problems such as image deblurring.

INTRODUCTION

The Total Variation (TV) model was introduced as a regularizing criterion for image denoising by Rudin, Osher and Fatemi [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF], and has been shown to be quite efficient for smoothing images while preserving contours. Moreover, a major advantage is that TV minimization is a convex problem, although coping with the non-differentiability of the involved objective function required the development of specific optimization techniques, particularly with respect to speed and efficiency. Much progress has been achieved by employing primal-dual approaches [START_REF] Chan | A nonlinear primaldual method for total variation-based image restoration[END_REF][START_REF] Carter | Dual Methods for Total Variation-based Image Restoration[END_REF], and more recently by the dual approach of Chambolle [START_REF] Chambolle | An algorithm for total variation minimization and applications[END_REF]. Qualitative improvements have also been obtained by introducing a weighted model [START_REF] Bougleux | Discrete regularization on weighted graphs for image and mesh filtering[END_REF][START_REF] Giboa | Nonlocal operators with application to image processing[END_REF]. In this model, the discretization of the TV energy is performed by the use of an edge-weighted graph. Gilboa and Osher provided an efficient dual algorithm in [START_REF] Giboa | Nonlocal operators with application to image processing[END_REF], based on Chambolle's [START_REF] Chambolle | An algorithm for total variation minimization and applications[END_REF] method to address this problem. Currently, one of the fastest method for optimizing weighted TV is proposed by Zhang et. al [START_REF] Zhang | Bregmanized nonlocal regularization for deconvolution and sparse reconstruction[END_REF].

The projection algorithms in [START_REF] Chambolle | An algorithm for total variation minimization and applications[END_REF][START_REF] Giboa | Nonlocal operators with application to image processing[END_REF] are based on a relatively simple local constraint on the norm of the projection variable. In this work, we extend the constraint on this variable, as this allows us to better adapt the optimization procedure to the local information. We name this new approach "Dual-Constrained Total Variation" (DCTV) regularization.

The paper is organized as follows: in Section 2, we generalize TV models on a graph by extending the dual formulation under a constrained form. In Section 3, we demonstrate how our constrained TV-based optimization problem can be efficiently solved by using a parallel proximal algorithm. Finally, results obtained with the proposed approach are presented in three applications and compared to classical weighted TV.

GRAPH EXTENSION OF TV MODELS

Given an original corrupted image f , the purpose of variational methods for image restoration is to deduce a restored image u close to the observed image f under the assumption of smooth variations of intensity values inside objects. Let λ ∈]0, +∞[ be a real positive value representing a regularization parameter. In a continuous setting, given a plane domain Ω, and denoting by x and y two arbitrary points of Ω, the weighted anisotropic TV model [START_REF] Giboa | Nonlocal operators with application to image processing[END_REF] is given by

min u Z Ω " Z Ω wx,y(uy -ux) 2 dy " 1/2 dx + 1 2λ Z Ω (ux -fx) 2 dx, ( 1
)
where w is a nonnegative valued function defined on Ω 2 . As shown by Chan et al in [START_REF] Chan | A nonlinear primaldual method for total variation-based image restoration[END_REF] the TV minimization problem (1) is equivalent to the min-max problem

min u max ||p||∞≤1 Z Ω 2 w 1/2
x,y (uyux)px,ydxdy+

1 2λ Z Ω (ux -fx) 2 dx, (2)
with p a projection vector field, that include weights on both the edges and nodes. An edge weight is a value assigned to each edge e i,j , and is denoted by w i,j . We assume that w i,j ∈]0, +∞[. The weight of a node v i is denoted by g i . We also assume that g i ∈ R * + . We note A the incidence matrix of the graph which is known to define the discrete calculus analogue of the gradient, while A similarly defines the discrete calculus analogue of the divergence (see [START_REF] Grady | Discrete Calculus: Applied Analysis on Graphs for Computational Science[END_REF] and the references therein). The incidence matrix A ∈ R m×n is defined, for every vertex v k and edge e ij , as

p ∞ = sup x∈Ω " R Ω p 2 x,
A eij ,v k = 1 if i = k, A eij ,v k = -1 if j = k,
and A eij ,v k = 0 otherwise. For any matrix M , we use |M | for the matrix composed of the absolute value of each entry individually, and denote • the Hadamard product. As the step size of the graph tends toward zero, the discrete version of (1) and its dual (2) approximate the continuous versions. Now, f and u are vectors in R n . The discrete TV model as defined by [START_REF] Bougleux | Discrete regularization on weighted graphs for image and mesh filtering[END_REF][START_REF] Giboa | Nonlocal operators with application to image processing[END_REF] is

min u n i=1 j∈Ni w i,j (u j -u i ) 2 1/2 + 1 2λ n i=1 (u i -f i ) 2
(3) where N i = {j ∈ {1, . . . , n} | e i,j ∈ E} and the dual formulation optimized by a projection algorithm [START_REF] Giboa | Nonlocal operators with application to image processing[END_REF] is given by

min u max p ∞ ≤1 p ((Au) • √ w) + 1 2λ u -f 2 (4)
where • denotes the Euclidean norm, w ∈ R m is a vector with components (w i,j ) i,j and denoting by (p i,j ) i,j the com-

ponents of p ∈ R m , p ∞ = max i∈{1,...,n} j∈Ni p 2 i,j 1/2
.

By introducing a vector F = (F i,j ) i,j ∈ R m such that, for every i ∈ {1, . . . , n} and j ∈ N i , F i,j = p i,j √ w i,j , the problem can be reformulated as

min u max F ∈C F Au + 1 2λ u -f 2 , ( 5
)
C = {(F i,j ) i,j | (∀i ∈ {1, . . . , n}) j∈Ni F 2 i,j w i,j ≤ 1}. (6)
The objective of this work is to extend the discrete weighted TV variational formulation in ( 5) by investigating the following optimization problem

min u max F ∈C F Au + 1 2 (Hu -f ) Λ -1 (Hu -f ) (7)
where f ∈ R q is an observed vector of data, H ∈ R q×n and Λ is weighting symmetric definite-positive matrix in R q×q . H may simply be the identity matrix for image denoising; a convolution operator in restoration tasks; or a projection matrix in reconstruction problems. Λ may be a matrix proportional to the covariance matrix of the noise corrupting the data as commonly used in weighted least squares approaches. In addition, the main contribution of this work is to consider more general convex sets C than those given by ( 6). More precisely, the proposed optimization approach allows us to address nonempty convex sets C which can be decomposed as an intersection of closed convex subsets (C r ) 1≤r≤s of R m , the projections onto which take closed forms. An example of a set C of interest is

C = {F | (∀i ∈ {1, . . . , n}) θ i • F α ≤ g i } (8)
where • α is the α norm of R m with α ∈ [1, +∞] and, for every i ∈ {1, . . . , n}, θ i ∈]0, +∞[ m is a vector of multiplicative constants. The form of C in (8) obviously includes (6) as a particular case while offering much more flexibility. In this paper, we are mostly interested in the case when θ i is the i-th line vector of |A | and the i-th node weight g i is a decreasing function of the image gradient. Given positive reals and β, we suggest using

g i = exp(-β ∇u i 2 ) + , ( 9
)
where u is some reference image (wich corresponds to some rough estimate of u) and ∇u i 2 is the Euclidean norm of its discrete gradient ∇u i at node i. In addition to intensity information, g i may be used to penalize changes in other relevant image quantities such as color or texture. In the absence of a contour, g i takes large values, so are the components of F corresponding to nonzero values of θ i , preventing large local variations of u in the minimization [START_REF] Giboa | Nonlocal operators with application to image processing[END_REF]. Conversely, in the presence of a contour, g i , and the components of F corresponding to nonzero values of θ i are small, thus allowing large local variations of u.

PROPOSED ALGORITHM

We show in this section that it is possible to efficiently solve Problem (7) by proximal methods [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF]. To do so, define the support function σ C of the closed convex constraint set C as

σ C : R m →] -∞, +∞]: a → sup F ∈C F a. ( 10
)
This is a proper lower-semicontinuous convex function, the conjuguate of which is the indicator function of C,

ı C : F → 0 if F ∈ C, +∞ otherwise. ( 11
)
This leads us to consider the following optimization problem:

min u σ C (Au) + 1 2 (Hu -f ) Λ -1 (Hu -f ) + η 2 Ku 2
(12) where η ∈]0, +∞[ and K ∈ R n×n is the projection matrix onto the nullspace of H. When H is injective (rank H = n), the last term vanishes and ( 12) is strictly equivalent to [START_REF] Giboa | Nonlocal operators with application to image processing[END_REF]. The term u → η Ku 2 /2 thus aims at introducing an additional regularization when H is not injective, so that the objective function remains strictly convex. The following holds: Proposition 3.1 Problem (12) admits a unique solution u.

The dual Fenchel-Rockafellar form of the problem is min

F ϕ(F ) + ı C (F ), ( 13
)
where ϕ : F → 1 2 F AΓA F -F AΓH Λ -1 f and Γ = (H Λ -1 H + ηK) -1 . The optimal solution to the primal problem (12) is deduced from any optimal solution F of the dual problem by the relation

u = Γ H Λ -1 f -A F . ( 14
)
Note that the dual forward-backward algorithm proposed in [START_REF] Combettes | Dualization of signal recovery problems[END_REF] is not applicable here since the projection onto C is not explicit. In order to numerically solve (13), recall that C = ∩ s r=1 C r , so that ı C can be decomposed into the sum of the indicator functions of the convex subsets (C r ) 1≤r≤s . Hence, the problem is equivalent to solving

min F s r=1 ı Cr (F ) + ϕ(F ). ( 15
)
The above sum of (s + 1) convex functions can be efficiently optimized by resorting to a parallel proximal algorithm [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF].

As shown in Algorithm 1, this requires to compute in parallel projections onto each set C r with r ∈ {1, . . . , s}, which are defined ∀F ∈ R m as P Cr (F ) = arg min Φ∈Cr Φ -F . Note that the convergence of the sequence (F k ) k generated by this algorithm to a solution F of ( 13) is guaranteed, which allows us to deduce a solution to (12) by using Relation (14).

PC r (y r,k ) if r ≤ s (γAΓA + I) -1 (γAΓH Λ -1 f + y s+1,k ) otherwise z k = 2 s+1 (π 1,k + • • • + π s+1,k ) -F k for r = 1, . . . , s + 1 do in parallel y r,k+1 = y r,k + ν(z k -p r,k ) F k+1 = F k + ν 2 (z k -F k ) until convergence
For a 4-connected lattice, when using [START_REF] Grady | Discrete Calculus: Applied Analysis on Graphs for Computational Science[END_REF] where, for every i ∈ {1, . . . , n}, θ i is the i-th line of matrix |A |, the projection onto C is not explicit. However, a decomposition of C can be performed by setting s = 2 and, for every r ∈ {1, 2},

C r = {F | (∀i ∈ S r ) θ i • F α ≤ g i }.
The partition (S r ) 1≤r≤2 corresponds to two spatial disjoint sets each with a checkered pattern. The projection onto C r for α = 2 (resp. α = 1 or α = ∞) reduces to simple projections onto hyperspheres (resp. hypercubes [START_REF] Van Den | Probing the pareto frontier for basis pursuit solutions[END_REF]). Note that the computation of π s+1,k requires a matrix inversion which can be efficiently performed by Fast Fourier Transform when H and Λ are (or can be approximated by) circulant-block circulant matrices.

RESULTS

We now demonstrate the performance of DCTV with respect to weighted TV. In our experiments, we compare solutions to the TV problem (3) given by the state-of-the-art augmented Lagrangian algorithm (a.k.a. split Bregman) of [START_REF] Zhang | Bregmanized nonlocal regularization for deconvolution and sparse reconstruction[END_REF], to solutions optimizing DCTV with our parallel proximal algorithm. We used four standard test images that we corrupted with synthetic Gaussian noise of variance σ 2 . The weights have been set according to [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF] with β = 0.04 for DCTV. The value of the image fidelity parameter λ was set according to an empirical rule depending of the variance of the noise. The value of was set to λ, and the stopping criteria for both algorithms is ||u k -u k-1 || ||u k || < 5.10 -3 . Signal to Noise Ratio (SNR) is used as performance measure in our quantitative evaluation. Table 1 reports SNR values for DCTV and weighted TV results obtained on each image corrupted with different values of noise variance. Examples of results are shown in Figure 1. Those experiments show that DCTV leads to improved results when the variance of the noise is lower than 50. Visually, DCTV are sharper and feature better contrast than the weighted TV results. This also explains the slight degradation of performance in presence of heavy noise (variances 50-100). Geometrically, the improvement of DCTV over the weighted TV can be interpreted in the following way: the In terms of computation time, DCTV is competitive with the most efficient weighted TV algorithm. Denoising the 512 × 512 Lena image corrupted with Gaussian noise (σ 2 = 15) requires 0.38 seconds for split Bregman, versus 0.7 seconds for PPXA to converge on an Intel Xeon 2.5GHz 8-core system. Figures 2 and3 show the ability of DCTV to easily generalize to applications beyond denoising.

CONCLUSION

In this paper we have extended existing TV models by generalizing the constraint on the projection variable of the dual TV formulation. This new approach yields improved results compared with the weighted TV approach in image restoration applications. More generally, the proposed algorithm makes it possible to efficiently solve convex minimization problems involving the support function of an intersection of convex sets as a penalty term. It is also worth emphasizing that this approach can be applied on graph data structures such as those frequently employed in 3D modelling. Future work will further improve our results by using image patches in the weight computation before applying DCTV, and applications on embedded manifolds such as triangulated surfaces.

Algorithm 1 :

 1 Parallel proximal algorithm solving (15) Fix γ > 0 and ν ∈]0, 2[. Set k = 0. Choose y 1,0 = y 2,0 = . . . = y s+1,0 ∈ R m and F 0 . repeat for r = 1, . . . , s + 1 do in parallel π r,k =

Fig. 1 .

 1 Fig. 1. Denoising 'Barbara' image corrupted with Gaussian noise of variance σ 2 = 15

Fig. 2 .

 2 Fig. 2. Image fusion from a an image corrupted with Gaussian noise of variance σ 2 = 20 (b) and a blurry image (c).

  (a) Original image (b) Corrupt SNR=12.3dB (c) DCTV SNR=17.2dB

Fig. 3 .

 3 Fig. 3. Denoising and deblurring an MRI image corrupted with synthetic uniform 7×7 blur and Gaussian noise (σ 2 = 10)

Table 1 .

 1 Quantitative denoising experiment on standard images corrupted with additive Gaussian noise of variance σ 2 .

	σ 2	5	10	15	20	25	50	100
	SNR values obtained by optimizing weighted TV with Split Bregman
	house	22.7	18.9	17.1	15.8	15.1	11.8	7.8
	man	23.3	19.6	17.8	16.5	15.5	12.3	8.3
	lena	22.7	18.9	17.2	16.2	15.4	12.5	8.3
	barbara 20.3	16.0	13.4	11.7	10.5	8.5	6.2
	mean	22.3	18.4	16.4	15.1	14.1	11.3	7.6
		SNR values obtained by optimizing DCTV
	house	22.7	19.3	17.4	15.8	15.2	11.8	7.4
	man	23.7	20.3	18.2	16.7	15.6	12.1	8.1
	lena	23.2	19.7	17.7	16.3	15.5	12.4	8.0
	barbara 20.3	16.2	13.8	12.2	10.9	8.4	6.0
	mean	22.5	18.9	16.8	15.3	14.3	11.2	7.4

convex we use for projection adapts itself to the local neighborhood, and this is reducing blur as a result.


Feng Li
email: feng.li@inria.fr

Antoniu Pop
email: antoniu.pop@mines-paristech.fr

Albert Cohen
email: albert.cohen@inria.fr

Advances in Parallel-Stage Decoupled Software Pipelining Leveraging Loop Distribution, Stream-Computing and the SSA Form

Keywords: D.3.4 [Programming Languages]: Processors-Compilers, Optimization General Terms optimization automatic parallelization, stream-computing, loop distribution

Decoupled Software Pipelining (DSWP) is a program partitioning method enabling compilers to extract pipeline parallelism from sequential programs. Parallel Stage DSWP (PS-DSWP) is an extension that also exploits the data parallelism within pipeline filters.

This paper presents the preliminary design of a new PS-DSWP method capable of handling arbitrary structured control flow, a slightly better algorithmic complexity, the natural exploitation of nested parallelism with communications across arbitrary levels, with a seamless integration with data-flow parallel programming environments. It is inspired by loop-distribution and supports nested/structured partitioning along with the hierarchy of control dependences. The method relies on a data-flow streaming extension of OpenMP.

These advances are made possible thanks to progresses in compiler intermediate representation. We describe our usage of the Static Single Assignment (SSA) form, how we extend it to the context of concurrent streaming tasks, and we discuss the benefits and challenges for PS-DSWP.

Introduction

In recent years, the CPU manufacturers have embraced chip multiprocessors because of technology, power consumption and thermal dissipation constraints, and because of diminishing returns in instruction-level parallelism. The amount of performance gained by the use of multicore processor depends highly on the software fraction that can be parallelized to run on multiple cores simultaneously. Multiprocessor programming leaves the burden to programmer who faces the extra complexity, heisenbugs, deadlocks and other problems associated with parallel programming. The situation is worse when dealing with the migration of legacy code.

Decoupled Software Pipelining (DSWP) is an automatic thread partitioning method which could partition a sequential program to run on multiple cores, and Parallel-Stage DSWP (PS-DSWP) exposes data parallelism into task pipelines extracted by DSWP. These automatic thread partitioning methods free the programmer from manual parallelization. They also promise much wider flexibility than data-parallelism-centric methods for processors, aiming for the effective parallelization of general-purpose applications.

In this paper, we provide another method to decouple controlflow regions of serial programs into concurrent tasks, exposing pipeline and data parallelism. The power and simplicity of the method rely on the restriction that all streams should retain a synchronous semantics [START_REF] Halbwachs | The synchronous dataflow programming language lustre[END_REF]. It amounts to checking the sufficient condition that the source and target of any decoupled dependence are control-dependent on the same node in the control dependence tree (this assumes structured control flow). This restriction may appear as a severe one for experienced parallel programmers; but at the potential expense of adding extra levels of nested parallelism, it does not restrict the degree of pipeline parallelism. In fact, any pair of computational statements can be decoupled and assigned to different concurrent tasks. The partitioning algorithms also handle DOALL parallelization within task pipelines, and arbitrarily nested data-parallel pipelines following the control dependence tree of a structured control flow graph. Unlike existing DSWP algorithms, our method does not explicitly copy conditional expressions and can handle arbitrary backward data and control dependences.

We are using two intermediate representations.

• A conventional SSA-based representation, annotated with natural loop and control dependence trees (for structured control flow).

• And a streaming data-flow extension of the latter representation as a backend for our partitioning algorithm, still in SSA form but with explicit task boundaries (for single-entry single-exit regions) and multi-producer multi-consumer streams to communicate across tasks.

The backend representation streamlines the decoupling of multi-producer multi-consumer data flow through explicit, compilercontrolled sampling and merging stages. Multi-producer multiconsumer semantics is absolutely essential to handle general decoupling patterns where data-parallel stages feature an unbalance in the number of worker threads. Sampling is handled transparently by nesting tasks into enclosing control flow. Merging is captured by Φ functions at task boundaries, introducing a minor variant of the SSA form satisfying the so-called task-closed property that multiple incoming flows targeting the same use in a given task should be explicitly merged by a dedicated Φ function at the task entry point.

Relying on SSA avoids building the complete program dependence graph: with the exception of the array dependence graph, our method only processes linear-size data structures, as opposed to the worst-case quadratic program dependence graph in DSWP.

Related Work

The most closely related work to this paper is decoupled software pipelining and loop distribution. We recall the state-of-theart in both and present the original finding at the source of this work: by extending loop distribution with pipelining and asserting a synchronous concurrency hypothesis, arbitrary data and control dependences can be decoupled very naturally with only minor changes to existing algorithms that have been proposed for loop distribution [START_REF] Kennedy | Loop distribution with arbitrary control flow[END_REF].

Decoupled software pipelining

Decoupled Software Pipelining (DSWP) [START_REF] Ottoni | Automatic thread extraction with decoupled software pipelining[END_REF] is one approach to automatically extract threads from loops. It partitions loops into long-running threads that communicate via inter-core queues. DSWP builds a Program Dependence Graph (PDG) [START_REF] Ferrante | The program dependence graph and its use in optimization[END_REF], combining control and data dependences (scalar and memory). Then DSWP introduces a load-balancing heuristic to partition the graph according to the number of cores, making sure no recurrence spans across multiple partitions. In contrast to DOALL and DOACROSS [START_REF] Cytron | Doacross: Beyond vectorization for multiprocessors[END_REF] methods which partition the iteration space into threads, DSWP partitions the loop body into several stages connected with pipelining to achieve parallelism. It exposes parallelism in cases where DOACROSS is limited by loop-carried dependences on the critical path. And generally speaking, DSWP partitioning algorithms handles uncounted loops, complex control flow and irregular pointerbased memory accesses.

Parallel-Stage Decoupled Software Pipelining [START_REF] Raman | Parallel-stage decoupled software pipelining[END_REF] (PS-DSWP) is an extension to combine pipeline parallelism with some stages executed in a DOALL, data-parallel fashion. For example, when there are no dependences between loop iterations of a DSWP stage, the incoming data can be distributed over multiple data-parallel worker threads dedicated to this stage, while the outgoing data can be merged to proceed with downstream pipeline stages.

These techniques have a few caveats however. They offer limited support for decoupling along backward control and data dependences. They provide a complex code generation method to decouple dependences among source and target statements governed by different control flow, but despite its complexity, this method remains somewhat conservative.

By building the PDG, DSWP also incurs a higher algorithmic complexity than typical SSA-based optimizations. Indeed, although traditional loop pipelining for ILP focuses on innermost loops of limited size, DSWP is aimed at processing large control flow graphs after aggressive inter-procedural analysis optimization. In addition, the loops in DSWP are handled by the standard algorithm as ordinary control flow, missing potential benefits of treating them as a special case. To address these caveats, we turned our analysis to the state of the art in loop distribution.

Loop distribution

Loop distribution is a fundamental transformation in program restructuring systems designed to extract data parallelism for vector or SIMD architectures [START_REF] Kennedy | Loop distribution with arbitrary control flow[END_REF].

In its simplest form, loop distribution consists of breaking up a single loop into two or more consecutive loops. When aligning loop distribution to the strongly connected components of the datadependence graph, one or more of the resulting loops expose iterations that can be run in parallel, exposing data parallelism. Barriers are inserted after the parallel loops to enforce precedence constraints with the rest of the program. An example is presented in Figure 1.

for (i = 1; i < N; i++) { S1 A[i] = B[i] + 1; S2 C[i] = A[i-1] + 1; } for (i = 1; i < N; i++) S1 A[i] = B[i] + 1; <barriers inserted here> for (i = 1; i < N; i++) S2 C[i] = A[i-1] + 1

OpenMP Extension for Stream-Computing as a Code Generation Target

A recently proposed stream-computing extension to OpenMP [START_REF] Pop | A stream-comptuting extension to OpenMP[END_REF] allows the expression of pipeline parallelism by making explicit the flow dependences, or producer-consumer patterns, between OpenMP tasks. It provides a simple way for explicitly building dynamic task graphs, where tasks are connected through streams that transparently privatize the data. The extension consists of two additional clauses, input and output to the task construct, that define the producer-consumer relationships between tasks. The OpenMP language, with this extension, is a natural fit as a target for our code generation. It provides for dynamic task creation and connection in the task graph, it handles arbitrary nesting of pipelined tasks in control-flow, and it allows the hierarchical nesting of tasks.

The task construct is extended with input and output clauses as presented on Figure 2. Both clauses take a list of items, each of which describes a stream and its behavior w.r.t. the task to which the clause applies. In the abbreviated item form, stream, the stream can only be accessed one element a time through the same variable s. In the second form, stream >> window, the programmer uses the C++ flavoured << >> stream operators to connect a sliding window to a stream, gaining access, within the body of the task, to horizon elements in the stream.

One of the main issues that needs to be addressed in order to distribute a PDG to the OpenMP stream-computing extension is that, in the latter, the data flow bypasses the control flow. In other words, when a task produces values on an output stream, these values will all reach the consumers of the stream, even if, in the serial semantics, the values would have been overwritten before reaching the consumers. This means that the only case where a direct annotation scheme will work is if all tasks are in the same control flow. There are multiple ways this issue can be handled, the most systematic one being to always ensure that every producerconsumer pair share the same control dependence. This is achieved by sinking all control flow surrounding the tasks, and not shared by both producer and consumer, in the tasks. To avoid the loss of parallelization opportunities, each task's body can be further partitioned into nested pipelines.

The GCC implementation of the OpenMP extension for streamcomputing has been shown to be efficient to exploit mixed pipelineand data-parallelism, even in dynamic task graphs [START_REF] Pop | A stream-comptuting extension to OpenMP[END_REF]. It relies on compiler and runtime optimizations to improve cache locality and relies on a highly efficient lock-free and atomic operation-free synchronization algorithm for streams.

Observations

It is quite intuitive that the typical synchronization barriers in between distributed data-parallel loops can be weakened, resulting into data-parallel pipelines. We aim to provide a comprehensive treatment of this transformation, generalizing PS-DSWP in the process.

Replacing loops and barriers with a task pipeline

In the previous example, we could remove the barriers between two distributed loops with pipelining so that the two loops could run in parallel.

(i = 1; i < N; i++) { S1 A[i] = B[i] + 1; produce(stream, A[i]); } for (i = 1; i < N; i++) { tmp = consume(stream); S2 C[i] = tmp + 1; } Figure 3.
Pipelining inserted between distributed loops. Initialize the stream (left), producer and consumer thread (right).

Figure3 shows that pipelined execution is possible: the INIT_STREAM function inserts one delay into a communication stream; the produce/consume primitives implement a FIFO, enforcing the precedence constraint of the data dependence on array A and communicating the value in case the hardware needs this information.

When distributing loops, scalar and array expansion (privatization) is generally required to eliminate memory-based dependences. The conversion to a task pipeline avoids this complication through the usage of communication streams. This transformation can be seen as an optimized version of scalar/array expansion in bounded memory and with improved locality [START_REF] Pop | Automatic streamization in GCC[END_REF].

Extending loop distribution to PS-DSWP

The similarity between DSWP and distributed loops with dataparallel pipelines is striking. First, both of them partition the loop into multiple threads. Second, both of them avoid partitioning the loop iteration space: they partition the instructions of the loop body instead. But four arguments push in favor of refining DSWP in terms of loop distribution.

1. Loop distribution leverages the natural loop structure, where the granularity of thread partitioning can be easily controlled. Moreover, it is useful to have a loop control node to which to attach information about the iteration of the loop, including closed forms of induction variables; this node can also be used to represent the loop in additional transformations.

2. Using a combination of loop distribution and fusion, then replacing barriers with pipelining leads to an incremental path in compiler construction. This path leverages existing intermediate representations and loop nest optimizers, while DSWP relies on new algorithms and a program dependence graph.

3. Considering the handling of control dependences, a robust and general algorithm already exists for loop distribution. McKinley and Kennedy's technique handles arbitrary control flow [START_REF] Kennedy | Loop distribution with arbitrary control flow[END_REF] and provides a comprehensive solution. The same methods could be applied for DSWP, transforming control dependences into data dependences, and storing boolean predicates into stream. After restructuring the code, updating the control dependence graph and data dependence graph, the code generation algorithm for PDGs [START_REF] Baxter | The program dependence graph and vectorization[END_REF][START_REF] Ferrante | On linearizing parallel code[END_REF][START_REF] Ferrante | Generating sequential code from parallel code[END_REF] can be used to generate parallel code. This solution would handle all cases where the current DSWP algorithm fails to clone a control condition.

4. Since loop distribution does not partition the iteration space, it can also be applied to uncounted loops. Unfortunately, the termination condition needs to be propagated to downstream loops. This problem disappears through the usage of a conventional communication stream when building task pipelines.

From this high-level analysis, it appears possible to extend loop distribution with pipelining to implement PS-DSWP and handle arbitrary control dependences. Yet the method still seems rather complex, especially the if-conversion of control dependences and the code generation step from the PDG. We go one step further and propose a new algorithm adapted from loop distribution but avoiding these complexities.

Motivating example

Our method makes one more assumption to reduce complexity and limit risks of overhead. It amounts to enforcing the synchronous hypothesis on all communicating tasks in the partition [START_REF] Halbwachs | The synchronous dataflow programming language lustre[END_REF]. A sufficient condition is to check if the source and target of any decoupled dependence is dependent on the same control node.

Consider the example in Figure 4. S1 and S7 implement the loop control condition and induction variable, respectively. S2, S3 and S6 are control dependent on S1. S3 is a conditional node, S4, S5 and L1 are control dependent on it. In the inner loop, L2 and L3 are control dependent on L1. When we apply DSWP to the outer loop, the control dependences originating from S1 must be if-converted by creating several streams (the number of streams depends on the number of partitions). When decoupling along the control dependence originating from S3, a copy of the conditional node must be created as well as another stream. partitioning technique will build a stream to communicate this condition from its definition site to the cond-Φ node's task. We build on the concept of treegion, a single-entry multiple-exit control-flow region induced by a sub-tree of the control dependence graph. In the following, we assume the control flow is structured, which guarantees that the control dependence graph forms a tree. Every sub-tree can be partitioned into concurrent tasks according to the control dependences originating from its root. Any data dependence connecting a pair of such tasks induces communication over a dedicated stream. We call taskM_N the N-th task at level M of the control flow tree.

S1 while (p != NULL) { S2 x = p->value; S3 if(c1) { S4 x = p-
S1 while (p1 = Φ loop (p0,p2)) { S2 x1 = p1->value; S3 if(c1) { S4 x2 = p1->value/2; S5 ip1 = p1->inner_loop; L1 while (ip2 = Φ loop (ip1, ip3)) { L2 do_something(ip2); L3 ip3 = ip2->next; } } x3 = Φ cond c1 (x1, x2); S6 ... = x3; S7 p2 = p1->next; } Figure 5. Uncounted nested loop in SSA form. //task0-0(main task) S1 while (p1 = Φ loop (p0, p2)) { //persistent-task1-1 #pragma task firstprivate (p1) output(x1) { S2 x1 = p1->value; } //persistent-task1-2 #pragma task firstprivate (p1) output(c1, x2) { S3 if(c1) { //persistent-task2-1 #pragma task firstprivate (p1) output(ip1) lastprivate(x2) { S4 x2 = p1->value/2; S5 ip1 = p1->inner_loop; } //persistent-task2-2 #pragma task input(ip1) { L1 while (ip2 = Φ loop (ip1, ip3)) { //parallel -task3-1 #pragma omp task firstprivate (ip2) { L2 do_something(ip2); } L3 ip3 = ip2->next; } } } } //persistent-task1-3 #pragma task input(c1, x1, x2) { x3 = Φ cond c1 (x1 ,
In Figure 5, after building the control dependence tree, one may partition it into 3 tasks (task1_1, task1_2 and task1_3) at the root level, and for task1_2, one may further partition this task into inner nested tasks task2_1 and task2_2. One may then check for data parallelism in the inner loops; if they do not carry any dependence, one may isolate them in additional data-parallel tasks, such as task3_1 in this example.

Figure 6 shows the task and stream-annotated code using an OpenMP syntax. Figure 7 shows the nested pipelining and data parallelization corresponding to the partitioned code. The main task will be executed first, and a pipeline will be created for the main task and its inner tasks three task1_1, task1_2 and task1_3. Among these, the same variable x used to be defined in the control flow regions of both task1_1 and task1_2, to be used in task1_3. This output dependence must be eliminated prior to partitioning into tasks, so that task1_1 and task1_2 could be decoupled, while task1_3 may decide which value to use internally.

Nested tasks are introduced to provide fine grained parallelism. It is of course possible to adapt the partition and the number of nesting levels according to the load balancing and synchronization overhead. The generated code will be well structured, and simple top-down heuristics can be used.

In the execution model of OpenMP 3.0, a task instance is created whenever the execution flow of a thread encounters a task construct; no ordering of tasks can be assumed. Such an execution model is well suited for unbalanced loads, but the overhead of creating tasks is significantly more expensive than synchronizing persistent tasks. To improve performance, we use the persistent task model for pipelining, in which a single instance will handle the full iteration space, consuming data on the input stream and producing on the output stream [START_REF] Pop | A stream-comptuting extension to OpenMP[END_REF]. In Figure 7, all the tasks except task3_1 use the persistent model to reduce the overhead of task creation; task3_1 is an ordinary task following the execution model of OpenMP 3.0 (instances will be spawned every time the control flow encounters the task directive). All these tasks will be scheduled by the OpenMP runtime.

One problem with the partitioning algorithms is the fact that the def-use edges (scalar dependences) can become very large, sometimes quadratic with respect to the number of nodes [START_REF] Kennedy | Optimizing compilers for modern architectures: a dependence-based approach[END_REF]. Figure 8 (left) presents an example that illustrates this problem, Statements S1, S2 define the variable x. These definitions all reach the uses in the statements S3, S4 by passing through S5. Because each definition could reach every use, the number of definition-use edges is proportional to the square of the number of statements. These dependences constitute the majority of the edges in a PDG. SSA provide a solution to this problem. In SSA form, each assignment creates a different variable name and at point where control flow joins, a special operation is inserted to merge different incarnations of the same variable. The merge nodes are inserted just at the place where control flow joins. Figure 8 (right) is the original program under SSA form. A merge node (Φ) is inserted at S5, and killed the definition of S1 and S2. We could see here, in the SSA form, we could reduce the definition-use edges from quadratic to linear. The systematic elimination of output dependences is also facilitated by the SSA form, with a Φ node in task3_1. Notice that the conditional expression from which this Φ node selects one or another input also needs to be communicated through a data stream.

When modifying loop distribution to rely on tasks and pipelining rather than barriers, it is not necessary to distribute the loop control node and one may run it all in the master task, which in turn will activate tasks for the inner partitions. The statements inside each partition form a treegion whose root is the statement that is dependent on the loop control node. With pipelining inserted, distributed loops could be connected with pipelining when there are data dependences.

One concern here is that loop distribution with task pipelines may not provide expressiveness to extract pipeline parallelism. This is not a problem however, since we may apply the same method to every conditional statement rooted treegion, with some special care to the nested tasks, we could get fine grained parallelism without explicitly decoupling the control dependences. Considering again the example in Figure 4, its control dependence tree is given in Figure 9. The root treegion includes all the nodes in the control dependence graph, treegion1_2 represents the treegion at conditional level 1 and its root is node 2, treegion1_3 is at conditional level 1 and includes nodes (S3,S4,S5,L1,L2,L3). treegion2_1 is in conditional level 2 and its root is node (L1), which is the loop control node of the inner loop.

So following our approach, we may start from the treegion at conditional level 0, which is the whole loop, an implicit task will be created as the master task. For the treegions at level 1, we could create them as sub-tasks running at the context of the main task. If there are data dependences between the treegions at the same level and without recurrence, we will connect them with communication streams. If there is a dependence from the master task to one inner task, the value from the enclosing context can be forwarded to the inner task like in a firstprivate clause of OpenMP. Dependences from an inner task to the master task are also supported, although lastprivate is not natively supported for OpenMP3.0 tasks, it is a necessary component of our streaming task representation. lastprivate(x) is associated with a synchronization point at the end of the task and makes the value of x available to the enclosing context. The same algorithms could be recursively applied to the treegion at the next inner level. e.g. For treegion1_3 at level 1, the sub treegion at level 2 is treegion2_4, treegion2_5 and treegion2_1, we could create sub-tasks by merging treegion2_4 and treegion2_5 as one sub-task and treegion2_1 (which is also the inner loop) as one sub-task, or just for part of them. To reveal data parallelism, we can reuse the typed fusion algorithm introduced by McKinley and Kennedy [START_REF] Kennedy | Typed fusion with applications to parallel and sequential code generation[END_REF]: it is possible to fuse communicating data-parallel nodes to increase the synchronization grain or improve the load balancing. In this example, the loop in node L2 does not carry any dependence, and we need to decouple it from its enclosing task to expose data-parallelism.

Partitioning Algorithm

In this section, we present our partitioning algorithm, based on the SSA and treegion representations. We define our model and the important constructs that will be used by our algorithm, then we present and describe our algorithm.

Definitions

In this work, we are only targeting natural structured loops [START_REF] Böhm | Flow diagrams, turing machines and languages with only two formation rules[END_REF]. Such loops are single-entry single-exit CFG sub-graphs with one entry block and possibly several back edges leading to the header from inside of the loop. break and continue statements can be preprocessed to comply with this restriction, but we plan to lift it altogether in the future.

Treegion The canonical definition of a treegion is a non-linear, single-entry multiple-exit region of code containing basic blocks that constitute a sub-graph of the CFG. We alter this definition to bear on the Control-Dependence Graph (CDG) instead, so we will be looking at single-entry multiple-exit sub-graphs of the CDG.

Loop Control Node

In the representation we employ later, we will use the loop control node to represent the loop. The loop control node include statements which will evaluate the loop control expression and determines the next iteration.

Although control dependences in loops can be handled by the standard algorithm by converting them to a control flow graph, there are advantages in treating them as a special case with coalescing them in a single node (loop control node): not only the backward dependence is removed by building the loop control node so that the control dependence graph will form a tree, but also, this node can be used to represent the loop in all sort of transformations.

Conditional Level

The control dependence graph of the structured code is a tree after building the loop control node. The root of the tree is the loop control node at the loop's outermost level. We define the conditional level for every node in the control dependence graph as the depth of the node in the tree. The root of the tree with depth 0 has conditional level 0.

We define the conditional level for the treegion is the conditional level of the root node of the treegion (subtree). We define treegionN_M to identify a treegion where N is the conditional level of the treegion and M is the root node number of the treegion.

The algorithm

The algorithm takes an SSA representation of a single function, and returns a concurrent representation annotated with tasks and communication streams.

Step 1: Transform Conditional Statements to Conditional Variables To achieve fine-grained pipelining, conditional statements are split to conditional variables. As showed in Figure 10. Full conversion to three-address SSA form is also possible (as it is performed in GCC or LLVM, for example).

if (condition(i)) //is transformed to c1 = condition(i) if (c1)
Figure 10. Split conditional statements to expose finer grained pipelining.

Step 2: Build the Program Dependence Graph under SSA By building the program dependence graph, the control dependence graph, data dependence graph (through memory) and scalar dependence graph (through registers) are built together.

The control dependence graph for the structured code is a tree, the root of the tree is the loop control node. The leaves of the tree are non-conditional statements and the other nodes inside the tree are the conditional statements or the loop control node of the inner loops. We start from building the control dependence graph, and evaluate the conditional level for each node in the graph. Every node inside the control dependence graph is an statement from the compiler's intermediate representation of the loop except for the loop control node. The loop control node will be built by searching the strongly connect component started from the loop header node (at each loop nest level) in the program dependence graph.

The data dependence graph could be built by the array dependence analysis [START_REF] Kennedy | Optimizing compilers for modern architectures: a dependence-based approach[END_REF] for the loop. We should analyse every pair of data dependences to mark the irreducible edges in a later step if there are recurrence.

Step 3: Marking the Irreducible Edges A partition can preserve all dependences if and only if there exists no dependence cycle spanning more than one output loop [START_REF] Allen | Automatic translation of fortran programs to vector form[END_REF][START_REF] Kuck | Dependence graphs and compiler optimizations[END_REF]. In our case, for the treegion at the same conditional level, if there are dependences that form a cycle, we mark the edges in between as irreducible. If we have statements in different conditional level, we promote the inner one to its ancestor until both of them are in the same treegion, mark the promoted root node and the other root node as irreducible. The algorithms is presented in Figure 11.

Step 4: Structured Typed Fusion Before partitioning, to reveal data parallelism, we type every node in the dependences graph as parallel or !parallel. If there are loop-carried dependence inside this node, then it should be typed as !parallel, otherwise, typed as parallel.

The parallel type nodes are candidates for data parallelization. The goal is to merge this type of nodes to create the largest parallel loop, reducing synchronization overhead and (generally) improving data locality. Further partitioning can happen in the following step, starting from this maximally type-fused configuration. Given a DAG with edges representing dependences and the vertices representing statements in the loop body, we want to produce an equivalent program with minimal number of parallel loops. We want it to be as large as possible to balance the synchronization

SSA representation

We are using the Static Single Assignment (SSA) form as an intermediate representation for the source code. A program in SSA form if every variable used in the program appears a single time in the left hand side of an assignment. We are using the SSA form to eliminate the output dependences in the code, and to disambiguate the flow of data across tasks over multiple producer configurations. Considering the example in Figure 15, if we partition the statements into (S1), (S2,S3), (S4), we need to implement precedence constraints for the output dependence between partition (S1) and (S2,S3), which decreases the degree of parallelism and induces synchronization overhead.

Eliminating the output dependences with the SSA form leads to the introduction of multiple streams in the partitioned code. In order to merge the information coming from different control flow branches, a Φ node is introduced in the SSA form. The Φ function is not normally implemented directly, after the optimizations are completed the SSA representation will be transformed back to ordinary one with additional copies inserted at incoming edges of (some) Φ functions. We need to handle the case where multiple producers in a given partition reach a single consumer in a different partition. When decoupling a dependence whose sink is a Φ node, the exact conditional control flow leading to the Φ node is not accessible for the out-of-SSA algorithm to generate ordinary code.

Task-closed Φ node In SSA loop optimization, there is a concept called loop-closed Φ node, which implements the additional property that no SSA name is used outside of loop where it is defined. When enforcing this property, Φ nodes must be inserted at the loop exit node to catch the variables that will be used outside of the loop.

Here we give a similar definition for task-closed Φ node: if multiple SSA variables are defined in one partition and used in another, a phi node will be created at the end of the partition for this variable. This is the place where we join/split the stream. We need to make sure that different definitions of the variable will be merged in this partition before it continues to a downstream one. This node will be removed when converting back from SSA.

Task-closed stream Our partitioning algorithms generate nested pipelining code to guarantee that all communications follow the synchronous hypothesis. For each boundary, if there are one or more definitions of a variable coming through from different partitions, we insert a consumer at this boundary to merge the incoming data, and immediately insert a producer to forward the merged data at the rate of the downstream control flow.

1. When partitioning from a boundary, if inside the treegion, there are multiple definitions of a scalar and it will be used in other treegions which has the same conditional level, we create a Φ node at the end of this partition to merge all the definitions, and also update the SSA variable in later partitions.

2. If there is a Φ node at the end of a partition, insert a stream named with the left-hand side variable of the Φ node.

3. At the place where this variable is used, which is also a Φ node, add a special stream-Φ node to consume.

4.

To generate code for the stream-Φ, use the boolean condition associated with the conditional phi node it originates from.

Let us consider the SSA-form example in Figure 15 where we partition the code into (S1,S2,S3) and (S4,S5). A Φ node will be inserted at the end of the first partition, r1 4 = phi(r1 1, r1 2), the Φ node in a later partition should be updated from r1 3 = Φ(r1 1, r1 2) to r1 5 = Φ(r1 4). In the second step, we find out that in partition (S1,S2,S3), there is a Φ node at the end, so we insert a stream to produce there. And in partition (S4,S5), after the Φ node there is a use of the variable, so we insert a stream consume. The generated code will look like Figure 16. This example illustrates the generality of our method and shows how fine-grain pipelines can be built in presence of complex, multilevel control flow.

If we decide to partition the statements into (S1), (S2,S3), (S4,S5), which is the case for multiple producers, the generated code will look like in Figure 17. For multiple consumers, the stream extension of OpenMP will broadcast to its consumers, which is appropriate for our case.

Conclusion

In this paper, we propose a method to decouple independent tasks in serial programs, to extract scalable pipelining and data-parallelism. Our method leverages a recent proposition of a stream-processing extension of OpenMP, with a persistent task semantics to eliminate the overhead of scheduling task instances each time a pair of tasks need to communicate. Our method is inspired by the synchronous hypothesis: communicating concurrent tasks share the same control flow. This hypothesis simplifies the coordination of communicating tasks over nested levels of parallelism. Synchrony also facilitates the definition of generalized, structured typed fusion and partition algorithms preserving the loop structure information. These algorithms have been proven to be essential to the adaptation of the grain of parallelism to the target and to the effectiveness of compile-time load balancing. These partitioning algorithms also handle DOALL parallelization inside a task pipeline. We are using a combination of SSA, control dependence tree and (nonscalar) dependence graph as an IR. With the support of SSA, our method eliminates the nested multiple producer and multiple consumer problems of PS-DSWP. SSA also provides additional applicability, elegance and complexity benefits. This work is currently under development in a development branch of GCC, the partitioning algorithms is partially developed. For the code generation part, we first need to migrate the existing OpenMP expansion pass of GCC to work under SSA form, which has been a long-running challenge. When this work is complete, our method will leverage the array data-flow analysis of the Graphite polyhedral compilation pass of GCC to provide more precise data dependence information in loop nests with regular control flow.

Figure 1 .

 1 Figure 1. Barriers inserted after loop distribution.

Figure 2 .

 2 Figure 2. Syntax for input and output clauses.

Figure 4 .

 4 Figure 4. Uncounted nested loop before partitioning.

Figure 6 .

 6 Figure 6. Loops after partitioning and annotated with OpenMP stream extension.

Figure 5 Figure 7 .

 57 Figure5shows the conversion to SSA form. Just like GCC, we use a loop-closed SSA form distinguishing between loop-Φ and cond-Φ nodes. The latter take an additional condition argument, appearing as a subscript, to explicit the selection condition. The

Figure 8 .

 8 Figure 8. Definition and use edges in the presence of control flow.

Figure 9 .

 9 Figure 9. Control dependence graph of Figure 4. Express the definition of treegion.

/Figure 15 .

 15 Figure 15. Normal form of code (left) and SSA form of the code (right).

Figure 16 .

 16 Figure 16. Apply our algorithm to generate the parallel code. Producer thread (left) and consumer thread (right).

Figure 17 .

 17 Figure17. Multiple producers with applied our algorithm, the generated code.

Acknowledgments This work was partly funded by the European FP7 project TERAFLUX id. 249013, http://www.teraflux.eu

Workshop on Intermediate Representations

Chamonix, 2011

overhead. Even when we don't want that coarse grained parallel loops, we could also partition between iterations if possible.

In our case, we need a structured typed loop fusion algorithm. We revisit McKinley and Kennedy's fast typed fusion [START_REF] Kennedy | Typed fusion with applications to parallel and sequential code generation[END_REF] into a recursive algorithm traversing the control dependence tree. Starting from the treegion at conditional level 0, which is the whole loop, we will check if there are loop carried dependences between iterations. If there are no loop carried dependence, we stop here by annotating the whole loop as parallel. If there are, we are going into each inner treegion, identifying those that have no loop carried dependences. If some of them carried no loop carried dependence, mark the nodes as parallel and try to merge them. There are some constraints when we fuse the nodes: (1) parallelization-inhibiting constraints;

(2) ordering constraints. The parallelization-inhibiting fusion is that there are no loop-carried dependences before fusion, but will have the loop carried dependences after. So we should skip this kind of fusion which will degrades data parallelism. The ordering constraints describe that two loops cannot be validly fused if there exists a path of loop-independent dependences between them that contains a loop or statement that is not being fused with them.

The time complexity of the typed fusion algorithms is O(E+V) [START_REF] Kennedy | Typed fusion with applications to parallel and sequential code generation[END_REF], and our structured extension has the same complexity. Step 5: Structured Partitioning Algorithms Updating the CDG after typed fusion, start from the treegion which has conditional level 0 for our partitioning algorithms, and for all of its child treegions at conditional level 1, we should decide where to partition. The partition point could be any point between each of these treegions at the same level except the irreducible edges that we have created in step 3. The algorithm may decide at every step if it is desirable to further partition any given task into several sub-tasks.

Look at the example Figure 13: The code in Figure 13 (left) is partitioned into 2 tasks, and one task (task1_2) is partitioned further into 3 sub-tasks.

Code Generation

After the partitioning algorithms, we have decided the partition point between the original treegions, with the support of the stream extension of OpenMP. We ought to generate the code by inserting the input output directives. With the support of nested tasks, relying on the downstream, extended OpenMP compilation algorithm (called OpenMP expansion). But some challenges remain, especially in presence of multiple producers and consumers. We are using SSA form as an intermediate representation and generating the streaming code.

Decoupling dependences across tasks belonging to different treegions

Clearly if we decouple a dependence between tasks in the same treegion, the appropriate input and output clauses can be naturally inserted. But what about the communication between tasks at different level? Considering the example in Figure 14, if we decide to partition the loop to 3 main tasks: task1 1 with S1, task1 2 with (S2,S3), and task1 3 with S4, task1 2 is further divided to task2 1 with S3. If we insert the produce and consume directly into the loop, unmatched production and consumption will result. The answer comes from following the synchronous hypothesis and slightly modifying the construction of the SSA form in presence of concurrent streaming tasks.