Damián Vicino
email: damian.vicino@carleton.ca

Olivier Dalle
email: olivier.dalle@unice.fr

Gabriel Wainer
email: gwainer@sce.carleton.ca

A data type for discretized time representation in DEVS

Keywords: I.6.2 [SIMULATION AND MODELING]: Simulation Languages Simulation, Data Type, Time, DEVS

This paper addresses the problems related to data types used for time representation in DEVS, a formalism for the specification and simulation of discrete-event systems. When evaluating a DEVS simulation model into an actual computer simulation program, a data type is required to hold the virtual time of the simulation and the time elapsed in the model of the simulated system. We review the commonly data types used, and discuss the problems that each of them induce. In the case of floating point we show how, under certain conditions, the simulation can break causality relations, treat simultaneous events as non simultaneous or treat non simultaneous events as simultaneous. In the case of integers using fixed unit we list a number of problems arising when composing models operating at different timescales. In the case of structures that combine several fields, we show that, at the cost of a lower performance, most of the previous problems can be avoided, although not totally. Finally, we describe an alternative representation data type we developed to cope with the data type problems.

INTRODUCTION

Discrete-event simulation is a technique in which the simulation engine plays an history following a chronology of events. The technique is called "discrete-event" because the processing of each event of the chronology takes place at discrete points of a time-line and takes no time with respect to virtual simulated time (even though it may take a nonnegligible time to compute with respect to the real wall-clock time). Unless the system being modeled already follows a discretized schedule (e. g. basic traffic lights), the process of modeling a phenomenon using a discrete-event simulator can be considered as a discretization process.

We will discuss some implementation issues related to the representation of the time in a discrete-event simulation engine. Since this problem is closely dependent on the implementation details, we choose to restrict our discussion to the particular case of the DEVS discrete-event simulation engines [START_REF] Chow | Revised DEVS: a parallel, hierarchical, modular modeling formalism[END_REF] [START_REF] Zeigler | Theory of modeling and simulation[END_REF]. The DEVS formalism provides a theoretical framework to think about Modeling using a hierarchical, modular approach.

The modeling hierarchy has two kinds of components: basic models and the coupled models. The basic models (with ports) are defined as a tuple: A =< S, X, Y, δint, δext, λ, ta > where: S is the set of states, X = {(p, v)|p ∈ IP orts, v ∈ Xp} is the set of input ports and values, Y = {(p, v)|p ∈ OP orts, v ∈ Yp} is the set of output ports and values, δint : S → S is the internal transition function, δext : Q × X → S is the external transition function, Q = {(s, e)|s ∈ S, 0 ≤ e ≤ ta(s)} is the total state set (with e is the time elapsed since last transition), λ : S → Y is the output function and ta : S → R + 0,∞ is the time advance function.

At any time of the simulation the basic model is in a specific state awaiting to complete the lifespan delay returned by the ta function, unless an input of a new external event occurs.

If no external event is received during the lifespan delay, the output function λ is called, and the state is changed according to the value returned by the δint function. If an external event is received, then the state is changed to the value returned by the δext function, but no output is generated.

Coupled models define a network structure in which nodes are any DEVS models (coupled or basic) and links represent the routing of events between outputs and inputs or to/from upper level. Formally, the Coupled Models are also represented by a tuple: C =< X, Y, D, M, EIC, EOC, IC, SELECT > where: X is the set of input events, Y is the set of output events, D is an index for the components of the coupled model, M = {M d |d ∈ D} is a tuple of basic models as previously defined, EIC, EOC, and IC are respectively the External Input Coupling, External Output Couplings and Internal Couplings that explicit the connections and port associations respectively from external inputs to internal inputs, from internal outputs to external outputs, and from internal outputs to internal inputs, SELECT : 2 D -{} → D is the tie-breaker function that sets priority in case of simultaneous events.

Zeigler [START_REF] Zeigler | Theory of modeling and simulation[END_REF] also described the algorithm of a simple abstract simulator, able to execute the DEVS models. When it executes a model, this simulation algorithm exchanges only two types of data with the models: events and time.

In the last four decades, many DEVS simulation engines have been implemented and studied using different computation models, such as Sequential, Parallel, Distributed, Cloud, etc. All these DEVS engines handle the passing of events between models following the principles of the abstract simulator: Compute the time of next events (minimum of the times remaining in each model until the next internal transition), compute resulting event outputs, route events in the network to their destination inputs, and recompute the time remaining until the next round (subtracting the time elapsed from the scheduled time advance).

Notice that time is never exchanged directly between models, but only between the models and the simulation engine. Therefore, a common time representation is needed between the models and the simulation algorithm, but this time representation can possibly vary from one model to another, provided that the simulation algorithm is able to make the translation.

Simultaneous events correspond to the cases in which the discretization that results from the discrete-event modeling produces identical time values. Simultaneous events occur in DEVS in two ways: either a model receives input events from multiple models at the same time, or it receives events from other models at the same time as it is reaching the end of the time advance delay. Some means are provided in each DEVS variants to ensure deterministic processing of simultaneous events. For example, in classic DEVS, a tie-breaker function named SELECT is used to resolve this situation.

The SELECT function is hard to define, and having an error in its definition can lead to unpredictably large modeling errors, as will be discussed in Section 3. However, assuming the SELECT function used is valid, we may still experience simulation errors due to the way the time is stored internally by the simulation engine. Indeed, in the DEVS formalism, the time is a real value, but it is well known that it is not possible to find a computer representation for all real numbers. Therefore, an approximation is needed for the values that have no exact representation. This approximation is the result of a quantization process. Unfortunately, contrary to modeling errors, simulation errors due to the quantization of time are hidden to the modeler even though they may have similar impact on the simulation results.

The focus of this paper is on this particular problem of selecting an internal representation and the possible effects of such a selection on the simulation results.

BACKGROUND

To contextualize the problem, we reviewed the source code and documentation of the eight following DEVS simulator projects:

• ADEVS [START_REF] Nutaro | ADEVS project[END_REF]: "A Discrete EVent System simulator" is developed at ORNL in C++ since 2001; • CD++ [START_REF] Lòpez | Improved Cell-DEVS model definition in CD++[END_REF]: "Cell-DEVS++" is developed at Carleton University, in C++, since 1998; • DEVSJava [START_REF] Zeigler | Introduction to DEVS modeling & simulation with JAVATM: Developing component-based simulation models[END_REF]: "DEVJava" is developed at ACIMS in Java since 1997; it is based on its predecessor DEVS++ [START_REF] Hwang | DEVS++: C++ Open Source Library of DEVS Formalism[END_REF] (implemented in C++). Currently it is part of the DEVS-Suite project [START_REF] Hu | The Architecture of GenDevs: Distributed Simulation in DEVSJAVA[END_REF] A common feature in all of them is to embed the passive state representation in the time representation data type using some infinity representation.

All simulators, but ADEVS and CD++, represent Time using the double precision floating point data type1 provided by the programming language. All but one use the internal representation of infinity provided by the data type. The exception is ODEVSPP that uses MAX DOUBLE constant in place of internal infinity.

In the case of ADEVS the time is represented using a C++ class template. The class works as a type wrapper adding infinity representation to the provided type. The operators are defined in a way that, if one of the parameters is an infinity the defined algorithms are used, else the wrapped type operators are used.

In the case of CD++ the time is represented using a class with five attributes: four integers for the hours, minutes, seconds and milliseconds, and a double precision floating point for sub-milliseconds timings. There is not an explicit representation of infinity, but the passive state is evaluated comparing the Hour field to 32767.

Something to keep in mind is that all DEVS simulators execute a single Model (sometimes composed of several others). This model is called the top or root model.Regardless of the simulator, a global variable is needed to keep track of the global Simulation Time. Some simulator may try to reduce the scope of the synchronization in sub-models wherever it is possible, while others simply use the root model in a systematic way to centralize all timing information. This is important to remember when evaluating Time representation range.

TIME REPRESENTATION PROBLEMS

