C T A Schmidt
email: archiveouverte@ensam.eu

Keywords: Artefacts, Difference, Epistemology, Diversity (human), Identity, Relation, Personhood, Philosophy (analytical), Transhumanism

Many have bowed before the recently acquired powers of 'new technologies'. However, in the shift from tekhne ¯to tekhnologia, it seems we have lost human values. These values are communicative in nature as technological progress has placed barriers like distance, web pages and 'miscellaneous extras' between individuals. Certain values, like the interpersonal pleasures of rendering service, have been lost as their domain of predilection has for many become fully commercially oriented, dominated by the cadence of profitability. Though the popular cultures of the artificial have surged forth to deliver us from the twentieth century, they have enabled some very superfluous dreaming-Man has succumbed to the Godly role of simulating himself and creating other beings. Communication is replaced by machines, services are rendered via many automated devices, procreation has entered the public sphere, robots and entertainment agents educate our youth and mesmerising screen-integrating 'forms of intelligence' even think for us. As such, this so-called culture threatens the very values Man constructed in the nineteenth and twentieth centuries to guide himself into the future. But what if the phenomena mentioned just reflect our new values? The author presents an investigation into this cultural shift, its impact on human practices with regards the mind and the body and evokes some pros and cons of generally accepting the 'Culture of the Artificial'.

Introduction

When it comes to the production of artefacts, our environment-rich with technical objects-turns us into artefacts; in the current context, our pathological obsession on technological conformism may lead one to think of oneself as an artefact (during an invited (2011) talk in French entitled ''La communication transhumaniste? Discernements et de ´fis'' at the 1st Conference of the French Transhumanistic Society, I gave examples of this logic which actually stem from [START_REF] Schmidt | Of robots and believing, minds and machines[END_REF][START_REF] Schmidt | Robots, Dennett and THE Autonomous. A terminological investigation, minds and machines[END_REF]Schmidt and 2007)). This may sound strange but people do indeed generally try to remodel their daily lives, capabilities and appearances to keep up with others around them, as just anyone may notice. Question: Should biodiversity studies include the redesigning of Man? Regardless of the way one decides to answer here, it would seem obvious that man's simulating himself has become a culture, or even a religion. The scientific references for new future situations do not often directly portray this activity as being ritualised. My personal observation from over the last 15 years allows me to confirm this tendency about Man's self-simulation. Cultural variations on such activity do exist and are too abundant to treat here.

The way I see things, the artefactualisation of the human species may be considered part of the evolutionary process. Some may not agree with this statement; for many, the matter dealt with herein appears ideological or perhaps irrelevant to down-to-earth daily concerns… for the time being. I intend to explore the types of arguments put forth on both sides of the equation in order to illustrate the incomprehension that exists between the protagonists responding from opposing viewpoints with regard to the question of whether or not we should redesign the human being, a fast growing issue.

More specifically, I will address issues involving advanced Artificial Intelligence, Humanoid Robotics and Transhumanism without wishing to put forth a decision why we should or should not accept new forms of life similar to our present state or those that deviate from it. This said, uprising nations should particularly be wary of deviant grassroots industrial practices (such as organ farming or operations that do not reflect the standards of the Establishment), now and in the future. Since the nature of body or personhood change is limited in some movements to behavioural or physical aspects, the actual change seems spectacular (i.e. a smart silicon implant) in comparison with changing one's body through exercise. The reader will not find any statistics here, just facts and principles about modifications aiming towards new creations (i.e. ethically impacted material such as female feticide are left out as it does not constitute an improvement on the being in question, but rather just an 'existential preference'.

Techniques for human modification

What this means is that there are basically two approaches to going about the artefactualisation of the human species. In fact, the evolutionary process itself has changed and possibly further diversification of it may come about especially if humans play a role in guiding evolution. The current two approaches are separate in that their starting points are not the same. On the one hand, the roboticsbased approach generally uses many components that are mechanical in nature (traditional hardware) though there is a growing tendency to accept organic elements into these constructions. The reasons for using organic stuffs in the robotics sphere of intervention in society are various: either they are less costly, increase functionality, render the resulting 'machine' more lifelike, less harmful to the environment, or they provide jobs to the local workforce. On the other hand, the transhumanist approach begins rebuilding man with one single very familiar component, the human body-just like yours or mine-assuming it is fully natural. The idea is to use technological advances to modify the body or brain in hopes of creating a desired effect. This could entail introducing various entities into the body for a variety of reasons: molecules (i.e. using metabolic control for 'slimming', anti-ageing medicine to stay young or live long), electronic chips (i.e. in the brain to help one understand better or remember more, in the eyes to improve sight, etc.), bionics for increased power, or all other implants.

Perhaps a minor detail would be the difference between implants and transplants. The former generally take the currently-used-to state of the individual to one she/he does not know in order to experience it-picture the average person having Steve Austin's bionic ability to lift and throw heavy objects!-whereas the latter aims at bringing one back to a state she/he has always been used to but has recently lost-an elderly person having a hip replacement. The only similarity between the two is that they both augment the person's present state.

Let us get back to the robotics versus transhumanism distinction. Although different, it is important to point out that there are similarities: for both approaches, it is the desired effect that leads to the design of a new or even novel being, which means there is a certain wilfulness driving towards a new world. I do not think this drive is new; it is just the techniques used that may surprise people. Change is a concept that the notion of Man has always integrated, this is the reason we are part of the world's evolutionary cycle.

But it would seem that this short-term aspect of evolution is mainly behaviour based, thus there is limited change to the identity of what it means to be human. The concept of being human does entail a highly social element to it, as well as a cultural one: it is not in individually modifying the bodies of the members of society that one can change the relationships they enjoy or detest. This said sustained corporal change over time could well have an effect on relations in society.

Difference and the concept of Man

The concept of Man would of course become different, but to what extent? Perhaps the thing that society is calling out for here is a concept of humans that is more material in nature when compared to the current idea of what it is to be human. The belief that we could/should/must modify our very own physical existence may mean that the immaterial-social, psychological, cultural and spiritual-aspects of our lives have become less important to us. Would such a statement be too simplistic or is it part of our new reality? Those working in advanced Artificial Intelligence, Cognitive Robotics, Neuro-evolution, and transhumanistic-type technologies generally do not delve into the intricate questions of love, faith or the respect of the Other in society, all of which directly concern the human immaterial sphere. They are not supposed to be intimately concerned with such matters, nor are other scientific fields. One could nevertheless be very mistaken in saying that these matters are not on scientists' agendas (cf. infra my mention of the advances towards this in 1. Presuppositions); how could they ever hope to do better than man if they cannot even copy these facettes of the human race? We could indeed conclude for the time being that the concept of being human today means being more physically human than 100 years ago.

Relation and the concept of Man

So the concept of Man has evolved. Does this operate on the relations this concept has with other concepts in its vicinity? The concepts of Nature and Artefact would be pertinent to explore here. The fact that we today accept to tamper with mother nature's 'products' is not new but the acceptance to direct this action on our physical and cognitive existences has increased exponentially. The place of nature-highly recognised as essential-in traditional cultures would have diminished. However, one is entitled to disagree with strong statements, in which case it could only be said that our relation to nature would have been altered only slightly. The important matter is in asking why this change suddenly became necessary and what our new relation with nature, however, close it is to the last version, means to us in the future.

As for the link to the concept of artefact, the shift does seem more radical, although Man has always accepted to handle his own destiny; hence the expression 'he is a selfmade man'. The tie between Man and being man-made has been strengthened in the consciousness of members of society, perhaps paradoxically. A 'self-made man' always referred to the self-assurance, aspirations, intellectual stamina and other paraphernalia of the purely psychological composition of the individual whereas now, we are able to apply it to his physical composition. The imminence of the merger between robotic products and the physically oriented propositions of transhumanism can now be addressed as more than just an eye-opening thought. If one prefers lesser-alarming realities for measuring one's consciousness of the conceptual relation between men and artefacts, one could examine the simple layman's example of the use of steroids: first they were used practically joyfully, then considered cheating, now they are amalgamated to deadliness-and this shift happened over a relatively short period of time whilst the effects of their use remained stable. Will our judgement on what can and what cannot be considered an artefact be this rapidly affected too?

Identity and the concept of Man

If the relations among the adjacent concepts in the system are modified because of human modification becoming abundant, what does this bring forth for the identity of man in the ecosystem? Imagine that we welcome this slide toward the increased physical definition of Man: some serious interrogations suddenly become visible. The fact that man would have the opportunity to change the very concept of himself in this manner, and that this would have a real effect on his surroundings, proves that homo sapiens would control his own 'conceptual environment' and that the techniques discussed here would be a mere side effect of his existence (i.e. other techniques could be used to sustain the developments sought). This would mean that individuals really would have obtained an overwhelming level of power vis-a`-vis their past and vis-a`-vis their counterparts.

3 Shouldn't we be against greater human diversity?

In the hypothetical system just described, the weaker are bound to suffer more. Is this the type of homo sapiens we wish to become?

The identity of others (and thus Selfhood too) would be heavily affected in such a world. The identity of the Self would be equated to a very heady position-practically Godhood. But today, we do have the 'magnificent opportunity' to actually apply this ill-formed logic to our lives ourselves.

So should biodiversity include the redesign of man? As I said, the key to strengthening the argument against modifying Man has to put forth practical ideas on how and why not to modify man in such an era-many are doing it (eugenics, implants…), though perhaps not to the point of becoming cyborgs. The undemocratic nature of such practice is not stopping people because such an argument remains quite ideological for most or at least irrelevant to down-to-earth daily concerns.

The way in which they-the artificial or modified beings-would seem (too) different from the average human today is in the values they would, perceivably, be able to share and apply. Because of the hypothetical differences in the importance of (traditional) human values and their application between original men and the nonorganic originating persons, one may not wish to see these latter caring for one's children or for the elderly. One may have difficulty trusting the moral judgements of a nonnatural neighbour or artificial person. Imagine one's multiply implanted neighbour is the judge in a widely mediated homicide case; this could be a real media 'scoop', which would help picture the situation here. The same would go for our discovering the judge was a machine.