In DEVS modeling, the domain of the time variables is R + . This choice results in a quantization problem at the time of implementing the model. Depending on the data type chosen for implementation, different approximations or restrictions may be observed.

Approximation effects belong to two categories: time shifting, and event reordering. In the case of time shifting, we have an error on the precision of the computation, but the time line of events is not affected. In the second case, the approximation leads to a different time line of events than the one expected in theory using exact arithmetic on real numbers, which may, in turn, break the causality chain. In case the causality chain is broken, the error can be arbitrarily large, such that no error bound can be found in the general case.

A characteristic of reordering errors is that they cannot be explicit in formal proofs on the paper. Furthermore, reordering errors are practically impossible to predict, and may occur following irregular frequency patterns, which, in the worst case, can fool the validation tests.

The following example shows that even the simplest model may be subject to reordering errors. In this example, we model a counter with 2 inputs, one to increase the count and the other to reset and output the current counter result.

The root coupled model is:

C =< X, Y, D, M, I, Z, SELECT > where: The output of this simulation is, in theory, the value 10 on the "out" port of C every 1 second.

X = ∅, Y ⊆ {out} × N, D = {1,
We implemented this in DEVSJAVA 3.1 and pyDEVS 1.1, and kept it running for a while. In DEVSJava we got a different result than expected: the output consisted of a majority of 10 values, but also some 9 and 11 values. In pyDEVS we first obtained the expected results for this experiment. Nevertheless, after changing the parameter of M2 from 1.0 to 100.0 seconds, we observed similar discrepancies, with a few occurrences of 999 and 1001 values on C's "out" port instead of only 1000 as theoretically expected. Interestingly, the errors of this model occur following a regular pattern (a geometric law of factor 4 in both simulators).

This result can be explained knowing that both simulators use single or double precision Floating Point numbers for time representation. Indeed, the value 0.1 is well-known bad quantization point in the Floating Point standard. Therefore, every time an algorithm accumulates many times the value 0.1, as the model M1 does, the rounding error accumulates, which leads to the observed errors. We discuss in more details the case of Floating Point numbers in section 3.1.

Looking at other data types, such as Integers for example, does not fully solve the problem. We further discuss the case of integers in section 3.2.

In Section 2, we also reviewed implementations that use a structured or object-oriented data-type. However, as further explained in section 3.3, these data-type don't have better results than using floats or integers because they work mostly as wrappers of these types or small extensions to their range.

Floating point data types

Adopted by a majority of simulators, this data type leads to several variants due to the selection of different precision levels (single, double, or extended) and how the passive state of a model is represented (the semantics of which being to wait for an infinite time). For example, infinity can be mapped onto a reserved value or an additional variable, e. g. a boolean, can be use to handle the special case of infinity, in a structure or wrapper data type.

Strengths of floating point

The floating point data type was engineered to represent an approximated real number and to support a wide range of values. The basic structure is the use of a fixed length mantissa and a fixed length exponent.

The main strengths of using floating points are:

• A compact representation, usually between 32 and 128 bits; • Implemented in almost every processor; • A large spectrum between max and min representable numbers; • An internal representation for infinities;

• A widely used standard, making the simulator code more portable; • Largely studied mechanics of its arithmetic approximations.

On the other hand, floating point has well known limitations [START_REF] Goldberg | What every computer scientist should know about floating-point arithmetic[END_REF], but are still often considered as an acceptable tradeoff, due to their intuitive use. We review those problems and limitations in the following sections, in the particular case of simulation. It is important for the analysis to notice that the only operations used internally by the DEVS Simulator on time values are additions, subtractions and comparisons.

Quantization errors

The floating point arithmetic rounds the results continuously, not only with divisions (as it happen with integer), but also with additions and subtractions. Rounding errors result in two kind of errors: shifted values, and artificially coincidental values.

We don't elaborate more on the case of shifted values, as this case was already discussed with the example in the previous section, that showed this kind of error can easily break the causality chain, e. g.. by accumulation.

The rounding of values due to quantization may also produce artificially coincidental values. As mentioned in the introduction, the case of simultaneous events must be dealt with care, which explains the need for a SELECT function. Unfortunately, when the coincidence is the result of a hidden process, the SELECT function may not be able to make a proper decision. Indeed without the rounding error, two distinct time values have an implicit natural order, and this natural order may not be the one produced by the SELECT function.

Non fixed step

The floating point representation uses a variable step size between consecutive represented values. The reasons for having variable step come from the nature of the incremental sequence of powers: As the absolute value of the exponent value increases, more bits are needed for coding the exponent value, and these bit are taken from the least significant side of the mantissa.

Some assertion code can be provided to detect when an addition of non zero numbers is done and to check if the result is equal to one of the operands.

Cancellation

This kind of problem happens when subtracting 2 close numbers. The effect is that the mantissa loses precision. For instance, when subtracting 0.999 from 1.000, we get 0.001:

The mantissa digits are shifted to obtain 1.0E-3, and the result is the loss of the 3 significant digits in the mantissa. This cancellation effect makes the rounding problem even worse.

A solution sometimes proposed to solve some of the floating point issues, since the time always advances, is to use a framing mechanism. This is a way to keep the computations in the dense area around the origin. This solution is implemented using a threshold value t0: Once the models have reached this value, t0 is subtracted from all the clocks in the system.

A downside of this solution is that it increases the number of subtractions, which in turn may lead to an increase in Cancellation errors.

Error detection can be implemented by keeping track in every operation of how the mantissa significance was affected. This does not solve the problem, but may still give a confidence indicator as the assertion in the non fixed step one.

Compiler optimization dangers

In 1991 Goldberg [START_REF] Goldberg | What every computer scientist should know about floating-point arithmetic[END_REF] shows a set of optimizations that when used with floating points can affect the precision of the computation. Changing compilers (for updating or portability) can change results if any problematic optimization is added or removed, which might be hard or impossible to detect.Also, at the time of distributing the simulation in an heterogeneous cloud or grid, different architectures, operating systems and compilers may be used making the problems scale together with the infrastructure.

Standard complexity

If the floating implementation follows the IEEE754, which is common, there are more issues to consider. IEEE754 gives special meaning to some quantities (NAN, Infinities, positive and negative zeros); it defines exceptions, flags, trap handlers and the option to choose a rounding models. These mechanisms need to be taken into account when implementing the Simulator, or they could result in unexpected behaviors.

Integer data types

Using integers is the first idea that comes up when trying to avoid the issues related to rounding and precision with floating point numbers. For example, in the field of networking, well-known simulators such as NS and Omnet++ have gone through major rewriting in order to change their representation from floating point to integers: NS3 has chosen a 128 bits integer representation, while Omnet++ v4.x has chosen a 64 bits integer representation.

Strengths of integers

The Integer type is used to represent a subset of consecutive numbers in Z.

Some interesting characteristics of the Integer type are:

• Its compact representation, usually between 8 and 128 bits; • Its generalized support by all processors;

• Its fast and exact arithmetic, except for divisions;

• Its fixed step between any two consecutive values;

Using integers in place of floats can be seen as a trade-off in which the large range and the ability to approximate the (dense) real numbers is traded against the accuracy of an exact representation, but only for a limited range of the (sparse) natural numbers. In many situations, when the range is not that much a concern, integers offer a better trade-off.

Overall range

For a given number of bits, the absolute range of values supported by an integer type using that amount of bits is much narrower than the absolute range of a floating point using the same amount of bits. Indeed, starting with the single precision floating point (32 bits), the largest absolute value is around 3 × 10 38 while it is only around 2 × 10 9 for a 32 bits integer; At the other end, the smallest single precision floating point value is around 1 × 10 -38 while for any positive integer the smallest (absolute) value is 1.

Trying to cover the same range as a floating point with in-tegers, using the smallest step offered by the floating points, is feasible but very space-inefficient for single precision (76 bits vs. 32) and unreasonable for higher precision.

However, since there is no approximation with the additive arithmetic on integers, a frame shifting mechanism as described in previous section is safe to implement and will not incur problems such the cancellation mentioned about floating points.

Quantity of unused representation

Since the smaller representable value is 1, the unit associated to the variable is the smallest representable time of a model. Some simulator use a fixed unit value, e. g. nano-seconds. This fixed unit implies that all timings are aligned on multiples of one nanosecond. Therefore, a model of a generator that would output ticks every second, would waste 99.9% of an already limited representation range.

To avoid this waste, an analysis of the model can help to adjust the unit. In this case, the unit does not need to be a standard one: a rational, eg. 1 7 s, can even be chosen if it allows to capture all the possible timings needed for the simulation of the model time line.

Multiple scales and model composition

One of the strengths of the DEVS formalism is its ability to compose models. This favors the reuse of previously developed models to create more complex ones and the division of the system into smaller parts.

Assuming the time unit is hard-coded within the Simulator implementation requires a unit to be found that covers all possible models. This approach can be used for simulators dedicated to a specific purpose, eg. micro-controllers electronics or planets dynamics. However, if the simulator is intended to be multi-purpose, it becomes impossible to find a compact and efficient representation that covers the full spectrum of possible requirements.

An alternative is to use time units defined locally within each model. For instance, this approach allows to pick months as the time unit for building a gravitational model of planets dynamics, and pico-seconds for micro-controllers as suggested by each model analysis.

When we have different time scales in different (Sub-)models, we need to have extra computation to compose them as DEVS Coupled Models. In case the time-step in one representation is divisible by the one in the other, the smallest step can be used as the common step, but the time-values in the larger-step model have to be systematically adjusted, at the cost of an extra multiplication operation.

When using time-scales that don't divide each other, e. g.

Structures and containers

Complex data types that come with each programming language can also be used be used in place of the afore-mentioned ones, in order to avoid some of the problems we have surveyed. These data types include: classes, objects, structs, tuples, arrays, vectors and other containers.

However, using such data types often results in new tradeoffs, such as user-friendliness vs. memory consumption, or vs. performance, or vs. complexity of time arithmetic.

Developer friendliness

In some implementations of time, a structure with an integer fields to handle different units of a time value is used. This can be done using an array, vector, class or struct, with the main objective of expanding the range. This kind of structure is used in the Time class in CD++ (which has also an extra float field) in which the fields are aligned on standard time sub-units (days, hours, minutes, nanoseconds). The use of human-friendly units, even when convenient for the model definition, makes the arithmetic more complicated, and is space inefficient. For example, using a byte-long integer to represent minutes on 8 bits, wastes 76% of the values.

Memory consumption and access

Integers and floating point numbers can both use either 16, 32, 64 or 128 bits depending on the processor architecture. A struct having multiple integers as previously defined, when implemented as a Class, Object or Struct, requires a pointer to the structure and offsets to each component. This gets worst if we combine with previous problem and we use human friendly representation. For efficiency purpose, it is advisable to use integer that match the processor register size, hence at least 64 bits on computing architectures currently available. An structure for [HH:MM:SS] using such integers for each field would use 224 bits and represent less range than a single 64 bits integer coding the same range using the seconds unit. In some languages, like C++, this can be reduced choosing an 8 bit integer data type, anyway this is not a possibility in some other languages like PHP or Python.

In case a vector-type collection is used, the memory footprint is increased to define pointers needed for iteration or other general purpose use of the data type as memory allocators.

In some architectures with few available registers, accessing multiple integers generates a higher number of cache miss and results in significant performance degradation.

Memory Alignment

In most architectures, the use of smaller than bus size or misaligned data types will require extra processing to be read and write. For example, if the bus width is 32 bits, 2 consecutive integer values of 16 bits can be read using a single memory access, but extra operations are needed for the actual values to be unpacked prior operation. Similarly, having a 16 bits integer followed by a 32 bits one ends up either with the second one being misaligned and requiring two memory access each time it is read/write, or with the waste of 16 bits to realign the representation.

Operation encapsulation

If we want this data type to be reusable, the operations giving the semantic to the data type and its representation need to be encapsulated. The only data types providing the mechanisms to embed operations in a data type are classes and objects, and in most languages one depends on the other. This encapsulation has to be compared with the usability of collections. In case the data type is first designed as a class or object, adding the operators to it doesn't increase the memory access to its internal fields. On contrary, in case the data-type is designed as the encapsulation of a collection, using a class or object wrapper, a new indirection level is added, which may result in increasing the number of memory access, with possible added degradation due to cache miss.

A NEW TIME DATA TYPE

The standard way to represent the measure of a time lapse is the product between its magnitude and the time unit: t = m × u. Our proposed data-type aims as finding how to best represent both of these quantities, the magnitude and the unit.

In Section 3 we described a set of problems and limitation with the representation and performance of the different data types used in DEVS implementations. Looking for a more adequate data type for the time computation in DEVS, we describe in the following subsections a set of representation requirements, performance requirements and we propose a data type covering most of those requirements to avoid the previously exposed problems.

Representation requirements

From previous analysis, we can derive the following requirements:

• For the unit part of the representation, different units must be supported; • For the magnitude part, we need to address four concerns:

1. Avoid (hidden) quantization as much as possible; 2. Avoid additive arithmetic errors; 3. Provide a large enough range for any time representation; 4. Support seamlessly the coupling of models possibly operating using different units; • Optionally, we may want to have an internal representation for infinity. (Some simulators don't need it.)

Performance requirements

We expect the use of the new data type not to impact significantly the performances of the simulation. Time is advanced in every model, and every coordinator has to operate with it.

We should keep our representation compact and avoid wasting representation space.

Also, when trade-offs for performance are possible, we prefer to optimize the additive operations and comaparisons, which are the only ones used internally by the simulator.

Since time scaling operations that result from the compositiopn of models may be expensive, we want to avoid them as much as possible.

The proposed data type

In this subsection, we first describe a data type that covers the representation requirements and enumerate it limitations. Afterward, we describe how to, doing a few trade-offs, we can adapt it to a more performant implementation.

For the magnitude, we need a representation that allows us to represent large subset of numbers without rounding, and that all operations result in values into the same subset. A large subset of R covering this 2 conditions is Q, any of its values q can be represented as q = qn q d where qn, q d ∈ Z.

A representation of numbers in Q includes the numbers of the Integer type representation and the numbers of the Floating Point type representation (but infinity). In addition, it covers new cases, such as the numbers required for representing frequencies that are not powers of 2. I.e. two generators, one sending three events per second and one sending one per second can be modeled and have a simultaneous event every second.

Unfortunately, a simple implementation of this data type requires more space and operations time than previous ones. The space complexity comes from having several pairs of integers that represent the same value. The time complexity comes from the need to compute reductions (Greatest Common Divisor) in the representation in order to keep it as compact as possible.

Recall from earlier that we note t = m × u.

To ensure that the range, m, can be large enough, a dynamic size representation can be used for the numerator and denominator. In the following, we assume m is a rational number q defined by a duple of integers that are large enough.

For u, the unit, we can use an enumeration of units and define the scaling function to convert the operands of binary arithmetics operations to the smallest common unit before proceeding.

The previous data types covers all the representation requirements.

The restriction is that magnitude values had to be represented by a subset of Q delimited by the size used to represent the 2 integers used to represent q and the units need to be a finite predefined set with relations well defined between them.

Unfortunately, using this data in a real simulation would waste a lot of representation space, and induce a big computation overhead when coupling models havoing different units.

A few optimization can help reduce this overhead in the general case and limit the increase of complexity to specific scenarios that where not supported by previous data types.

A first change is to replace the unit enumeration representation by an integer scale factor e and a patron unit p, so that any unit can be thought as p × 2 e . Using the Second as patron unit may be reasonable for general purpose if the e is represented with a 32 bit integer.

The second change would be to use a custom quantum qu ∈ Q and a scale i ∈ Z to represent the time inside each model as: qu×i. This quantum is a divisor of all the possible times that will be used by the model, modelers usually think about it at the time of modeling. In the case the modeler doesn't set the quantum, if the model has regularity in its behavior, qu can be readjusted by the Simulator automatically (by finding the Greatest Common Divisor between last used qu and the one from current operation), until a stable state is reached. An example of a case where it never reaches stability is when the frequency of events is continuously increasing, eg. with ta = ( 1 2 ) n where n is the quantity of previous outputs.

Thanks to this quantum, all operations that apply to operands having the same quantum, which is expected to be a most common case, become Integer computations. In case operands have different quantum, a quantum adjustment is needed, which may be an expensive operation used few times. The case in which an adjustment is needed very frequently corresponds to a model that was not properly handled by previous representations.

We propose then to reduce the qu representation to two common integers as large as the architecture bus. This trade-off of keeping a small qu can increase significantly the optimization possibilities of the data type and keep it compact. This reduction is not as restrictive as it may look, when combining to the unit representation all float and integer space is covered and we still have another component i to scale the values.

If we look the whole representation (magnitude and unit), with current changes a time lapse is represented as: t = i × qu × 2 e × p. In 32 bits architecture, the qu takes 64 bits, the e 32 bits, p is implicit and i we can fix it to 32 bits too and have a representation of 128 bits which most the time the operations are just comparison of the first fields and add/substract/compare on the last one, all possible with 32 bits Integer operations.

To simplify the notation, we define the unit u as: u = qu × 2 e × p and the time lapse t as: t = i × u.

When working with binary operators, if both operands have the same unit, the computation only involves i part of each operand. The problem of having the i part limited to 32 bits is the coupling of models using different units, which makes necessary to have binary operators supporting different units in their operands. If the two unit differ in scale (different value of e), the computation applies only to the overlapping segments of the qu part, but we need to check for possible carry to expand the representation.

In case the quantums of each time value are different, a common quantum needs to be found and then i and e adjusted accordingly. Since this is a more expensive operation, a good selection by the modeler of the quantum of each model is important to quickly reach a stable state with a common quantum.

A way to agree on the value of qu can be added either to the Simulator logic so the models negotiate it at initialization, or it can be encapsulated in the Time data type itself if using an Object Orient Language.

The problem with the initial negotiation approach is that it extends the logic of the simulator and doesn't really solve the problem when the time unit can change during simulation (as the model ticking half as previous tick showed before). Anyway, if the trade-off is worth it between predictability on the performance and flexibility to allow the change of the quantum in simulation time, it can be done.

In the dynamic case, a simple implementation of the negotiation mechanism that agrees in a global qu common to all models participating can be implemented just checking qu local = qu global at the time of operate, if they are equal then just operate with no worry, in the case they are not equal a new qu gloabl that divides both is computed, the local i and e are adjusted and the new qu is assigned to both. All models defined with different q l ocal will do this until reach consensus. Obtaining consensus using this approach is lineal in the quantity of basic models participating in the simulation, anyway stabilization of the representation in every model can require square of the quantity of Basic Models adjustments to internal representation variables. This approach after reach consensus is as performant as the previous one, but it doesn't exclude the use of Models that can change its qu. If a model changes it qu a new consensus needs to be reached, and automatically will be done.

The header of a Time c++ class with a mechanism to agree the q encapsulated may look like this: This c++ class implementing the operations as described earlier in this section covers fairly the requirements.

Another representation can be done changing the i and e data types to vector of pair of integers where each position in the vector has a segment of i and an associated e. This way, when operating between numbers that don't have digits in common, the addition or subtraction can be done as just a concatenation. Also, this representation saves space when several 0s are in the middle after an operation, this are just not represented in the vector, any segment of 0s can be removed from the vector. About complexity of the operations in the alternative representation. For add/subtract operations once the qu is stable the operations need to do a check of the compatibility and then the operation as described before. If we simulating something that was possible to simulate using Integers, the overhead is 3 integer compares for equality, if it was not possible to simulate the worst case is linear to the size of the vector used to represent i. If the vectors don't have gaps and the e matches, complexity is equivalent to use BigInt, in the case there is any gap or the e is miss-aligned a displacement of all segments may be done in the worst case and this is a linear operation in the size of the i vector. For the compare operations worst case is linear to the size of the i vector too.

Communicating values

From previous analysis, we can think about how to communicate models in different simulator modules to handle interaction with web server based simulator between different technologies.

An XML, JSON, or any serialization can be provided based in the class header described before.

An message using xml, may look like this: <? xml version = " 1.0 " ? > < time type = " Time " > < unit numerator = " 58 " denominator = " 57 " patron = " Second " / > < magnitude offset = " -64 " > 123812937129837129 </ magnitude > </ time > The magnitude can be broke in pieces and send the offsets and content of the segments and avoid completing the gaps with 0s. Anyway, given the time it takes to parse the XML the difference may not be significant and more validation is needed to check the overlaps when parsing.

PERFORMANCE COMPARISON

In this section, we show the results of 4 experiments we realized to measure the performance impact of choosing the proposed Data Type compared against other possible Time data types implementation.

The experiments are implemented in C++11 compiled with clang 500.2.79. Each experiment was run without compiler optimization and with -O2 compiler optimization. Five data types are used for comparison:

• Float: We use the native float implementation.

• Integer: We use the native integer wrapped in a class to make it also provide infinity arithmetic. • Rational : We use the boost [START_REF] Demming | Introduction to the Boost C++ Libraries; Volume I-Foundations[END_REF] implementation of rational (1.55) wrapped in the same class as Integer to make it also provide infinity arithmetic. • DTime: A naive implementation of the proposed data type using global agreement with statics as described in previous section. • FDTime: A naive alternative implementation to the proposed one. This one doesn't intent to reach global agreement, in place the agreement is produced between operands in each operation.

Experiment: series y=1

This experiment implements: In this experiment all types but float provide exact results. The problem with float is that after accumulating enough, the exponent is increased and the increment starts being discarded.

In Fig. 1 all data types show similar behavior, the worst implementation appears to be the Rational, while both FD-Time and Time appear to be similar to just using an Integer with Infinity operations.

In Fig. 2 Integer and Rational are as good as float after being optimized, while FDTime and DTime get almost no improve from the compiler optimization.

Experiment: series y=1/x

This experiment implements:

N x=1 1 x
Computing this series using integers makes no sense and computing it using float is very imprecise, it rounds most of the values. In Fig. 3 we obtain very similar results in DTime, FDTime and Rational and a float that goes 100x times faster (but returns the wrong results).

In Fig. 4 after run compiler optimizations, we see that float flattens getting a great improve, but the other data types are mostly unaffected by the optimizations, just doing slightly better at the start.

Experiment: series 1/2 x

This experiment implements: This experiment returns all representable values in float until the exponent is increased. In numeric methods, it may be suggested to start the Sum from the smaller number to reduce the error. Anyway, in DEVS, we don't have a way to reorder the operations, so it is a valid use case doing it this way. In Fig. 5 we obtain similar results to previous experiment. In Fig. 6 something changes, this time the float curse doesn't flatten.

Experiment: series 1/p

This experiment implements: Here, all numbers but 1/2 are not representable by float, so the result is very imprecise, also the prime numbers force to fail the simplification of the fractions in the internal representation of the other data types. In both, Fig. 7 and8, we obtain similar result to those in the N x=1 1 x experiment.

CONCLUSIONS

We surveyed a set of Simulators used to run DEVS formalism simulations and we analyzed the implications of the data type selected in those to characterize the need of a new Time data type. From this characterization we described requirements for the new data type and a possible implementation which covers them at the expense of a slight increase in spacial and time complexity compared to its predecessors when running in the same scenarios, but increases significantly the accuracy and the range of the data type compared to previous ones. We compared naive implementations of the proposed Data Type against other Data Types used in the surveyed simulators and against other data types, the worst penalty obtained was in the order of 100X. We think the trade-off is worth it even when work need to be done in optimizing and compacting the data type yet. We will focus our future work in the complexity aspects, the implementation and evaluation of the data type in the CD++ simulator and then we will study how this work applies to other Simulation formalisms in discrete events.

1 7 s and 1 5 s

 15 , finding a common time-step requires to find a common denominator, and adds a systematic multiplication on both sides.

  typedef pair < int , int > fraction ; template < typename T > class Time { public : // c o n s t r u c t o r s Time () = delete ; Time ( const fraction q , const int exp , const int m ) ; // o p e r a t o r s Time & operator =( const Time &) ; const Time operator +( const Time &) ; Time & operator +=( const Time &) ; // n e g o t i a t i o n static fraction q_global ; o m p a r i s o n o p e r a t o r s bool operator ==( const Time & , const Time &) ; bool operator !=( const Time & , const Time &) ; bool operator <( const Time & , const Time &) ; bool operator >( const Time & , const Time &) ; bool operator <=( const Time & , const Time &) ; bool operator >=( const Time & , const Time &) ;

Figure 1 :

 1 Figure 1: Series y = 1 with no compiler optimization.

Figure 2 :

 2 Figure 2: Series y = 1 with -O2 compiler optimization.

Figure 3 :

 3 Figure 3: Series y = 1/x with no compiler optimization.

Figure 4 :

 4 Figure 4: Series y = 1/x with -O2 compiler optimization.

Figure 5 :

 5 Figure 5: Series y = 1/2 x with no compiler optimization.

Figure 6 :

 6 Figure 6: Series y = 1/2 x with -O2 compiler optimization.

N x=1 1 p

 1 i with pi the i-esieme prime number in N.

Figure 7 :

 7 Figure 7: Series y = 1/pi where pi is the i-esieme prime with no compiler optimization.

Figure 8 :

 8 Figure 8: Series y = 1/pi where pi is the i-esieme prime with -O2 compiler optimization.

In the case of PyDEVS, the dynamic type is forced to Floating Point by systematically using Python casting syntax.


Oleksandr Zinenko
email: oleksandr.zinenko@inria.fr

Stéphane Huot
email: stephane.huot@lri.fr

Cédric Bastoul
email: cedric.bastoul@unistra.fr

Keywords: parallel programming, polyhedral model, software visualization, direct manipulation interface

Parallel systems are now omnipresent and their effective use requires significant effort and expertise from software developers. Multitude of languages and libraries offer convenient ways to express parallelism, but fall short at helping programmers to find parallelism in existing programs. To address this issue, we introduce Clint, a direct manipulation tool aimed to ease both the extraction and the expression of parallelism. Clint builds on polyhedral representation of programs to convey dynamic behavior, to perform automatic data dependence analysis and to ensure code correctness. It can be used to rework and improve automatically generated optimizations and to make manual program transformation faster, safer and more efficient.

I. INTRODUCTION

The large scale adoption of parallel architectures in modern electronic devices urges the software industry to provide developers with tools for building efficient parallel applications. Despite tremendous effort, parallel program development remains challenging, especially when porting legacy sequential applications or targeting new parallel architectures. Several languages and libraries, such as Cilk, X10, Chapel or UPC, offer high-level parallelism expression constructs, but parallelism extraction is left under the responsibility of the developer.

Most existing approaches to automate parallelism extraction focus on compiler techniques for loop vectorization and thread-level parallelization and rely on an algebraic representation of programs called the polyhedral model [START_REF] Feautrier | Polyhedron model[END_REF]- [START_REF] Bastoul | Contributions to High-Level Program Optimization[END_REF]. It captures the dynamic behavior of a class of compute-intensive loops and allows to perform accurate data dependence analysis to ensure semantics preservation, and to restructure the input program to expose parallelism. However, polyhedral compilers are black-boxes that do not help users if the proposed restructuring does not solve their problem, e.g. has worse dynamic characteristics or does not expose parallelism due to the lack of semantics preservation hints. Semi-automatic tools based on directive scripts were proposed [START_REF] Kelly | A framework for unifying reordering transformations[END_REF], [START_REF] Girbal | Semi-automatic composition of loop transformations for deep parallelism and memory hierarchies[END_REF], but this method has seen little interest as writing such scripts remains a tedious task, requiring expertise in a complex underlying theory.

Although the polyhedral model has a direct graphical representation featuring its geometric nature, it has not been explored yet as a way to leverage the power of the model through interactive visualizations. We have thus designed Clint, a tool to support parallelization through interactive direct manipulation of loop nest visualization (Fig. 1). It presents loop iteration spaces as geometric shapes, dynamic statement executions as solid dots and data dependences as arrows, together with continuous feedback during manipulation. Clint maintains correspondence between the valid source code and the visualization using polyhedral compilation techniques. This interactive approach may ease the design and exploration of program transformations comparing to manual or even semiautomatical code restructuring thanks to immediate feedback, editable history and undoability, allowing programmers to focus on parallelism given that code semantic equivalence is guaranteed by the tool.

II. POLYHEDRAL REPRESENTATION OF PROGRAMS

The polyhedral model is an algebraic representation of a subset of imperative programs that encodes dynamic executions of statements [START_REF] Feautrier | Polyhedron model[END_REF]. It was designed to represent static control parts (SCoPs), which are loops with affine control and memory accesses, i. e., such that conditions and array subscripts are linear forms of outer loop counters and constants. Even with these restrictions, the model captures a broad range of compute-intensive, execution time-consuming loops that can be optimized [START_REF] Wolfe | High Performance Compilers for Parallel Computing[END_REF]. Moreover, extensions to the model make it applicable up to entire functions [START_REF] Benabderrahmane | The polyhedral model is more widely applicable than you think[END_REF].

The model considers dynamic statement instances, which are modeled for each statement as integer points inside a polyhedron called the iteration domain of the statement, that captures loops and branches surrounding it. For the polynomial multiplication loop nest presented in Fig. 2a, it contains twodimensional vectors where the first dimension corresponds to the outer loop on i and the second to the inner loop on j (Fig. 2b). Each vector refers to a particular execution of the statement S. Several compilers have the ability to raise SCoPs to a polyhedral form such as GNU GCC 1 or LLVM 2 .

f o r ( i = 0 ; i < N; i ++) f o r ( j = 0 ; j < N; j ++) S : The statement instances can be executed following a specific order defined by associating each point in the iteration domain with a logical execution date by means of a mapping relation. In case of parallel systems, multiple points having equal dates can be executed simultaneously given that a location of execution (e. g. processor core) is also associated with each instance. Mapping relations are a convenient and expressive way to define a complex composition of program transformations including, e.g., loop interchange, fusion, fission, skewing, tiling, index-set-splitting [START_REF] Girbal | Semi-automatic composition of loop transformations for deep parallelism and memory hierarchies[END_REF]. To ensure semantic equivalence, any mapping relation must not violate data dependences, i.e. preserve the relative order of instances accessing the same memory location.

z [ i + j ] += x [ i ] * y [ j ] ; D S (N ) = i j 0 ≤ i < N 0 ≤ j < N (
A powerful property of the polyhedral approach is to make it possible to compute exact and instance-wise data dependences, and to check the mapping correctness against them [START_REF] Feautrier | Dataflow analysis of array and scalar references[END_REF]. It allows polyhedral compilers to construct correct mapping relations for various objectives automatically, e.g. to efficiently use cache memory or to expose parallelism. For the polynomial multiplication example, a mapping defined in Fig. 3a groups accesses to the same memory address z[i+j] in a single iteration of the outer loop, removing the data dependence between iterations, and thus making this loop parallelizable. The mapping reshapes the iteration domain to the form of Fig. 3b. Once the transformed domains have been computed in the model, the resulting code is generated by building a program that enumerates the integer points within the domain with respect to lexicographic ordering of their coordinates. Several efficient algorithms and tools exist for that purpose including CLooG [START_REF] Bastoul | Contributions to High-Level Program Optimization[END_REF] and CodeGen+ [START_REF] Chen | Polyhedra scanning revisited[END_REF].

θ S (N ) = i j → t1 t2 t1 = i + j t2 = i (
#pragma omp p a r a l l e l f o r p r i v a t e ( t 2 ) f o r ( t 1 = 0 ; t 1 <= 2 * N-2; t 1 ++) f o r ( t 2 = max ( 0 , t 1-N+ 1 ) ; t 2 <= min ( t 1 , N-1); t 2 ++) S :

z [ t 1 ] += x [ t 2 ] * y [ t 1-t 2 ] ;
Fig. 4. Generated code from the transformed domain T S (N ).

1 http://gcc.gnu.org/wiki/Graphite 2 http://polly.llvm.org Fig. 4 shows the final code generated from the transformed domain in Fig. 3b with a parallel loop on t1. The parallelizing optimization corresponds here to a classic loop transformation called skewing [START_REF] Wolfe | High Performance Compilers for Parallel Computing[END_REF]. Although several automatic mapping algorithms were proposed, they are not necessarily resulting in the best parallelizing solution because of the multiple heuristics they use to handle complex underlying processor architectures [START_REF] Bondhugula | A practical automatic polyhedral parallelizer and locality optimizer[END_REF], [START_REF] Bastoul | Contributions to High-Level Program Optimization[END_REF]. Furthermore, since they heavily transform the original program structure, the result is often unreadable for a programmer and thus uncorrectable. Semiautomatic approaches allow programmers to construct and modify a mapping relation freely, but they do not help in finding an appropriate transformation and still require significant expertise to describe it in a special syntax. Our objective with Clint is thus to provide programmers with an interactive tool to ease both parallelism extraction through direct manipulation of visual dependence abstraction and parallelism expression through automatic code generation.

III. THE CLINT INTERACTIVE CODE VISUALIZATION

Design rationale -Clint is a direct manipulation interface designed to (i) help programmers with parallelizing computeintensive programs parts; (ii) ease the exploration of possible transformations; and (iii) guarantee the correctness of the final code. Clint leverages the geometric nature of the polyhedral model by presenting code statements, their instances, and dependences in a scatterplot-like visualization of iteration domains similar to those used in the literature on polyhedral compilation. By making the visualization interactive, it reduces parallelism extraction and code transformation tasks to visual pattern recognition and geometrical manipulations, giving a way to manage the complexity of the underlying model.

During our preliminary design studies, interviewed domain experts expressed the following reasons for showing little interest in semi-automatic approaches: (i) the lack of immediate feedback on the effect of each transformation, (ii) the tight coupling between the code structure and the transformation specification language that complicates composition since each subsequent transformation potentially changes the code structure. Clint addresses these concerns by maintaining the consistency of the visualization independently of the underlying program structure, and by postponing the final code generation, making it possible to assess the results after each step of the transformation. For example, in Fig. 6, both blue and green statements are duplicated after transformation. Therefore subsequent transformations of the "original blue" statement require two operations with a semi-automatic tool while it can still be manipulated as one with Clint.

Finally, Clint introduces a navigable history view of transformations to favor try-and-fail strategy in the search for the best optimization. It can also save and reuse transformation scripts in order to share best practices in parallelizing optimizations. Overall, we believe that this approach may enable finer reasoning in terms of statements and instances rather than in terms of loops. is represented in Fig. 5-left: axes represent iterators; the green shape delimits the iteration domain of the statement; the points inside the shape represent the statement instances, linked with arrows to denote data dependences between them. To convey information about the original or transformed execution order of statement instances, Clint visualizes a mapped iteration domain: a set of points obtained by applying the mapping relation to the integer points of the iteration domain like in Fig. 5-right for the polynomial multiplication example. Colors and thickness adjustments provide visual feedback on the mapping result: thick red arrows depict violated dependences while thick green axes represent parallel loops.

Visualizing Iteration Domains of Program Statements

Multiple statements in the same loop are represented as different colored polygons of the same shape slightly displaced according to the execution order (Fig. 6). They share an axis that corresponds to their common loop. Multi-dimensional iteration domains are split up into two-dimensional projections and displayed as a scatterplot matrix [START_REF] Heer | A tour through the visualization zoo[END_REF]. Multiple statement instances, which differ only in the dimensions that were projected out, are depicted as a single point with a shade of gray representing the relative number of points. Two-dimensional visualization allows direct manipulation with a standard 2D input device (e.g. mouse) [START_REF] Beaudouin-Lafon | Designing interaction, not interfaces[END_REF] and guarantees the consistency of the interface if the number of dimensions outgrows three.

Direct Manipulation of Statements and Instances -This geometrical representation of program statements affords direct manipulation of either the entire iteration domain (polygons) or particular statement instances (points), making parallelization just a matter of aligning dependence lines. Clint's direct manipulation builds on the geometry-related vocabulary used in the polyhedral compilation community: points or polygons can be dragged and thus shifted (Fig. 6); polygons can be fused or split apart (Fig. 7), or skewed in space (Fig. 5) to reshape the respective iteration domains. Each graphical action is mapped to a sequence of program transformations that, if applied to the code, would change the program structure so that its visualization corresponds to the displayed one.

Clint Interface -As shown in Fig. 1, Clint combines three components: (i) the interactive visual representation of the program statements; (ii) an editable and navigable transformation history view; and (iii) the source code editor, in which the

f o r ( i = 0 ; i < 3 ; i ++) f o r ( j = 0 ; j < 3 ; j ++) { A[ i ] [ j ] = 1 / 2 * (A[ i ] [ j ] + A[ i ] [ j + 1 ] ) ; B[ i ] [ j ] += A[ i -1 ] [ j ] ; } f o r ( j = 0 ; j <=2; j ++) B [ 0 ] [ j ] += A [ ( 0 ) -1 ] [ j ] ; #pragma omp p a r a l l e l f o r \ p r i v a t e ( j ) f o r ( i = 1 ; i <=2; i ++) f o r ( j = 0 ; j <=2; j ++) { A[ i -1][ j ] = 1 / 2 * (A[ i -1][ j ] + A[ i -1][ j + 1 ] ) ; B[ i ] [ j ] += A[ i -1][ j ] ; } f o r ( j = 0 ; j <=2; j ++) A [ 2 ] [ j ] = 1 / 2 * (A [ 2 ] [ j ] + A [ 2 ] [ j + 1 ] ) ;
Fig. 6.

Interactive shift transformation: the lighter polygon is dragged right to make dependence lines vertical so that they do not span between different iterations on i. Thanks to consistent visualization, statements can be manipulated as if they were not split between two loops. program structure can be changed manually. The three views are synchronized and updated automatically to reflect actual changes. The source code view can however be "locked" to the original version as the transformed code tends to quickly become complex and unreadable. The history is represented as a transformation script (using Clay [3] language) with coloring and grouping to reflect visual manipulations that resulted in multiple transformations.

f o r ( i = 0 ; i < 3 ; i ++) f o r ( j = 0 ; j < 3 ; j ++) { sum += A[ i ] [ j ] ; A[ i ] [ j ] = 0 ; } f o r ( i = 0 ; i < 3 ; i ++) f o r ( j = 0 ; j < 3 ; j ++) sum += A[ i ] [ j ] ; #pragma
Architecture and Implementation -Clint's implementation relies on a collection of research tools and libraries that are the building blocks for polyhedral compilers and that usually work together as a single black-box [START_REF] Bastoul | Contributions to High-Level Program Optimization[END_REF]. Clint is a standalone application that consolidates editing facilities, feedback from the tools and output into a single consistent interface. It takes the original source code and uses Clan to raise it to the polyhedral model, it then maps direct manipulations to Clay transformation scripts and represents dependence violation feedback from Candl. Finally, it relies on CLooG to generate the C code with OpenMP pragmas for parallelism (Fig. 8).

IV. RELATED WORK

Visual Representations for Polyhedral Model -Spatial projections of the polyhedra are widely used in the domain literature. Several polyhedral libraries and automatic tools include visualization functionality for iteration domains and dependences, such as VisualPolyLib [START_REF] Loechner | Polylib: A library for manipulating parameterized polyhedra[END_REF] or LooPo [START_REF] Griebl | Automatic Parallelization of Loop Programs for Distributed Memory Architectures[END_REF]. 3D iteration space visualizer [START_REF] Yu | Loop parallelization using the 3d iteration space visualizer[END_REF] relies on an automatic search for parallelizing transformation. Tulipse [START_REF] Wong | Tulipse: a visualization framework for user-guided parallelization[END_REF] integrates dependence visualization into the Eclipse IDE. Our work extends these static visualizations by providing a mapping between direct manipulation of the visualization and code transformation.

Semi-Automatic Program Transformations -Several frameworks expose a high-level interface on top of polyhedral libraries through scripting languages, for example UTF [START_REF] Kelly | A framework for unifying reordering transformations[END_REF] or URUK [START_REF] Girbal | Semi-automatic composition of loop transformations for deep parallelism and memory hierarchies[END_REF], more recent ones focusing on transformation sequencing like Clay [START_REF] Bastoul | Contributions to High-Level Program Optimization[END_REF]. Clint relies on it to convert visual actions to mapping relations and to provide programmers with a graphical interface that facilitates analysis tasks.

V. CONCLUSION This paper introduced Clint, a direct manipulation visualization tool for parallelism extraction and expression. It builds on polyhedral compilation techniques to feature automatic data dependence analysis with real-time feedback, and to ensure correctness of the generated code. Clint leverages the geometric nature of the polyhedral model to complement conventional code editing with interactive visualization.

The underlying model limits the applicability of Clint to programs with static control flow, but still covers a broad range of compute-intensive loops, in which most of execution time is spent. Extending the visualization to capture entire programs would however lead to visual cluttering issues, especially in the case of numerous data dependences. Zoomable interfaces and on-demand animated transitions [START_REF] Dragicevic | Gliimpse: Animating from markup code to rendered documents and vice versa[END_REF] between any combination of code and visualization may help to reduce cluttering as well as to provide a smoother link between the original and transformed code. Another future work direction is to investigate the use of interactive visualizations for learning parallelization by featuring dynamic feedback and interactive guidance through manipulation restrictions/enhancements (e.g. pseudo-haptic feedback [START_REF] Lécuyer | Feeling bumps and holes without a haptic interface: The perception of pseudo-haptic textures[END_REF] or semantic pointing [START_REF] Blanch | Semantic pointing: Improving target acquisition with control-display ratio adaptation[END_REF]).

Finally, despite being widely used in the literature, the visualization of programs used in Clint was never evaluated. We conducted a preliminary study whose results suggest that users are able to reliably reconstruct code structure from the visualization and vice versa, which in turn suggests that the visualization is a complete representation of a program. We also asked some potential users of Clint to parallelize codes of different complexities, which they did faster and with greater accuracy than with usual methods (code editing and semiautomatic tools). While these results should be considered as preliminary, our interactive visual approach is a promising solution to ease the construction and adjustment of complex parallel optimizations of imperative programs.

Fig. 1 .

 1 Fig. 1. Clint interface with synchronized views: (1) interactive visualization of loop iteration spaces; (2) editable transformation script history (3) code editor switching between transformed and original code to maintain context.

Fig. 2 .

 2 Fig. 2. Iteration domain of a statement.

2 TFig. 3 .

 23 Fig. 3. Time mapping example and resulting transformed iteration domain.

  -Clint's visualization of a program statement is a twodimensional projection of its iteration domain as a polyhedron. The statement S of the polynomial multiply example of Fig.2a

Fig. 5 .

 5 Fig. 5. Performing a skew transformation to parallelize polynomial multiplication loop by deforming the polygon. The code is transformed from its original form (left) to the skewed one (right) automatically with OpenMP pragmas generated for parallel loops depicted as thick green axes.

Fig. 7 .

 7 Fig. 7. Interactive complementary transformations to split statements between loops (left to right) or fuse them in a single loop (right to left). Each loop nest is represented as a separate coordinate system.

  Clint Software Architecture and Interaction Loop. Previously the user interacted with all tools (in black); Clint provides a single interface (in blue).

	Clan code	Candl polyhedra	Clay dependencies	polyhedra transformations	CLooG code
	direct manipulation visual feedback	Clint
	code
	Fig. 8.


L Grenet

M Jamet

P Noé

V Calvo

J.-M Hartmann

L E Nistor

B Rodmacq

S Auffret

P Warin

Y Samson

Spin injection in silicon at zero magnetic field

In this letter, we show efficient electrical spin injection into a SiGe based p-i-n light emitting diode from the remanent state of a perpendicularly magnetized ferromagnetic contact. Electron spin injection is carried out through an alumina tunnel barrier from a Co/Pt thin film exhibiting a strong out-of-plane anisotropy. The electron spin polarization is then analyzed through the circular polarization of emitted light. All the light polarization measurements are performed without an external applied magnetic field, i.e., in remanent magnetic states. The light polarization as a function of the magnetic field closely traces the out-of-plane magnetization of the Co/Pt injector. We could achieve a circular polarization degree of the emitted light of 3% at 5 K. Moreover this light polarization remains almost constant at least up to 200 K. © 2009 American Institute of Physics. ͓DOI: 10.1063/1.3064135͔ Spintronics aims at manipulating both charge and spin of carriers and has been growing tremendously fast for the past 10 years. 1 Indeed it has been predicted as an alternative solution to the conventional complementary metal oxide semiconductor technology that will reach physical limitations in the near future. 2 This subject relies on efficient spin injection into semiconductor heterostructures. Moreover for integration purposes, spin injection should be carried out electrically. So far, electrical spin injection has been widely demonstrated in III-V semiconductors using both ferromagnetic metals and magnetic semiconductors. The spin polarization of injected carriers is then given by the degree of circular polarization of the light emitted by a GaAs based spin-light emitting diode ͑LED͒ ͑Ref. 3͒ through dipolar selection rules. The most spectacular spin injection and optical detection were performed from BeMnZnSe, 4 Fe, 5 or CoFe ͑Ref. 6͒ injectors. Only recently, spin injection in silicon has been achieved using ferromagnetic metals. Spin polarized hot electrons were injected in a silicon wafer, and spin detection was achieved by electrical means using two ferromagnetic electrodes in a giant magnetoresistancelike geometry. 7 In parallel, Jonker et al. 8 demonstrated spin injection in the conduction band of silicon by measuring the circular polarization of a Si p-i-n LED. It should be stressed that silicon appears as the most promising candidate for spintronic circuitry and quantum computing systems since very long diffusion lengths and coherence times are predicted in this material, owing to its low atomic mass ͑low spin-orbit coupling͒ and inversion symmetry of the crystal itself. Moreover the dominant naturally occurring isotope has no nuclear spin, ruling out spin relaxation due to hyperfine interactions. 9,10 In this letter, we demonstrate spin injection into silicon using the remanent state of a perpendicularly magnetized ferromagnetic semitransparent contact. The electron spin polarization is detected optically using a Si based spin-LED. The active region is a fully strained Si 0.7 Ge 0.3 quantum well ͑QW͒. Compared to pure silicon spin-LEDs, the presence of germanium enhances the spin-orbit coupling, which is required to observe circularly polarized light from the recombination of spin polarized injected electrons with holes. Spin injection from the Co/Pt ferromagnetic metal is carried out in the tunneling regime through an alumina barrier to overcome the conductivity mismatch obstacle. 11 Dipolar selection rules in spin-LEDs require perpendicular magnetization. In most experiments on spin injection, ferromagnetic films exhibit in-plane magnetization due to shape anisotropy. Hence strong external magnetic fields ͑2.2 T in the case of iron͒ are applied to saturate the film magnetization perpendicular to the QW plane. This may be the cause of many spurious effects such as magnetic circular dichroism ͑MCD͒ in the experimental setup itself or in the spin-LED heterostructure. Moreover the Zeeman effect in the active region artificially increases spin lifetimes. Very few experiments have been carried out in remanent magnetic states ͑i.e., with no applied field͒ into GaAs based spin-LEDs. Adelmann et al. 12 reported a remanent spin injection of about 5% at 2 K using the ␦-MnGa Schottky contact. Gerhardt et al. 13 demonstrated remanent spin injection using Fe/Tb multilayers with a circular polarization of 0.75% at 90 K. Finally Sinsarp et al. 14 used a FePt/MgO spin injector and obtained a lower estimate of spin injection of 1.5% at room temperature. In this work, the Co/Pt ferromagnetic film exhibits strong perpendicular magnetic anisotropy, allowing us to make spin injection in different remanent states.

We used the reduced pressure chemical vapor deposition in order to grow at 26.6 mbar the following stack on a slightly p-type doped Si͑001͒ substrate ͑resistivity in the 7-10 ⍀ cm range͒: boron-doped Si ͑10 19 cm -3 ͒ 500 nm/ intrinsic Si ͑50 nm͒/SiGe 30% ͑10 nm͒/intrinsic Si ͑50 nm͒/ phosphorous-doped Si ͑10 18 cm -3 ͒ 50 nm. The growth temperature was around 650 °C in order to avoid any significant surface roughening of the Si 0.7 Ge 0.3 layer during its capping with Si. The gaseous precursors used were silane a͒ Author to whom correspondence should be addressed. Electronic mail: louis.grenet@cea.fr.

͑for the Si layers͒ and dichlorosilane+ diluted germane ͑for the SiGe layer͒. For the boron ͑phosphorus͒ doping, we have used diborane ͑phosphine͒. High purity hydrogen with a flow of a few tens of SLM ͑SLM denotes standard liters per minute͒ was used as the carrier gas. This SiGe based spin-LED is subsequently protected with a 2 nm thick thermal SiO 2 layer for air transfer. Prior to the deposition of the tunnel barrier and ferromagnetic layer, this protective layer is etched using hydrofluoric acid and rinsed in de-ionized water. A 2 nm thick alumina tunnel barrier is then deposited at room temperature using radio frequency magnetron sputtering. The ferromagnetic layer is grown by direct current magnetron sputtering at room temperature through a mechanical mask defining 900 m wide mesas. It consists of a 1.6 nm of Co covered with 3 nm of Pt. The as-grown ferromagnetic layer exhibits an in-plane magnetization. After annealing the whole heterostructure at 420 °C for 1 h and 30 min the Co layer magnetization points out of plane. The whole structure as well as the band structure 15 is sketched in Fig. 1.

Superconducting quantum interference device ͑SQUID͒ measurements with the applied field perpendicular to the film plane are displayed in Fig. 2. In Fig. 2͑a͒ are hysteresis loops recorded at different temperatures. As shown in Fig. 2͑b͒, the coercive field strongly decreases with temperature by 88.3%. However low coercive fields at room temperature are required to switch easily up and down the magnetization of the spin injector in future spintronic applications. In the inset of Fig. 2͑b͒, we have plotted the temperature dependence of the saturation magnetization. The saturation magnetization is the one of bulk cobalt at low temperature. It however decreases by 30.5% when temperature increases from 5 to 400 K. This temperature dependence follows a Bloch law and we could roughly estimate a Curie temperature of 900 K, which is much lower than bulk cobalt ͑T C = 1400 K͒. This may be due to dimensionality effects or more probably to Co-Pt intermixing during sample annealing. Since the coercive field is correlated with the magnetization, its sharp decrease with temperature can be explained by the magnetization temperature evolution. However, we could not find a simple power law for H C versus M S .

We performed both photoluminescence ͑not shown͒ and electroluminescence ͑EL͒ measurements on the SiGe spin-LED at various temperatures using a cooled InGaAs photomultiplier tube. The EL was excited with square electrical pulses at low frequency ͑63.7 Hz͒ with a 1:1 duty cycle, and the signal was recorded using a lock-in technique. The circular polarization P = ͑I + -I -͒ / ͑I + + I -͒ is determined using a rotatable quarter waveplate and a fixed linear polarizer placed in front of the detector. In Fig. 3͑a͒ are reported EL spectra recorded at 5 and 77 K for a saturated remanent state of the ferromagnetic contact. In addition to the silicon EL ͑not shown͒, the spectra are dominated by two main features at ϳ950 meV corresponding to the no phonon ͑NP͒ recombination line and ϳ900 meV corresponding to its TO phonon replica. 15 Both lines are clearly circularly polarized as a consequence of momentum transfer between the spin polarized injected electrons and the emitted photons. In the following, we only consider the highest circular polarization of the NP line. We could reach a maximum of P =3% at 5 K. Moreover this circular polarization remains almost constant up to 200 K. Above this temperature, the EL signal from the SiGe QW becomes too weak to measure the circular polarization. In order to rule out any extrinsic effects, we performed EL measurements using a gold semitransparent contact and inserting in the optical path a thin Si plate on which we deposited the Co/Pt ferromagnetic layer. Hence we could confirm that MCD in the Co/Pt layer in a saturated remanent state is negligible. Moreover, Zeeman effect in the Si 0.7 Ge 0.3 QW could be ruled out since the stray field from the saturated Co/Pt contact is negligible as pointed out by Gerhardt et al. 13 In Fig. 3͑b͒, we have plotted the EL circular polarization for different magnetic remanent states of the ferromagnetic contact. In this experiment, EL measurements are performed without applied magnetic field. The reported magnetic field in the graph corresponds to the magnetic field required to reach a given remanent state of the magnetization of the Co/Pt layer. The circular polarization exactly traces the SQUID measurements, indicating that the spin orientation of the injected electrons that radiatively recombine in the SiGe QW reflects that of the electron spin orientation in the Co/Pt film. The measured circular polarization corresponds to a lower bound of the actual spin injection in silicon. Indeed spin relaxation ͑with a characteristic time sf ͒ may take place during the transit time of electrons in the silicon top layer as well as during the radiative lifetime ͑ r ͒ in the SiGe QW. Neglecting spin relaxation in the silicon top layer, the actual electron spin polarization achieved in silicon is given by P͑1+ r / sf ͒, where the ratio r / sf strongly depends on temperature and the exact electronic structure of the QW. We estimate the total lifetime in the QW; it is tens of nanoseconds, which is very short. However we have no clear estimate of the spin relaxation time ͑ sf ͒ in the SiGe QW. The determination of sf requires time-resolved spin dynamics measurements. We can only mention electron spin resonance measurements performed by Jantsch et al. 10 on Si 1-x Ge x ͑x Ͻ 0.1͒ QWs grown pseudomorphically between Si 1-y Ge y ͑y Ͼ 0.2͒ barriers giving sf Ϸ 1 s. In our case, with higher Ge content and strain, we believe that spin relaxation time should be shorter, although we cannot give a lower bound. Further experiments are thus needed to estimate the actual spin injection efficiency into silicon.

In summary, we have shown tunneling spin injection into silicon from the remanent state of a ferromagnetic Co/Pt contact. A maximum circular polarization of 3% could be achieved at 5 K and remains almost constant up to 200 K. Moreover spin injection was performed for different magnetic remanent states of the ferromagnetic layer, and the associated circular polarization exactly traces the SQUID mea-surements as expected for spin polarized electrons injected from the Co/Pt contact.

͑Received 13

 13 November 2008; accepted 12 December 2008; published online 20 January 2009͒

FIG. 1 .

 1 FIG. 1. ͑Color online͒ ͑a͒ Schematic device structure of the SiGe spin-LED with the Co/Pt ferromagnetic top electrode and ͑b͒ associated simplified band diagram under bias in the injection regime.

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 132.168.65.215 On: Wed, 09 Apr 2014 06:30:33

The authors would like to thank L. Notin, A. Brenac, and D. Leroy for their technical assistance with the sample growth.