The practical measures necessary to supporting the view that biodiversity should not include the redesign of Man would entail, among other things, avoiding simulation in all its forms. This measure on simulation could be presented as just general advice, negotiations for special cases determined by a set of criteria, or out-right prohibition, according to the adopted political stance. The important aspect here is the urgency of the question as simulations are changing our vision of ourselves and our world. diversity If we take up the judge example once again, the fact that a human judge with silicon chip implants in the brain or a purely robotic judge can be hooked up to all the sources of jurisprudence possible in a continuous manner may one day comfort the average citizen in the decision he/it makes. It goes without saying that the sources the 'judge' has access to will have to be limited to a select official set. If accepted, this would be a very practical solution.

Those supporting the view that biodiversity should include redesigning humans have to develop strategies on the ideological and political levels in order to seriously further their cause, that of artefactualising humans, because after all, Man would be an 'artefactual object' if he was remodelled in the ways spoken of above. Furthermore, our judge example is simply that: an example. One cannot change society's strategies on a political level based on a mere practical example. Those that wish to promote the vision of a widened biodiversity in which homo sapiens would be one of the species involved would have to either directly modify the moral position of humans in the world equilibrium (difficult?) or show the strategic advantages to becoming robotic individuals (acceptable?), transhuman or posthuman (meaningless?). This may help people reexamine their traditional stances. It may help here if I make things explicit; the judge example would be a very practical solution, if accepted by all or by the majority, as long as average citizens are consulted.

Looking at the transhumanist movement will reveal that the values put forth, whether one sees them to be acceptable or not, are done so within the framework of specific basic conditions according to [START_REF] Bostrom | Transhumanist Values, Journal of Philosophical Research, Special Supplement on 'Ethical Issues for the Twenty-First Century[END_REF], of Oxford University: 'Global Security, Technological Progress, Wide Access'. 1 Any sensible being shares these conditions and would like to have them protected, so the movement is not so off-tilt as some might say it to be. Have we not collectively been tampering with nature for a long time already? Man has always had the tendency to 'diversify' in one way or another, more on moral grounds than on physical ones.

Accepting to go along with such change would be a strategic move if it were to be used to unite people, communities, etc. Allowing only weaker members of society to better themselves would enable them to gain back their dignity. But would creating laws prohibiting naturally endowed persons access to changes be unfair? It is clear that if the biodiversity of man is to be accepted by the average citizen, any sensible discourse on the matter would have use the legal/moral or non-physical level as its cornerstone.

When one considers the argumentation necessary to change things, it is tempting to just say that the physical aspects of human life are quite malleable in comparison with its non-manifest 'components'. But Bostrom gives us an indication of the little our considering transforming human beings entails and where to look for inspiration to change mindsets. In his own words, 2 he mentions that the necessary ideals are found outside of our bios. We must therefore act on our logos to better fathom the advent of change. It is only if we focus on human reason that we will be able to accept our own redesign.

5 Dialogue for dealing with the cultural shift towards the artificial

The reader may find that I fail to correctly transcend the practical aspects of modifying man to develop sound arguments on the ideological and political levels of why to expand human diversity. However, pulling one way or another was not the goal here. The current discussion in this area of thought reminds me of Ricoeur's [START_REF] Changeux | Dialoguer avec tous, et d'abord avec soi'', dossier 'Disparition: Paul Ricoeur, Philosophe de tous les dialogues[END_REF] stance on the impossible adjustment between our finite body and our reason being infinitely open. Although the two levels of discourse are complementary to one another, their refusal to blend is what leads to our mistakes, our miscalculations, etc. and renders the whole process of decision-making fallible. This partially because we find ourselves before dilemmas such as the one in this question: should the 'ideal' humanoid have all the weakness man has? Do people need to be more machinelike? Today it is clear that a certain number of values are perceived by at least some people as lost. According to Mukherjee (2007) of the Centre for Human Values at the Indian Institute of Management in Calcutta, the distinction between necessity and luxury has been greatly blurred by the effects of globalisation. He states that 'human consciousness is spellbound at the marvels of modern technology […]. Caught as if in a trance, in the maze of 1 Bostrom's article culminates in a ' 'Table of Transhumanist Values'', cf. Bostrom (2005), ''Transhumanist Values'', Journal of Philosophical Research, Special Supplement on 'Ethical Issues for the Twenty-First Century', p. 13. Transhumanism does not require us to say that we should favour posthuman beings or human beings, but that the right of way of favouring human beings is by enabling us to realise our ideals better and that some of our ideals may well be located outside the space of modes of being that are accessible to us with our current biological constitution''. Cf. BOSTROM N., ibidem, p. 8. ceaseless activity at an alarming pace, hardly do we feel the need or urgency to reflect on such deeper questions as the purpose of life and the meaning of our actions.' 3 The 'religion' of simulating or modifying oneself will be the inspiration of a grassroots industry, if it is not already. I do hope to have provided elements essential to engaging dialogue on the higher-order matters herein addressed; taking action in the future depends on the knowledge we are able to create through dialogue. Grassroots industry based on a religion of simulating or modifying oneself obviously entails that one transcends the opposing standpoints presented here: either being for or against modifications to elements that do not originate in the initial material used, for example, injecting human entities into the machine or injecting mechanistic entities into the human body (cf. the writings of [START_REF] Warwick | Cyborg[END_REF]). New forms of organisations need to be driven by vigilance either to protect the unchanged or to accompany the being, whatever he, she or it is through the modifications.

  2 ''The realm of posthuman values does not entail that we should forego our current values. The posthuman values can be our current values, albeit ones that we have not yet clearly comprehended.

C. T. A. Schmidt (&)Le Mans University & ENSAM-ParisTech Angers, Laval, France e-mail: Colin.Schmidt@univ-lemans.fr


Guillaume Maitrejean
email: guillaume.maitrejean@ujf-grenoble.fr

Amine Ammar
email: amine.ammar@ensam.eu

Francisco Chinesta
email: francisco.chinesta@ec-nantes.fr

Simulating microstructure evolution during passive mixing

Keywords: Passive mixing, Microstructural approach, Area tensor, Closure relation

The prediction of microstructure evolution during passive mixing is of major interest in order to qualify and quantify mixing devices as well as to predict the final morphology of the resulting blend. Direct numerical simulation fails because of the different characteristic lengths of the microstructure and the process itself. Micro-macro approaches could be a valuable alternative but the computational cost remains tremendous. For this reason many authors proposed the introduction of some microstructural variables able to qualify and quantify the mixing process at a mesoscale level. Some proposals considered only the effects induced by the flow kinematics, other introduced also the effects of shape relaxation due to the surface tension and coalescence. The most advanced integrate also the break-up process. However, the derivation of the evolution equations governing the evolution of such microstructural variables needs the introduction of some closure relations whose impact on the computed solution should be evaluated before applying it for simulating complex mixing flows. In this work we

Introduction

When two immiscible fluids are mixed, the induced properties of the mixture are strongly dependent on the resulting microstructure, whose length scale is much smaller than the one associated with the macroscopic flow. Indeed droplets of various size and shape are formed and complex phenomena such as deformation, break-up and coalescence occur. Passive mixing implies two immiscible fluids having the same viscosity and density.

In many engineering processes, the determination of the mixing rate is very important because it is directly linked to the prediction of the mixing time and that is why the study of liquid-liquid laminar mixing has been carried out for many years and many rheological tools have been developed.

Many engineering processes benefit from good bulk mixing. In such processes, determination of mixing rate is of importance in terms of understanding and prediction of mixing time. That understanding allows defining new flows and the associated processes for improving the mixing rate, but one could also try to optimize other parameters related to the microstructure as for example the characteristic length, the shape and orientation of the discrete phase, ... Mixing rates have often been quantified by a mixing time, i.e. the time required to reach a certain degree of homogeneity [START_REF] Rielly | Mixing processes for agricultural and food materials: Part 4, assessment and monitoring of mixing systems[END_REF]. Most studies have based on the concentration evolution at a small number of sample points [START_REF] Hall | The mixing rates of highly viscous Newtonian and non-Newtonian fluids in a laboratory sigma-blade mixer[END_REF]. However the sensibility of the results to the number and locations of the probes is its main limitation. An alternative way of quantifying the mixing rate can be derived from the consideration of only one of the basic mechanisms of mixing processes: the increase of the material interface due to the fluid mechanics in absence of interfacial tension and molecular diffusion. The molecular diffusion leads to smooth concentration gradients across the interface. However this mechanism only becomes significant when the interface has increased significantly. It is therefore expected that the overall mixing rate will be closely linked to the rate of mechanical stretching of the interfacial area. This mechanism has been widely analyzed [START_REF] Ottino | The kinematics of mixing: stretching, chaos and transport[END_REF][START_REF] Roberts | The simulation of stretch rates for the quantitative prediction and mapping of mixing within a channel flow[END_REF]. This approach was revisited in Mackley and Neves [START_REF] Mackley | The quantitative description of fluid mixing using Lagrangian-and concentration-based numerical approaches[END_REF] using a Lagrangian and concentration based numerical approach.

The approaches quantifying the mixing from the increase of the material interface have two important drawbacks: (i) the first one is related to the difficulty of introducing other additional physics, as the one related to the surface tension, and (ii) sometimes the microstructure description needs other information (morphology, characteristic length, shape and orientation of the discrete phase -inclusions-, ...) that the interface evolution cannot provide. The evolution of the morphology can then be modeled by predicting the change in the local morphological measure due to this velocity or deformation gradient. We will call such an approach, which treats some local characteristic morphological measure as a field variable, a micromixing analysis. From now on we focus on incompressible fluid flows and passive mixing.

The micromixing approach has been used extensively for modeling passive mixing, where interfacial energy is negligible and the two phases have identical viscosities. In passive mixing the global velocity field can be found independently of the microstructure and then used to evolve the mixture structure, described with some area tensor, very rich from the morphological and microstructural points of view. See [START_REF] Wetzel | Area tensors for modeling microstructure during laminar liquid-liquid mixing[END_REF] for an excellent description of this approach, in which the genesis of this kind of approaches is retraced precisely. The main difficulties noticed in that work [START_REF] Wetzel | Area tensors for modeling microstructure during laminar liquid-liquid mixing[END_REF] were related to the necessity of introducing a closure relation in order to derive an evolution equation for the second order area tensor. Different closure relations have been proposed, being their impact, in general, unpredictable.

A very effective way, widely used, of quantifying the mixing rate can be derived from the consideration of only one of the basic mechanisms of mixing processes: the increase of the material interface due to the flow kinematics in absence of surface tension and molecular diffusion. Thus, the higher is the interface per unit of volume, the more effective is the mixing. Because it is impossible from a practical point of view the accurate tracking of fluid interfaces, coarser descriptions have been successfully proposed.

For this purpose we can introduce the interface distribution function (x, t, n) given, at each point in the physical domain x ∈ and for any time t, the specific surface (i.e. interface per unit of volume) oriented in the direction n [START_REF] Doi | Dynamics and rheology of complex interfaces[END_REF]. If we assume that the interface evolution only depends on the flow kinematics, the evolution of (x, t, n) is given by

d dt = - ∂ ∂n • ( ṅ ) -(∇v : (n ⊗ n)) (1)
where d dt denotes the material derivative, i.e. d dt =

∂ ∂t + v • ∇ , being v the velocity field, ⊗ denotes the tensor product and : the tensor product twice contracted. The time evolution of the interface orientation ṅ induced by the flow velocity field v, is given by

ṅ = -(∇v) T n + (∇v : (n ⊗ n))n (2)
On the other hand the flow kinematics induces also a variation of the interface area. Thus, a surface element dS whose orientation is defined by the unit vector n evolves according to:

Ṡ = dS dt = -(∇v : (n ⊗ n)) dS (3)
The main difficulty related to such an approach making use of the interface distribution function (x, t, n) lies in its multidimensional character because the interface distribution function depends on the physical coordinates x, the time, and the conformation coordinates related to the interface orientation n.

To circumvent the curse of dimensionality related to such kinetic theory approaches coarser descriptions can be derived by introducing the so-called area tensor, that as described below contains very valuable information on the fluid morphology. Let be the domain in which the flow problem is defined. Points in will be referred by x which consists in a vector in 2D or 3D. In the numerical examples presented in this work only 2D models are considered because in such case reference solutions can be computed. In order to quantify the morphology at any point x ∈ a microscopic representative volume V(x) is considered centered at that point. This volume is small enough with respect to the macroscopic scale (related to the variation of the velocity field v in ) but large enough with respect to the characteristic size of the microstructure. Let S V (x, t) be the interface within V(x) at time t. The area tensor A(x, t) is then defined as:

A(x, t) = 1 V(x) S(x,t) n ⊗ n dS (4)
The equation governing the evolution of the area tensor can be obtained by taking time derivative of Eq. 4 and then using Eqs. 2 and 3. After some manipulations it results:

Ȧ = -(∇v) T A -A∇v + ∇v : A (5)
which involves a fourth order tensor

A A = 1 V(x) S(x,t) n ⊗ n ⊗ n ⊗ n dS (6)
and where the dependence of A, A and v on the space and time coordinates, x and t respectively, is omitted for the sake of clarity. The area tensor A just defined is symmetric and has following appealing properties:

1. The first property concerns its trace (sum of the diagonal components of A) that we symbolize by Tr(A) that taking into account the normality of n results:

Tr (A(x, t)) = 1 V(x) S(x,t) (n 2 1 + n 2 2 + n 2 3 ) dS = 1 V(x) S(x,t) dS = S(x, t) V(x) = S V (x, t) (7)
Thus, the trace of the area tensor gives the specific surface whose maximization is currently searched in mixing processes. 2. If we define the volume fraction of the disperse phase as φ, assumed uniformly distributed within , then we can define a characteristic length L(x, t) of the microstructure at point x from the volume of the discrete phase

V d (x) = φ • V(x): L(x, t) = V d (x) S(x, t) = φ • V(x) S(x, t) = φ S V (x, t) = φ Tr(A (x, t)) (8)
3. In passive mixing the strain rate varies linearly at the scale of V(x), so that an initially spherical discrete domain will deform into an ellipsoid for any state of strain. This characteristic ellipsoid provides a convenient way to interpret the area tensor.

The microstructure shape and orientation can be easily deduced from the normalized area tensor à defined as

Ã(x, t) = A(x, t) Tr(A (x, t)) (9)
and whose eigenvalues allow computing the length of the ellipsoid axes, being their orientation given by the associated eigenvectors. See [START_REF] Wetzel | Area tensors for modeling microstructure during laminar liquid-liquid mixing[END_REF] for more details concerning the relation between that ellipsoid and the area tensor.

It is easy to verify that the interface distribution function and the area tensor descriptions can be linked by taking into account the expression

A(x, t) = B (n ⊗ n) (x, t, n) dn ( 10
)
where B represents the surface of the unit sphere. The solution of Eq. 5 requires the introduction of a closure relation for expressing the fourth order area tensor A as a function of the second order one A. Different closure relations have been proposed, being the simplest ones, the quadratic and the orthotropic [START_REF] Wetzel | Area tensors for modeling microstructure during laminar liquid-liquid mixing[END_REF].

The quadratic closure, initially proposed by Doi and Ohta, writes:

A = 1 S V (A ⊗ A) (11)
The impact of the quadratic closure relation was deeply analyzed in [START_REF] Chinesta | Microstructure evolution during liquid-liquid laminar mixing: a kinetic theory approach[END_REF] where the solution computed by using Eq. 5 with the quadratic closure (Eq. 11) was compared to the area tensor obtained from the interface distribution function (x, t, n) according to Eq. 10, where the distribution function (x, t, n) was computed by integrating Eq. 1 without the necessity of introducing any closure.

The orthotropic relation proposed by [START_REF] Cintra | Orthotropic closure approximations for flow-induced fiber orientation[END_REF] postulate a different dependence of A on A [START_REF] Wetzel | Area tensors for modeling microstructure during laminar liquid-liquid mixing[END_REF].

Until now, the only mechanism affecting the interface evolution is the flow kinematics. However one could expect that other phenomena could participate actively in the interface evolution, as for example the surface tension, the droplets coalescence or the droplets break-up.

Enriched modeling

One could expect that in absence of flow the microstructure evolves toward an isotropic state, i.e. spherical droplets, and that the interface area decreases because the droplets coalescence. These mechanisms are expected to be dependent on the matrix viscosity η and the surface tension Ŵ.

Thus, Doi and Ohta considered the evolution of the specific surface, i.e. the trace of the area tensor, given by ṠV = -r 1 • S V [START_REF] Roberts | The simulation of stretch rates for the quantitative prediction and mapping of mixing within a channel flow[END_REF] whereas the deviatoric part of the area tensor tends to vanish in absence of flow

d dt A -Tr(A) 3 I S V = -r 2 • A -Tr(A) 3 I S V ( 13
)
where r 1 and r 2 represents two relaxation rates.

A dimensional analysis leads to:

r 1 = c 1 Ŵ • S V η (14) r 2 = c 2 Ŵ • S V η (15)
where c 1 and c 2 are two positive dimensionless coefficients.

Lee and Park considered a third mechanism related to the droplet break-up whose effect on the evolution of the specific surface was assumed depending on the deviatoric part of the area tensor according to [START_REF] Lee | Rheology and dynamics of immiscible polymer blends[END_REF]:

ṠV = -c 3 Ŵ η A - Tr(A) 3 I : A - Tr(A) 3 I (16)
Thus the effects of these extra-mechanics write:

d dt S V extra = -c 1 ŴS 2 V η -c 3 Ŵ η A - Tr(A) 3 I : A - Tr(A) 3 I (17)
and

d dt A -Tr(A) 3 I S V extra = -c 2 ŴS V η A -Tr(A) 3 I S V ( 18
)
where c 1 , c 2 and c 3 are three positive dimensionless coefficients which could depend on the volume fraction φ. By introducing the new dimensionless coefficients λ = c 1 + c 2 , µ = c 1 /(c 1 + c 2 ) and ν = c 3 /(c 1 + c 2 ), the equation governing the evolution of the area tensor writes:

dA dt = -(∇v) T A -A∇v + ∇v : A -λ Ŵ η S V A - S V 3 I -λµ Ŵ η S V 2 I 3 -λν Ŵ η A : A - S V 2 3 A S V (19)
Evaluating the impact of closure relations Chinesta and Mackley [START_REF] Chinesta | Microstructure evolution during liquid-liquid laminar mixing: a kinetic theory approach[END_REF] analyzed the impact of the quadratic closure in the evolution of the area tensor in absence of droplets break-up. In this section we are performing a similar analysis when shape relaxation, coalescence and break-up mechanisms are present. Inspired from the procedure described in [START_REF] Chinesta | Microstructure evolution during liquid-liquid laminar mixing: a kinetic theory approach[END_REF], it is easy to prove that the equation governing the evolution of the interface distribution function reads:

d dt = - ∂ ∂n ( ṅ ) -(∇v : (n ⊗ n)) - ∂ ∂n D( ) ∂ ∂n + F( ) (20)
where the diffusion coefficient and the source terms are given by:

D( ) = λ Ŵ η 1 6 S V + ν S V A : A - S V 2 3 (21)
and

F( ) = -λ Ŵ η 1 4π µS V 2 + ν A : A - S V 2 3 . (22)
with

A = B (n ⊗ n) dn (23)
The solution of Eqs. 20-23 allows computing the interface distribution function, and then the second order and the fourth order area tensor, A and A respectively, without the introduction of any closure relation. The price to be paid is the solution of a model defined in a multidimensional space. In general complex flows the solution of these equations is quite difficult despite the recent progresses accomplished in numerical strategies able to circumvent or at least alleviate the multidimensionality issue [START_REF] Ammar | A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modeling of complex fluids[END_REF][START_REF] Ammar | A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modeling of complex fluids. Part II: transient simulation using spacetime separated representations[END_REF][START_REF] Mokdad | A fully deterministic micro-macro simulation of complex flows involving reversible network fluid models[END_REF].

In the case of simple rheometric flows we can assume a homogeneous kinematics and then ignore the dependence of all the variables on the physical coordinates. In this case, the resulting kinetic theory model making use of the interface distribution function can be solved easily because it only needs an appropriate discretization of the surface of the unit sphere.

Numerical results

Area tensor evolution equation: Dimensionless form

The area tensor equation can be rewritten in the dimensionless form by introducing the dimensionless capillary number Ca defined by

Ca = η γ ŴS init V ( 24
)
where η is the viscosity, S init V the initial specific surface, Ŵ the surface tension and γ the equivalent strain rate defined by γ = √ ∇v : ∇v. The dimensionless form is obtained by affecting both the orientation tensor and the specific surface by the initial specific surface and time by the inverse of the shear rate. Thus, the dimensionless form of the equation governing the area tensor evolution reads Different values of the parameters λ, µ and ν where considered for analyzing their effect on the time evolution of the specific area.

dA dt = -(∇v) T A -A∇v + ∇v : A + -λµ 1 Ca S 2 V I 3 -λ 1 Ca S V A - S V 3 I -λν 1 Ca A : A - S 2 V 3 A S V ( 25
Figures 1 and2 depict the specific surface time evolution in the simple shear flow for Ca = 1 and Ca = 10 respectively, whereas Figs. 3 and4 show similar results in the case of the bi-elongational flow. In these figures we can appreciate the impact of the closure relations. Despite the fact that the quadratic closure relation introduces in some cases an appreciable deviation with respect to the reference solution computed within the kinetic theory framework, its implementation is the most efficient from the computational time viewpoint. In these figures we can also notice the strong sensibility of the solution to the model parameters Ca, λ, µ and ν.

Simulating complex flows

We consider the mixing device whose geometry is depicted in Fig. 5. The area tensor will be represented by using the ellipse associated to the rotation of π 2 of tensor A. Such representation is closer to the real microstructure morphology. The ellipse axes are given by the eigenvectors of tensor A rotated of π 2 whereas their lengths are related to the associated eigenvalues. We are also representing using a color map the value of the inverse of the specific surface that is a key parameter for quantifying the mixing. Initially, the droplets should be represented as spheres, i.e. area tensor is diagonal, all the components having the same value.

In the simulations that we address in this section we consider a coupling between the flow kinematics and the microstructure evolution:

-Flow kinematics

∇ • σ = 0 (28) ∇ • v = 0 (29) σ = -pI + 2ηD - Ŵ η A -S V I 3 -Microstructure evolution dA dt = -(∇v) T A -A∇v + ∇v : A + -λµ 1 Ca S 2 V I 3 -λ 1 Ca S V A - S V 3 I -λν 1 Ca A : A - S 2 V 3 A S V ( 30
)
where a quadratic closure relation was considered, because it represents a good compromise between accuracy and computational efficiency. Finally, Fig. 6 compares the streamlines for two different microstructural morphologies: the initial one perfectly isotropic and the one induced after a complete turn. It can be noticed that in the flow here considered the difference is not too significative.

Conclusion

The use of the area tensor allows describing the morphology evolution during mixing processes at the mi-croscopic scale and is very convenient in term of microstructural description (specific interface, shape, orientation, ...).

In this paper we considered a model taking into account the flow effects, surface tension, droplets coalescence and droplets break-up. The equation governing the evolution of the area tensor in that model involves a fourth order area tensor that must be written as a function of the second order orientation tensor by using an appropriate closure relation that introduces in general an unpredictable error. To evaluate the effects of a quadratic or an orthotropic closure relation in the area tensor evolution first we derive a kinetic theory model governing the evolution of the so-called interface distribution function that does not involve any closure relation, from which the different area tensors can be calculated. Then, we compared in a simple shear flow the reference solution (computed from the interface distribution function) and the ones obtained by using both closure relations. Both closure relations seem to be good, and the quadratic one was retained because its simplicity. Finally, the model was used for simulating complex flows in mixing devices.

)

  Evaluating the closure relations in a simple shear flowIn what follows we analyzed the evolution of the specific surface S V for two different capillary numbers Ca = 1 and Ca = 10 in the simple shear and bielongational flows defined respectively by:

1 Fig. 1 1 Fig. 2 10 0 1 Fig. 3 Fig. 4 10 (Fig. 5

 111210134105 Fig. 1 Specific surface evolution in the simple shear flow with Ca = 1

Fig. 6

 6 Fig. 6 Influence of the microstructure evolution on the flow streamlines


Michael Kruse

Introducing Molly: Distributed Memory Parallelization with LLVM

Keywords: Distributed memory, LLVM, Polyhedral model, HPC 1. Introduction

Programming for distributed memory machines has always been a tedious task, but necessary because compilers have not been sufficiently able to optimize for such machines themselves. Molly is an extension to the LLVM compiler toolchain that is able to distribute and reorganize workload and data if the program is organized in statically determined loop control-flows. These are represented as polyhedral integer-point sets that allow program transformations applied on them. Memory distribution and layout can be declared by the programmer as needed and the necessary asynchronous MPI communication is generated automatically. The primary motivation is to run Lattice QCD simulations on IBM Blue Gene/Q supercomputers, but since the implementation is not yet completed, this paper shows the capabilities on Conway's Game of Life.

for the IBM Blue Gene/Q supercomputer. The data layout was changed multiple times in order to find out which the fastest one is. Every time large portions had to be rewritten such that the overhead remains low for a fair comparison. The different versions were conceptually similar, which motivated to let the compiler do the main work using already existing techniques that represent a program as a polyhedron [START_REF] Feautrier | Some efficient solutions to the affine scheduling problem. I. One-dimensional time[END_REF][START_REF] Fautrier | Some efficient solutions to the affine scheduling problem. Part II. Multidimensional time[END_REF][START_REF] Darte | Scheduling and automatic Parallelization[END_REF].

The core of most Lattice QCD solvers is an 8-point stencil in 4 dimensions with 1320 floating-point operations and 2688 bytes to read per site update in double precision. It is the most optimized part and any overhead may significantly slow down the solver compared to hand-optimized code. A compiler-optimized version has to be similar to the equivalent manually optimized code in order to be useful or physicists will continue optimizing code by hand. Dynamic techniques at runtime such as task scheduling are therefore off the table.

Usually just the instructions are represented and reordered in the polyhedral model, although also the elements of an array can be modeled this way. The alpha language does this to compact array storage [START_REF] Mauras | Scheduling affine parameterized recurrences by means of variable dependent timing functions[END_REF][START_REF] Quiller | Optimizing memory usage in the polyhedral model[END_REF]. It is however rarely used to just change the location of data elements although it is a powerful representation. [START_REF] Größlinger | Precise management of scratchpad memories for localising array accesses in scientific codes[END_REF] suggests this for scratchpad memories, this paper for homogeneous distributed memory machines (DMM).

This work extends the Clang+LLVM [START_REF] Lattner | LLVM: A compilation framework for lifelong program analysis & transformation[END_REF][START_REF] Lattner | Clang: A C Language Family Frontend for LLVM[END_REF] compiler toolchain with a new pass called Molly. The choice fell to LLVM because [START_REF] Jansen | tmLQCD: A Program suite to simulate Wilson Twisted mass Lattice QCD[END_REF] it is open source, (2) relatively easily extendable, (3) promotes link-time optimization (4) it already contains many necessary components, especially Polly and (5) can act as a Just-In-Time (JIT) compiler.

Point (3) is important because reordering storage is an inherently global process. Different components cannot be translated inependently if they do not agree on the data layout of process-wide data. Polly [START_REF] Grosser | Polly -Performing polyhedral optimizations on a low-level intermediate representation[END_REF] is the polyhedral optimization subproject of LLVM. It provides detection of static control parts (SCoPs), is able to optimize them and generates code again. Molly intercepts Polly to modify the SCoPs before code re-generation.

The fifth point is relevant because the current implementation of Molly requires the size of distributed arrays and the geometry of compute nodes to be known at compile time. This is because blockcyclic (or just any interesting) data distribution is not a linear function. If N is the length of an array, P the number of processors, i an element of the array array, p a processor number, and l the index of element i on processor i, then the relation for block-distribution is i = p(N/P ) + l. The term p(N/P ) is not of the form "constant times variable" required to be an affine term. p naturally varies in an MPI/SPMD program, therefore N/P , the block size, must be constant. Future versions may precompute p(N/P ) and parametrize the SCoP using this newly introduced variable, such that at least block-distributions do not need to predefined the geometry of the cluster machine.

The Integer Set Library (ISL) [START_REF] Verdoolaege | isl: An integer set library for the polyhedral model[END_REF] is, despite its name, a library for sets of vectors of integers and rationals represented as unions of Z-polyhedra. It is already used by Polly and Molly uses it as well for indexsets.

This article is organized as follows. Section 2 gives an overview on other solutions on the same problem. Section 3 introduces the toolchain and what parts had to be modified for Molly. In Section 4 the main part about how generating communication code is presented and Section 5 shows this on the example of Conway's Game of Life. Experimental execution results are presented in section 6 and possible further developments discussed on Section 7. Finally, Section 8 summarizes this work.

Related Work

In the past, efforts to shift the responsibility of parallelization and data distribution to the compiler were not very successful. Probably the most noteworthy project to relief the programmer of the burden of explicit parallelization is High-Performance Fortran (HPF). Today, HPF plays a minor role because of many reasons [START_REF] Kennedy | The rise and fall of High Performance Fortran: an historical object lesson[END_REF]. One of them is that compilers promised to optimize code, but could only do so for specifically structured code. Every compiler could only optimize a different set of code patterns.

Probably the most similar work is [START_REF] Adve | Using integer sets for data-parallel program analysis and optimization[END_REF]. The authors also build polyhedral subsets of array indexes that have to be transferred to different nodes. They call a "communication event" is called a chunk in this paper. However, their set operations are very different and also does more advanced techniques like index set splitting and merge of communication events. Their tool dHPF is a HPF to Fortran77 with MPI source-to-source compiler and uses the Omega [START_REF] Kelly | The Omega calculator and library, version 1.1.0[END_REF] library for the set operations. It only allows blockcyclic data distribution and uses the owner-computes policy for statements. More work of this kind are the Last Write Tree [START_REF] Amarasinghe | Communication optimization and code generation for distributed memory machines[END_REF] and [START_REF] Dathathri | Generating efficient data movement code for heterogeneous architectures with distributed-memory[END_REF].

A different kind of optimization partitions the work into tiles which are then distributed to the nodes [START_REF] Bondhugula | Automatic Distributed-Memory Parallelization and Code Generation using the Polyhedral Framework[END_REF][START_REF] Classen | Automatic code generation for distributed memory architectures in the polytope model[END_REF][START_REF] Yuki | Memory Allocations for Tiled Uniform Dependence Programs[END_REF]. The tiles size is fixed at compile-time, but the number of nodes doesn't need to be fixed until runtime. Therefore most of the work uses a notion of "virtual" processors. Virtual processors are mapped to physical ones at runtime.

A novel technique has been presented in [START_REF] Ravishankar | Code generation for parallel execution of a class of irregular loops on distributed memory systems[END_REF]. It combines syntactical code analysis with a dynamic part that allows non-static code behavior such as indirect array accesses. These can normally not be represented in a SCoP.

Again another approach is to invent a language that has highlevel operators on multi-dimensional arrays. ZPL [START_REF] Snyder | A Programming Guide to ZPL[END_REF] belongs to this category. Instead of "accessing the element to the left" there is a shift operator that moves the entire field one index to the right. Such languages require rethinking of algorithms, but high-level operators can be pre-implemented and algebraically optimized efficiently.

The middle way are extensions to existing languages. To name a few, there is OpenACC, OmpSs, C++ AMP, Universal Parallel C (UPC), etc. The first three are language extensions that targeting accelerators like GPUs to offload computation from the main CPU. The programmer has to state explicitly the data region that must be available on the device memory. UPC, a representative of languages using PGAS (Partitioned Global Address Space), create an address space window for non-local data. Code flow still has to be explicitly organized and synchronized.

The Toolchain

The reduced variant of Conway's Game of Life shown in Figure 1 illustrates a program that Molly can optimize. The original Game of Life is a 9-point stencils whereas this variant has only 5 points. The example does 100 iterations and copies back the last iteration's result because pointer swaps cannot (yet) be handled in SCoPs. Such molly::array-typed variables that are distributed are called fields in this paper. This comes from physics' jargong of a quantity for every point in spacetime. They are a typical application for such arrays, notably in lattice QCD.

The Clang Part

Clang [START_REF] Lattner | Clang: A C Language Family Frontend for LLVM[END_REF] is the C, C++ and Objective-C frontend that compiles source files of those languages into LLVM [START_REF] Lattner | LLVM: A compilation framework for lifelong program analysis & transformation[END_REF] intermediate representation, also called LLVM IR. It has to be modified in order to pass information on distributed arrays to the optimizer as metadata.

First, there is a new header file molly.h which declares the molly::array variadic template class. It represents a distributed array of arbitrary size and dimensions. Its semantics are different than the standard C++ arrays. Pointer arithmetic beyond the current element are invalid. The compiler may change the elements' order, the value might be outdated or even not stored on the local node.

The header file uses an overloaded index operator in order to make field accesses like back[x][y] = living resemble bool *ptr = __builtin_molly_ptr(&back, x, y); *ptr = living; builtin molly ptr is a compiler-builtin and returns a pointer. The value of the pointer is meaningless if the value has no location on the current node. Later passes need to fix-up for this.

Clang also generates metadata for every distributed array. The information include what the element type is, as well as the shape of the array. The template arguments are not kept when translated to the respective LLVM type.

In principle, any language frontend can generate such kind of accesses and forward it to LLVM, but for the current project only C++ is supported.

The LLVM Part

The intermediate representation of the previous snippet that is passed to LLVM is shown in Figure 2. The newly introduced LLVM intrinsic llvm.molly.ptr semantically resembles the %1 = call i8* llvm.molly.ptr(%"class.molly::array"* @back, i64 %x, i64 %y) store i8 %living, i8* %1 GetElementPtr instruction: it returns the pointer to an array's element, but logically allows accessing remote memory.

Because pointer arithmetic is disallowed on values returned by llvm.molly.ptr, only load and store instructions can use it. Whenever the Molly pass encounters a load/store of a pointer returned by llvm.molly.ptr it may look up the accessed field and the coordinates to know what element is accessed.

Molly is a module-level pass, therefore it has access to all molly::array<> variables and all the location they are used. Hence, it is possible to align value distributions of the fields, but is currently not implemented. A simple block-distribution is used and instructions are executed on the node that owns the value they compute (owner computes rule).

The next step is to create an isolated basic block for every field access. The splits must happen because the data source domains where either local data of remote data is accessed do not match. For instance, the data for front[x-1][y] is local to node (0, 0) on iteration (x, y) = (512, 1), but on node (1, 0) for the access front[x][y+1] although both accesses must be executed on the same node because they both access the scalar neighbors.

Data between the statement is exchanged through variables on the stack, for instance the variable neighbors. There is no internode transfer generated for such non-field scalar variables and therefore must always be executed on the same node.

Polly: SCoP-Detection

At this stage Polly detects whether a loop structure is SCoPcompatible and if so, identifies the statements and their natural scheduling parameters (iteration domain and relative ordering). The isolation from the previous step ensures that every field access gets its own statement and can therefore scheduled independent of the rest.

Next, Molly introduces two virtual statements that are executed before and after all the statements in the SCoP, called prologue and epilogue. The prologue virtually writes everything and the epilogue reads every data location referred to in the SCoP. They are use to compute the data flow into and out of the SCoP.

The data dependencies between the statements are computed by ISL. Molly requires the flow dependencies (read-after-write) of data in fields and local memory independently. Other types of dependencies (write-after-write and write-after-read) are irrelevant because Molly does not change the ordering of instructions, but needs to know where data is used and generated.

ISL determines direct dependencies only, i.e. for every value that is read it is able to determine the unique statement instance that computed that value, as long as all array index expressions are affine.

Molly

Then, Molly has to decide on which node a statement is executed. By default, statements are executed preferably on the same node as one of the statements that produces its input. Statements that produce data that needs to be written back at the SCoP's end (i.e. has a dependency to the epilogue) is computed on that data's home location. Analogous for statements that read data from before the statements, i.e. dependencies from the prologue.

Statements that produce data written to a node-local variable (i.e. non-field flow dependency) need to be executed on at least all the nodes that potentially consume it.

The next section also belongs to the Molly-part, but explains the core message passing optimization and therefore has its own section.

Communication Code Generation

The main part is the generation of code that transfers the data between the nodes. We start off explaining some terminology used in polyhedral optimization and notation used in this paper.

Application Modeling

Molly needs information about six different spaces and the relation between them in order to transform a program. The spaces are

• Ω, the set of all statement in the SCoP. Typical names for statements are R, S1, S2, G, C etc.

• P ⊆ Z N P , the set of computing nodes or processors; also called the cluster space. A coordinate is typically named (a, b) or p.

• DS ⊆ Z N S , the domain of the statement S ∈ Ω; also called the iteration space. Typical variable names are i, j. (i, j) can also denote a single coordinate.

• F, the set of fields. The only specific fields on this paper are f ront and back, the arrays from Figure 1.

• Ie ⊆ Z Ne , the set of addressable indexes of a field e ∈ F, therefore called the indexset space. An index is typically named k.

• T ⊆ Z N T , the scheduled time of the statements in the current SCoP, also called the scatter space.

The NP,S,e,T describe the number of dimensions of the spaces. A tuple (S, i) ⊆ Ω × DS is called an instance of statement S. A tuple (S, i, p) ⊆ Ω × DS × P is called an execution of instance (S, i) on node p.

The previous analyses of a SCoP by Polly and Molly resulted in the following functions and relations between those spaces:

• θS : DS → T , the schedule of every statement representing the time of execution in a sequential program.

• λS : DS → 2 F ×I , the field elements a statement accesses, usually at most one, with the exception of the prologue and epilogue.

• (G, iG) δ (C, iC ) is the flow-dependence relation in which the consumer instance (C, iC ) reads a value computed by the producer/generator instance (G, iG).

• πe∈F : Ie → 2 P , the node that owns and stores a field element, called the home location. Multiple nodes may be responsible for the same element in which case all of them must store the same value.

• πS∈Ω : DS → 2 P , on which node a statement will be executed.

A single statement can also be executed on several nodes. This is called the where-mapping. An execution (S, i, p) is therefore restricted to the nodes p ∈ πS.

For any two statements G and C let (G, iG) δ (C, iC ) the data flow dependence between a data-generating statement G and a con-sumer C. They both access the field element λG( iG) = λC( iC ). In case the statement instances are executed on different nodes, some communication code must be generated between them.

To avoid transferring single values with immense overhead values are grouped into chunks, such that all data of one chunk can be written to the same buffer (per destination), then sent and received. Two values cannot be in the same chunk iff there is a direct or indirect dependency between the statement that produce them. Definition 1 (Chunking Function). A chunking function is an projection ϕ : (Ω, D) → (Ω, D) of statement instances to representative instances. By convention ϕ is idempotent, i.e. any representative instance projects to itself. The equivalence classes

[ϕ(R, j)] = {(S, i) | ϕ(R, j) = (S, i)} are called chunks.
A chunking function is a mean to identify chunks. One can think of the representative instance as the first statement of the chunk, although an order may not have been established yet. In fact, any image space can be used, but any other choice would be more arbitrary.

Definition 2 (Valid Chunking Function). Let δ * be the transitive closure of flow dependencies between all statements. A chunking function ϕ is valid if it does not violate any dependencies in δ * , i.e. ∃(S, i, S, i) ∈ ϕ(δ * ), given that δ * was irreflexive before.

The objective is to find a valid chunking function with as few elements as possible and preferably whose chunks are about the same size. The worst chunking, but always possible function is the identity function. This yields to every value transferred individually as mentioned before.

If the chunking function as well as the dependencies are (piecewise) affine, finding a chunking function with minimal chunks involves finding a solution of an inequality system.

Let (G, iG) → (C, iC ) be a dependency of δ * . iG = iC because otherwise the instance requires data that itself produces which is invalid. Let ϕ {G,C} ( i) = α0io + α1i1 . . . a possible chunking function with unknown α0,.... For validity, we have to find α0,..., s.t. ϕ( iG) = ϕ( iC ) for all flow dependencies. The goal is a function ϕ whose image spans minimally per dimension.

This problem is analogous to the scheduling problem to find a schedule with minimal span and therefore maximal parallelism. One can think of the chunking function collapsing a set of statement instances that are executed in parallel, but here their data is received at the same time. An algorithm is suggested in [START_REF] Bastoul | Improving data locality by chunking[END_REF].

Molly uses a much simpler algorithm which exploits that loops are typically already written in a way that allow chunking of innermost loops. The algorithm finds the outermost loop such that (1) all consumer instances of one chunk are executed after all of the generator instances and (2) the resulting chunking function is valid. ( 2) implies (1) because the consumer is always dependent on the generator and therefore is always scheduled later, but it is much faster to check and covers most practical cases. Also, this gives the direct scatter locations where to insert the send and receive statements and therefore guarantees schedulability. The algorithm is shown as Algorithm 1.

Once the chunking function id known, Molly can insert calls that trigger communication and redirect any accesses to remote data to the communication buffers.

First, Molly builds a giant relation between the representative instance of the chunk, generator execution, consumer execution, the field element they access:

T ransf ers := ((R, iR), (G, ig, pg), (C, iC , pC), (F, k f )) | ig ∈ DG, pg ∈ πG(ig), ic ∈ DC , pc ∈ πC (ig), (G, ig)δ(C, ic), λC (G, ig) = (F, k f ) = λC( ig), (R, iR) = ϕC ( iC )
A generator instance can have multiple executions on different nodes and therefore one has to be picked, preferably the one that executed on the same node as the consumer if there is one. Molly's algorithm for this is to pick the lexicographical minimal. The result, T ransf ers ′ , has a unique tuple (G, ig, pg) for every tuple (C, iC , pC), (F, k f ), and therefore can be defined as a function Source : (C, iC , pC, F, k f ) → (G, ig, pg).

To construct the communication events, the relation elements in T ransf ers ′ are grouped by the chunk they belong to (Equation 1).

For any such chunk Molly creates two static variables, one for the sending send and another for the receiving side, and six statements which replace the consumer executions.

1. A call to molly send wait, which waits for the previous chunk to finish sending data if still active and returns a pointer to the memory buffer into which the data of the current chunk must be written. The scatter function on node pG to node pC is Equation 2, i.e. one timestep before the first generator execution. To avoid overlapping with an already existing instance at this timestep one can simply multiply all the scatter function by 2 or more beforehand.

2. A new generator statement. Instead of writing to local memory, it writes the buffer returned by the previous statement at an index that is computed from the hull of field indices that are to be transferred between the sending and receiving node. The scatter function is the same as the original generator statement.

3. A call to molly send that starts the transfer. Its scatter is one step after the last generator execution of this chunk.

4.

A call to molly recv wait on every receiving node iC. Its scatter function is Equation 3, one step before the first consume execution on this node. This runtime function returns a pointer to the buffer that has been received.

5.

A replacement consumer statement which delegates the original return value of molly ptr to the receiving buffer at the index

Chunk : (R, iR) → ((G, ig, pg), (C, iC , pC), (F, k f )) ((R, iR), (G, ig, pg), (C, iC , pC), (F, k f )) ∈ T ransf ers ′ Chunks := Chunk(R, iR) Chunk(R, iR) = ∅ (1) θsend wait(R, iR) = min lex ({ θ(G, iG)|((G, ig, pg), (C, iC , pC), (F, k f )) ∈ T ransf ers ′ }) -(0 . . . , 1) (2)
θrecv wait(R, iR) = min lex ({ θ(C, iC )|((G, ig, pg), (C, iC , pC ), (F, k f )) ∈ T ransf ers ′ }) -(0 . . . , 1) (3)
θrecv(R, iR) = max lex ({ θ(C, iC )|((G, ig, pg), (C, iC , pC ), (F, k f )) ∈ T ransf ers ′ }) + (0 . . . , 1) (4)
computed the same way as the sendbuffer of the generator statement.

6. A call to molly recv that releases the receive-buffer and resets the communication object in order to be prepared for the next chunk. Its scatter function on node is Equation 3, that is, one timestep behind the last replacement consumer.

In addition, Molly generates calls initialization-and releasefunctions for every chunk and transfer source/destination . They setup (respectively free) the target node coordinate, the buffer size and an unique identifier (tag) for every communication event in advance. For the receiving nodes, it also resets the buffer to ready state.

Finally, all the original generator statements are removed. All the consumers now read from the communication buffers such that the original writes are unused. As a result, no calls to llvm.molly.ptr remain active in this SCoP.

Dataflows from the prologue and to the epilogue are special cases of the more general flow between arbitrary statements. The difference is that the reading part reads from, respectively writes to the local storage of the field access per node. The prologue and epilogue are singleton statements (per node) that touch all the data. Therefore, there is just a single chunk that transfers data to, respectively out of the SCoP. Flows from the prologure to the epilogue can be ignored.

Since Polly does not understand placements of statement instances (π) one has to reduce them to something more generic. The coordinate of the currently executing node is added as a parameter of the SCoP (the same way runtime parameters like the upper bound of a for-loop). Then, the iteration relation i, which implicitly uses these parameters, is modified such that the statement instances are removed that

Polly: Code Optimization and Generation

Once the communication statements and iteration domains are fixed, the set of statements are given to Polly again. Polly may apply additional optimizations that reorder statements by applying different scatter functions (θ), vectorize and parallelize using OpenMP. There are two different optimization engines that Polly supports: PLuTo [START_REF] Bondhugula | A practical automatic polyhedral parallelizer and locality optimizer[END_REF] and ISL [START_REF] Verdoolaege | isl: An integer set library for the polyhedral model[END_REF].

Polly's final part is to generate IR code again from the statement instances. Again, Polly supports two engines: Cloog [START_REF] Bastoul | Code generation in the polyhedral model is easier than you think[END_REF] and ISL [START_REF] Verdoolaege | isl: An integer set library for the polyhedral model[END_REF]. The original code will be deactivated and finally removed by the dead code elimination pass.

ISL is required by Polly in any case and can do optimization and code generation, the other two engines are therefore optional.

Finalization

A few smaller task remain to be done. There can be leftover llvm.molly.ptr remaining in the code that were not inside a SCoP because they potentially access non-local data and they are unknown to the LLVM machine code generator 1 . Accesses to its returned pointer are replaced by the runtime functions with the names molly value load and molly value store that handle single value transfers in a synchronous manner. This is slow, but necessary for correctness.

Moreover, a set of functions for accessing the local data storage is being generated. The data ordering is just known to the compiler, therefore the runtime needs a way to find the correct index inside local memory. They are needed, for instance, to implement the previous functions that transfer single values between nodes.

Finally the main function is wrapped by another runtime function which now serves as the entry point. The runtime main initializes the communication API and calls the transfer buffer initialization functions, calls the application's main, and releases any resources before exiting the application. As additional parameters it receives the geometry of nodes the application has been compiled for. It must abort if it doesn't match the geometry it is running an since correctness is not guaranteed anymore.

The runtime needs to call the application's main function by itself because some communication APIs may change the application's command line parameters. Noteably MPI Init is allowed to do this if the implementation submits the command line parameters by some other way than by calling the executable on the target node with said parameters in advance.

MollyRT

MollyRT is the runtime library for applications compiled with Molly. It provides implementations of functions that are called during the code generation as mentioned in the previous section.

The current library's only backend is the Message Passing Interface (MPI) which is available on virtually every cluster supercomputer these days. But other interfaces are thought of as well, especially the close-to-the-metal System Programming Interface (SPI) of the Blue Gene/Q. The MPI overhead on this machine is immense. Other possible backends are possible of as well: Transfers between threads/processes of a shared memory machine, SHMEM, IBM PAMI, PGAS, etc.

The mapping between the MollyRT functions and the MPI functions is obvious in ost cases. Communication buffers setup requests using MPI persistent transfers upon their initialization. The wait-functions call MPI Wait on these. molly send and molly recv issue an MPI Start.

Example: Conway's Game of Life

The first part of the Molly passes is the isolation of basic blocks that contain field accesses. Originally, the Game of Life code in Figure 1 has two blocks inside the loops that do not belong to a for-loop. We can ignore the return statement as is does not contain a field access nor does it belong to a SCoP. Isolation yields 9 basic blocks resembling Figure 3. θS1.α(i, x, y) = (0, i, 0, x, 0, y, α) θS2.α(i, x, y) = (0, i,

Molly reads the metadata generated by Clang to find which global variables are fields. In this example, there are two fields: f ront and back. The indexset domains are always zero-based and the dimension's sizes have been declared in the source files. The fields are accessed in the isolated statements with the elements accessed determined by the λ-function.

F = {f ront, back} I f ront,back = {(w, h) | 0 ≤ w, h < 1024} λS1.1(i, x, y) = {(f ront, (x -1, y))} λS1.2(i, x, y) = {(f ront, (x, y + 1))} λS1.3(i, x, y) = {(f ront, (x + 1, y))} λS1.4(i, x, y) = {(f ront, (x, y -1))} λS1.5(i, x, y) = {(f ront, (x, y))} λS1.6(i, x, y) = ∅ λS1.7(i, x, y) = {(back, (x, y))} λS2.1(i, x, y) = {(back, (x, y))} λS2.2(i, x, y) = {(f ront, (x, y))}
To compute the flow a data into and out of the SCoP, the virtual statements P rologue and Epilogue are added to the list of statements. They execute just once (Equation 5), upon entering the SCoP (Equation 6), respectively after the execution of all statements (Equation 7). The both touch all the data, the prologue is a virtual write access and the epilogue as a read access.

Ω ′ = Ω ∪ {P rologue, Epilogue} D P rologue,Epilogue = {()}
(5) θ P rologue () = (-1, 0, 0, 0, 0, 0, 0) (6) θ Epilogue () = (+1, 0, 0, 0, 0, 0, 0) (7) λ P rologue,Epilogue () = {(e, (w, h)) | e ∈ {f ront, back}, (w, g) ∈ Ie}

Now, Molly can compute the data flow between the statements, or def-use chains. For brevity, the exact set constraints are omitted here, but of course there can only be a flow between instances that are in the iteration domain D.

Equations 8 to 12 describe the direct data flow between the statements. These are data transfers within a SCoP whose chunk size will be computed later. Equations 13 to 17 are values that must be read from the local storage and potentially transferred to another node if the consumer statement is executed on another node. The data flows 18 to 19 are not overwritten by any other statement and therefore must be written back to the local storages.

Field flow deps:

(S1.7, (i, x, y)) δ (S2.1, (i, x, y)) (8) (S2.2, (i -1, x -1, y)) δ (S1.1, (i, x, y)) (9) (S2.2, (i -1, x, y + 1)) δ (S1.2, (i, x, y)) (10) (S2.2, (i -1, x + 1, y)) δ (S1.3, (i, x, y)) (11) (S2.2, (i -1, x, y -1)) δ (S1.4, (i, x, y)) (12)
(P rologue, ()) δ (S1.1, (0, x, y))

(13) (P rologue, ()) δ (S1.2, (0, x, y)) (14) (P rologue, ()) δ (S1.3, (0, x, y))

(15) (P rologue, ()) δ (S1.4, (0, x, y))

(16) (P rologue, ()) δ (S1.5, (0, x, y)) (17) (S1.7, (99, x, y)) δ (Epilogue, ())

(18) (S2.2, (99, x, y)) δ (Epilogue, ()) (19)
Equations 20 to 25 are non-field flows, i.e. they are stored and loaded from memory that is always local. Therefore, they force related instances to be executed on the same node.

(S1.1, (i, x, y)) δ (S1. 2,(i,x,y)) (20) (S1.2, (i, x, y)) δ (S1. 3,(i,x,y)) (21) (S1. 3,(i,x,y)) δ (S1.4, (i, x, y)) (22) (S1.4, (i, x, y)) δ (S1.6, (i, x, y)) (23) (S1.5, (i, x, y)) δ (S1.6, (i, x, y))

(24) (S2.1, (i, x, y)) δ (S2.2, (i, x, y))

The geometry of the nodes must be fixed here. In this example, the cluster consists of 64 nodes in an 8 × 8 mesh. The block distribution result in 128 × 128-sized tiles.

P = {(a, b) | 0 ≤ a, b < 8} π f ront,back (w, h) = (⌊w/128⌋, ⌊h/128⌋) (26)
Based on the owner-computes policy, one possible distribution of statement instances to executions is presented below. πS1.7(i, x, y) = π back (x, y) πS2.2(i, x, y) = π f ront (x, y)

The non-field dependencies 20 to 25 now force the rest of the statements to be computed on the same nodes as the consumers: πS1.1...S1.6(i, x, y) = πS1.7(i, x, y) = π back (x, y) πS1.2(i, x, y) = πS2.2(i, x, y) = π f ront (x, y)

Next, one has to find a chunking function for the intra-SCoP flow dependencies. This example shows the chunking of one of the flows, the one of Equation 9. The Algorithm 1 computes a level of independence of 2 since the for-loop on i serializes the two interior loops. In addition, no dependence is violated. The chunking function therefore is ϕS1.1(i, x, y) = (S1.1, (i, 1, 1)), or any other representative consumer with fixed i.

The transfer relation therefore becomes the Equation 27. Since the generator instances S2.2 are each executed on exactly one node, no unique producer has to be chosen and T ransf ers ′ = T ransf ers. The chunks become Equation 28.

For transferring the data between statements S2.2 and S1.1 two possibilities remain: The case π f ront (x -1, y) = π f ront (x, y), i.e. the value is used on the same node where it is computed, and π f ront (x -1, y) = π f ront (x, y) i.e. the value is computed on on the neighbor node. These are the elements {(f ront, (w, h)) ∈ I f ront | w = 128a -1, 128b ≤ h < 128(b + 1)} that need to be transferred from node (a -1, b) (where the producer executes) to (a, b) (where the consumer executes). Per source node/target node either 128 or 127 (at the extremities) values have to be transferred. A possible indexing function into the communication buffer can be f (w, h) = h -128b, where b = ⌊h/128⌋ as of Equation 26.

Therefore, on the producer node, Molly generates three statements to prepare the send buffer, write the value into the send buffer and then start the transfer. Analogously, there are three statements on the receiver side. For instance, Molly generates a statement Srecv on the receiving node (a, b) that calls molly recv((a -1, b)) with the properties

Drecv = {(i) | 1 ≤ i < 100} θrecv(i) = (0, i, 128a, 0, 128b, -0.5) πrecv(i) = {(a, b) ∈ P}

Experiments

The example program in Figure 1 has been executed on the Blue Gene/Q supercomputer in Jlich. The results are shown in Table 1, with varying field sizes (1024 in the Figure ) and on different number of nodes. The times are averaged over three executions of the SCoP with one warmup execution before.

The experiments show the approximately perfect weak scaling of the program and on this hardware. If the tile size per node is 128 2 , the 100 iterations take almost exactly one second, whatever the total number of nodes in the system. The nodes on the mesh surface have less work to do, therefore the execution with just 4 nodes takes just 0.9 seconds.

The Blue Gene/Q compute nodes have 16 physical cores with 4way SMT 2 each. Up to 64 ranks can be put onto one node, in which 3 Ranks per node (virtual MPI nodes per physical processor) 4 Tile size per rank 5 Million lattice-site updates per second (executions of hasLife) Of course, this serves more as a demonstration of the capabilities of Molly than to run the Game of Life. Any optimized Game of Life program reaches this performance without thousands of nodes with techniques that are not in Molly's scope. In addition, the code generated by Molly is not as optimal as it could. For instance, it generates a lot of unnecessary conversions between 32-and 64-bit integers.

Extensions

The current implementation is fairly limiting and yet unusable for practical uses. In the current state Molly is no more than a prototype. The most severe missing capabilities are

• Reductions are not supported. In the current implementation a reduction is detected as a non-field dependency that cases all statements to be executed on every node.

• Field elements must be scalars. As pointer arithmetic on element pointers are forbidden, no byte offset can be added before reading or writing to any member except the first.

• Support for periodic boundary conditions (arrays that are tori instead of meshes)

• Runtime information on the chosen data layout. This is necessary, for instance, to implement filling a field with data from a file, where each node loads a section of the file and stores it locally.

• Fields must be global variables. Fields local to a function or as a struct member do not work. Also, fields cannot be passed as arguments to functions.

• Conditionals inside the SCoP, even if statically determined, are unsupported. Conditionals on affine constraints can just be subtracted from the iteration domain. Other conditions might be handled as a no-op write of the condition is not met (F

[ k] = cond ? newValue : F[ k]).
At least the first four capabilities are required for LatticeQCD simulation and therefore planned as part of my PhD thesis. Some more visionary features are presented in the following.

T ransf ers = (S1.1, (i, 1, 1), S2.2, (i -1, x -1, y), (⌊(x -1)/128⌋, ⌊y/128⌋), S1.1, (i, x, y), (⌊x/128⌋, ⌊y/128⌋))

(27) Chunks = {(S2.2, (i -1, x -1, y), π f ront (x -1, y) producer , S1.1, (i, x, y), π f ront (x, y)) consumer , f ront, (x, y) field element } 0≤i<99
(28)

Explicit Data Distribution

The current implementation always uses a fixed block-distribution with equally sized blocks of data for every node. Data regions do not overlap. Further development will allows the programmer to define an arbitrary mapping by annotation the field declaration with a #pragma. For instance, the following annotation interchanges the index dimensions such that the data is stored in column-major format, instead of the default of row-major in C, by reversing the indexing. The home location(s) can be selected as well.

#pragma molly transform("{ [i,j] -> [j,i] }", 0) molly::array<bool,128> field;

Further ideas include tiling of memory areas, grouping elements for access using SIMD-style loads, using space-filling curves (Hilbert Curve, Z-Curve) for local indexing and compilergenerated indirect addressing (F[Idx[ k]]) if the index function becomes too complicated. One might even switch the data layout dynamically at runtime based on what the next operation on it will be 3 .

Furthermore, Molly might decide itself where to place data based on heuristics on the SCoPs that access them.

Explicit Statement Distribution

Instead of employing a policy like "owner computes", the user might explicitly state where to execute a statement. This gives more control to the programmer who might know better what gives the best performance or wants to try out different possibilities.

Conclusion

Molly is a compiler extension intended to take away the complexity of writing scientific programs for distributed memory machines by generating the necessary communication code. Such communication code often has the same patterns but is difficult to write because a lot of index computations, allocations, node synchronizations etc. are involved. The goal is not to take away control from the programmer who can still decide where to place and execute computations and data. Quite the contrary, the developer can try out different configurations to see which runs faster on a particular machine.

In contrast to previous attempts on automatic parallelization on distributed memory machines Molly does not try to be the one tool that fits it all but optimizes only special cases well. Nor does it introduce a new programming language that forces the developer to only use high-level constructs. The programmer does not need to write explicitly for a specific architecture but still has full control if needed for optimization. The integration into the general purpose compiler Clang+LLVM hopefully simplifies its use. It also has the advantage of taking away the necessity of reinventing the wheel for every language feature or optimization.

The first experiments give a promising outlook on what is possible for a large class of scientific programs. These include stencilcomputations, image processing, linear algebra, differential equations, quantum chemistry, quantum physics and generally anything 3 The author did this for the manual optimization of lattice QCD that can be expressed as a static control part (SCoP). It is not intended for other things like bioinformatics, data(base) lookups, or graphs computations. This work will be part of the author's PhD thesis that aims to automatically parallelize Lattice QCD programs for various types of cluster architectures. It grew from the author's experience of optimizing such a program manually with many different optimization where most involved writing just a lot of boilerplate code. A comparison of the manually and Molly-optimized program will appear in the PhD thesis.

Figure 1 .

 1 Figure 1. Conway's Game of Life

Figure 2 .

 2 Figure 2. Equivalent LLVM IR code of "back[x][y] = living"

1 :

 1 Input: parametric flow dependence (G, iG) δ (C, iC ) Result: chunking function ϕ for l=0 (entire SCoP) to NT -1 (innermost loop of C) do if not θ(G, i) 0,...,l <lex θ(C, i) 0,...,l then continue;; // (1)ϕ = lambda: (C, i) → (C,(i0, i1, . . . , i l , 0, . . . , 0)) ; δ ′ := apply rho on any instance (C, iC ) on the left or right side of δ (collapse) ; if ϕ(C, iC ) δ ′ * ϕ(C, iC ) then continue;; // (2) return ϕ ; end return lambda: (C, i) → (C, i) ; Algorithm Chunking function heuristic

Figure 3 .

 3 Figure 3. SCoP with isolated field accesses

#

  pragma molly where("{[i,x,y]->[floor(x/128),floor(y/128)]}") auto living = hasLife(front[x][y], neighbors);

Table 1 .

 1 Runtimes of the reduced Game of Life (Figure1) on Blue Gene/Q. case the execution time increases to just 1.7 seconds, but gets even a little faster with one rank per core due to in-memory transfers on the node.The program also scales as expected if the size of the tile per rank increases. Per 4-times larger time, the program takes 4 times longer to execute. The number of stencil executions per seconds remains approximately the same.

	2 Hardware threads

Unsurprisingly, there is no assembly instruction that can localize remote memory on DMMs. That's why they are termed that way.

Acknowledgments

Not yet written


Michael Kruse

Perfrewrite Program Complexity Analysis via Source Code Instrumentation

Keywords: profiling, computational complexity, source code instrumentation, LLVM Clang

Most program profiling methods output the execution time of one specific program execution, but not its computational complexity class in terms of the big-O notation. Perfrewrite is a tool based on LLVM's Clang compiler to rewrite a program such that it tracks semantic information while the program executes and uses it to guess memory usage, communication and computational complexity. While source code instrumentation is a standard technique for profiling, using it for deriving formulas is an uncommon approach.

Motivation

Several methods are common in order to deduce the performance characteristics of a program, usually used to identify a program's bottleneck that is most worthy to optimize. We are specially interested in HPC applications used by physicists. Our case study are the programs tmLQCD [START_REF] Jansen | tmLQCD: A Program Suite to Simulate Wilson Twisted Mass Lattice QCD[END_REF] and DD-HMC [START_REF] Lüscher | Schwarz-Preconditioned HMC Algorithm for Two-Flavour Lattice QCD[END_REF], both implementations of the lattice quantum chromodynamics (LQCD) simulation. The goal is to automatically derive the execution time in terms of the size of the input field.

The usual approach for program complexity analysis is done statically on the source code with tools such as PIPS [KAC + 96]. Unfortunately, challenges like pointer arithmetic make it near impossible for such tools to process the given programs. Another idea is to use modelling program such as PAMELA [vG93], but it requires to rewrite a program in its domain specific language and still does not support language constructs such as pointers.

Profiling with tools like gprof, oprofile, Tau and many more is the most mature technique, but do not return the complexity in big-O notation. This drawback can be coped with by running the program multiple times with different inputs and let a statistics tool (like GNU R) analyse it. To get a meaningful result also large input sizes have to be sampled. In case of shared cluster systems probably someone will complain about eating computation time just to find out how long a program takes to execute.

Approach

The main idea is to replace all relevant types in a C program by custom C++ classes. For instance, the double type gets replaced with a custom Double 2 class. Then, we rely on C++ operator overloading to call our method whenever an arithmetic operation is to be executed. In addition to just simply execute the operation, it increases a global FLOP counter. Similarly, malloc and free can be substituted to track memory usage as well as calls to MPI to track communication between cluster nodes.

By itself this does nothing more than standard profiling does. In addition, also loops are annotated using preprocessor macros. If the number of loop iterations is dependent on the program's input size, say n, it will multiply the counters in the loop body by n and divide by the number of loop iterations it actually executes in its configuration. In order to make this work, integral types are also replaced by custom C++ classes such that they carry annotations of their values' origin. The associated semantic information are the actual integer value, the value expressed as a term (a symbolic formula of the input size), and a typical integer value of a typical large input size.

For example, a program has two inputs N and M of size n and m. We execute the program using a configuration n = 8 and m = 4, a very small problem size. The symbolic representations are n term = "n" and m term = "m" respectively. When the values are used in an arithmetic operation, for instance they are multiplied, the annotated value becomes nm term = "n • m". The annotated large input size might be n large = 256 and m large = 128. The annotation after the multiplication is nm large = 32768. It is used whenever two values are compared since the formula representation has no total order, assuming that the large problem case is more relevant. One use is the update of the peak memory usage during the program's execution, avoiding a gigantic step function. In short, we execute the program for a small problem size while interpolating its characteristics for a large problem size.

As most profiling tools do, we log the begin and return of every function. This will allow to create a call tree with stats for every function. Of course this approach is not perfect. It works only with polynomial complexities. Loops must have the structure of for-loops with iteration count directly depending on the input. Data-dependent while-loops for instance are not possible, but can be worked-around by making the number of while-iterations an input size. But the advantage is that the techniques does not need to cope with logic unrelated to the complexity. No need to handle pointer arithmetic, polymorphism, the halting problem, etc. because the program is executed and these intractabilities are evaluated naturally.

The implemented approach executes every loop body at most twice. The first iteration may initialize global fields3 , the second iteration is assumed to represent all following iterations. The is sufficient for the two case study programs mentioned in the previous section. The program's result will be wrong, but we are not interested in it anyway.

Implementation

A simple textual insertion and replacement is not sufficient. We want to track memory blocks, therefore we also need to replace pointers to them. Pointers are typed, meaning the pointer replacement class needs to be a template whose variable type declaration syntax is different than a pointer's. To instrument manually is a tedious task to be avoided if possible. In addition, instrumenting only necessary parts of the code will safe us from incompatibilities in unaffected and working program sections. E. g. C++ classes cannot perfectly emulate the behaviour of C pointers. If necessary the user might need to adapt the code manually, like when a loop syntax does not exactly match any implemented patters. This is why the instrumented code must stay readable.

The Clang compiler [Cla] has the necessary facilities for a semantic analysis of C and C++ code. It also retains the source code locations in the internal representation so our tool can locate and replace the parts that have to be replaced. Clang already includes an Objective-C to C++ translator that works similarly. Also, a full C++ compiler already includes facilities to get the return type of an expression involving custom types. The principal tool working is shown in Figure 2. Intrinsic types are transformed into custom C++ types whenever it is assigned an expression that evaluates to a custom type or the types are incompatible in some context (e.g. when passing by reference). Types may influence themselves, so this is done until a fixpoint is reached. When types types declared extern are changed -this includes non-static functions -this change must be propagated to the other source files of the program. If the change is part of a macro, the macro has to be expanded before the instrumentation. Clang does not retain the complete information of how a macro was expanded so the complete parse must be repeated after macro expansion.

The instrumented code can now be compiled using any C++ compiler. The code changes require a custom runtime to be linked against it containing the implementation of the replacement classes. The implementation of the symbolic representation is from GiNaC [BFK + ], a computer algebra system using C++ as control language. When the compilee runs, it writes a call tree into a .dot file that can be processed by Graphviz [AT ]. Figure 3 shows such a call tree. For every function call it shows how often it has been called, the number of floating point operations, peak memory usage and MPI calls, all as complexities depending on the input parameters.

Figure 1 :

 1 Figure 1: An example program before (left) and after instrumentation (right)

Figure 2 :

 2 Figure 2: Processing pipeline.

Figure 3 :

 3 Figure 3: Graphviz rendering of an example program.

E-mail: michael.kruse@inria.fr

of the kind if (!g_field) g_field = malloc(...)