You are an expert at summarizing long articles. Proceed to summarize the following text: 
the need for rapid implementation of high performance , robust , and portable finite element methods has led to approaches based on automated code generation . this has been proven successful in the context of the fenics @xcite and firedrake @xcite projects . in these frameworks , the weak variational form of a problem is expressed in a high level mathematical syntax by means of the domain - specific language ufl @xcite . this mathematical specification is used by a domain - specific compiler , known as a form compiler , to generate low - level c or c++ code for the integration over a single element of the computational mesh of the variational problem s left and right hand side operators . the code for assembly operators must be carefully optimized : as the complexity of a variational form increases , in terms of number of derivatives , pre - multiplying functions , or polynomial order of the chosen function spaces , the operation count increases , with the result that assembly often accounts for a significant fraction of the overall runtime .    as demonstrated by the substantial body of research on the topic , automating the generation of such high performance implementations poses several challenges . this is a result of the complexity inherent in the mathematical expressions involved in the numerical integration , which varies from problem to problem , and the particular structure of the loop nests enclosing the integrals . general - purpose compilers , such as those by _ gnu _ and _ intel _ , fail to exploit the structure inherent in the expressions , thus producing sub - optimal code ( i.e. , code which performs more floating - point operations , or `` flops '' , than necessary ; we show this in section  [ sec : perf - results ] ) . research compilers , for instance those based on polyhedral analysis of loop nests , such as pluto @xcite , focus on parallelization and optimization for cache locality , treating issues orthogonal to the question of minimising flops . the lack of suitable third - party tools has led to the development of a number of domain - specific code transformation ( or synthesizer ) systems . show how automated code generation can be leveraged to introduce optimizations that a user should not be expected to write `` by hand '' . and employ mathematical reformulations of finite element integration with the aim of minimizing the operation count . in  , the effects and the interplay of generalized code motion and a set of low level optimizations are analysed . it is also worth mentioning two new new form compilers , uflacs @xcite and tsfc @xcite , which particularly target the compilation time challenges of the more complex variational forms . the performance evaluation in section  [ sec : perf - results ] includes most of these systems . however , in spite of such a considerable research effort , there is still no answer to one fundamental question : can we automatically generate an implementation of a form which is optimal in the number of flops executed ? in this paper , we formulate an approach that solves this problem for a particular class of forms and provides very good approximations in all other cases . in particular , we will define `` local optimality '' , which relates operation count with inner loops . in summary , our contributions are as follows :    * we formalize the class of finite element integration loop nests and we build the space of legal transformations impacting their operation count . * we provide an algorithm to select points in the transformation space . the algorithm uses a cost model to : ( i ) understand whether a transformation reduces or increases the operation count ; ( ii ) choose between different ( non - composable ) transformations . * we demonstrate that our approach systematically leads to a local optimum . we also explain under what conditions of the input problem global optimality is achieved . * we integrate our approach with a compiler , coffee , which is in use in the firedrake framework . * we experimentally evaluate using a broader suite of forms , discretizations , and code generation systems than has been used in prior research . this is essential to demonstrate that our optimality model holds in practice . in addition , in order to place coffee on the same level as other code generation systems from the viewpoint of low level optimization ( which is essential for a fair performance comparison ) :    * we introduce a transformation based on symbolic execution that allows irrelevant floating point operations to be skipped ( for example those involving zero - valued quantities ) .    after reviewing basic concepts in finite element integration , in section  [ sec : lnopt ] we introduce a set of definitions mapping mathematical properties to the level of loop nests . this step is an essential precursor to the definition of the two algorithms  sharing elimination ( section  [ sec : sharing - elimination ] ) and pre - evaluation ( section  [ sec : pre - evaluation ] )  through which we construct the space of legal transformations . the main transformation algorithm in section  [ sec : optimal - synthesis ] delivers the local optimality claim by using a cost model to coordinate the application of sharing elimination and pre - evaluation . we elaborate on the correctness of the methodology in section  [ sec : proof ] . the numerical experiments are showed in section  [ sec : perf - results ] . we conclude discussing the limitations of the algorithms presented and future work . we review finite element integration using the same notation and examples adopted in   and  . consider the weak formulation of a linear variational problem : @xmath1 where @xmath2 and @xmath3 are , respectively , a bilinear and a linear form . the set of _ trial _ functions @xmath4 and the set of _ test _ functions @xmath5 are suitable discrete function spaces . for simplicity , we assume @xmath6 . let @xmath7 be the set of basis functions spanning @xmath4 . the unknown solution @xmath8 can be approximated as a linear combination of the basis functions @xmath7 . from the solution of the following linear system it is possible to determine a set of coefficients to express @xmath8 : @xmath9 in which @xmath10 and @xmath11 discretize @xmath2 and @xmath3 respectively : @xmath12 the matrix @xmath10 and the vector @xmath11 are assembled and subsequently used to solve the linear system through ( typically ) an iterative method . we focus on the assembly phase , which is often characterized as a two - step procedure : _ local _ and _ global _ assembly . local assembly is the subject of this article . it consists of computing the contributions of a single element in the discretized domain to the equation s approximated solution . during global assembly , these local contributions are coupled by suitably inserting them into @xmath10 and @xmath11 . we illustrate local assembly in a concrete example , the evaluation of the local element matrix for a laplacian operator . consider the weighted poisson equation : @xmath13 in which @xmath8 is unknown , while @xmath14 is prescribed . the bilinear form associated with the weak variational form of the equation is : @xmath15 the domain @xmath16 of the equation is partitioned into a set of cells ( elements ) @xmath17 such that @xmath18 and @xmath19 . by defining @xmath20 as the set of basis functions with support on the element @xmath21 ( i.e. those which do not vanish on this element ) , we can express the local element matrix as @xmath22 the local element vector @xmath3 can be determined in an analogous way . it has been shown ( for example in  ) that local element tensors can be expressed as a sum of integrals over @xmath21 , each integral being the product of derivatives of functions from sets of discrete spaces and , possibly , functions of some spatially varying coefficients . an integral of this form is called _ monomial_.      quadrature schemes are typically used to numerically evaluate @xmath23 . for convenience , a reference element @xmath24 and an affine mapping @xmath25 to any element @xmath26 are introduced . this implies that a change of variables from reference coordinates @xmath27 to real coordinates @xmath28 is necessary any time a new element is evaluated . the basis functions @xmath29 are then replaced with local basis functions @xmath30 such that @xmath31 . the numerical integration of over an element @xmath21 can then be expressed as follows : @xmath32 where @xmath33 is the number of integration points , @xmath34 the quadrature weight at the integration point @xmath35 , @xmath36 the dimension of @xmath16 , @xmath37 the number of degrees of freedom associated to the local basis functions , and @xmath38 the determinant of the jacobian of the aforementioned change of coordinates .      by exploiting the linearity , associativity and distributivity of the relevant mathematical operators , we can rewrite as @xmath39 a generalization of this transformation was introduced in  @xcite . since it only involves reference element terms , the quadrature sum can be pre - evaluated and reused for each element . the evaluation of the local tensor can then be abstracted as @xmath40 in which the pre - evaluated _ reference tensor _ , @xmath41 , and the cell - dependent _ geometry tensor _ , @xmath42 , are exposed .      depending on form and discretization , the relative performance of the two modes , in terms of the operation count , can vary quite dramatically . the presence of derivatives or coefficient functions in the input form increases the rank of the geometry tensor , making the traditional quadrature mode preferable for sufficiently complex forms . on the other hand , speed - ups from adopting tensor mode can be significant in a wide class of forms in which the geometry tensor remains sufficiently small . the discretization , particularly the polynomial order of trial , test , and coefficient functions , also plays a key role in the resulting operation count . these two modes are implemented in the fenics form compiler @xcite . in this compiler , a heuristic is used to choose the most suitable mode for a given form . it consists of analysing each monomial in the form , counting the number of derivatives and coefficient functions , and checking if this number is greater than a constant found empirically @xcite . we will return to the efficacy of this approach in section [ sec : perf - results ] . one of the objectives of this paper is to produce a system that goes beyond the dichotomy between quadrature and tensor modes . we will reason in terms of loop nests , code motion , and code pre - evaluation , searching the entire implementation space for an optimal synthesis . in this section , we characterize global and local optimality for finite element integration as well as the space of legal transformations that needs be explored to achieve them . the method by which exploration is performed is discussed in section  [ sec : optimal - synthesis ] .      in order to make the article self - contained , we start with reviewing basic compiler terminology . a perfect loop nest is a loop whose body either 1 ) comprises only a sequence of non - loop statements or 2 ) is itself a perfect loop nest . if this condition does not hold , a loop nest is said to be imperfect . an independent basic block is a sequence of statements such that no data dependencies exist between statements in the block . we focus on perfect nests whose innermost loop body is an independent basic block . a straightforward property of this class is that hoisting invariant expressions from the innermost to any of the outer loops or the preheader ( i.e. , the block that precedes the entry point of the nest ) is always safe , as long as any dependencies on loop indices are honored . we will make use of this property . the results of this section could also be generalized to larger classes of loop nests , in which basic block independence does not hold , although this would require refinements beyond the scope of this paper .    by mapping mathematical properties to the loop nest level , we introduce the concepts of a _ linear loop _ and , more generally , a ( perfect ) multilinear loop nest . [ def : linear - loop ] a loop @xmath3 defining the iteration space @xmath43 through the iteration variable @xmath44 , or simply @xmath45 , is linear if in its body    1 . @xmath44 appears only as an array index , and 2 . whenever an array @xmath2 is indexed by @xmath44 ( @xmath46 $ ] ) , all expressions in which this appears are affine in @xmath46 $ ] . [ def : multi - linear - loop ] a multilinear loop nest of arity @xmath37 is a perfect nest composed of @xmath37 loops , in which all of the expressions appearing in the body of the innermost loop are affine in each loop @xmath45 separately . we will show that multilinear loop nests , which arise naturally when translating bilinear or linear forms into code , are important because they have a structure that we can take advantage of to reach a local optimum . we define two other classes of loops . [ def : i - loop ] a loop @xmath45 is said to be a reduction loop if in its body    1 . @xmath44 appears only as an array index , and 2 .   for each augmented assignment statement @xmath47 ( e.g. , an increment ) , arrays indexed by @xmath44 appear only on the right hand side of @xmath47 . [ def : e - loop ] a loop @xmath45 is said to be an order - free loop if its iterations can be executed in any arbitrary order . consider equation  [ eq : quadrature ] and the ( abstract ) loop nest implementing it illustrated in figure  [ code : loopnest ] . the imperfect nest @xmath48 $ ] comprises an order - free loop @xmath49 ( over elements in the mesh ) , a reduction loop @xmath45 ( performing numerical integration ) , and a multilinear loop nest @xmath50 $ ] ( over test and trial functions ) . in the body of @xmath51 , one or more statements evaluate the local tensor for the element @xmath52 . expressions ( the right hand side of a statement ) result from the translation of a form in high level matrix notation into code . in particular , @xmath53 is the number of monomials ( a form is a sum of monomials ) , @xmath54 ( @xmath55 ) represents the product of a coefficient function ( e.g. , the inverse jacobian matrix for the change of coordinates ) with test or trial functions , and @xmath56 is a function of coefficients and geometry . we do not pose any restrictions on function spaces ( e.g. , scalar- or vector - valued ) , coefficient expressions ( linear or non - linear ) , differential and vector operators , so @xmath56 can be arbitrarily complex . we say that such an expression is in _ normal form _ , because the algebraic structure of a variational form is intact : products have not yet been expanded , distinct monomials can still be identified , and so on . this brings us to formalize the class of loop nests that we aim to optimize . [ def : fem - loopnest ] a finite element integration loop nest is a loop nest in which the following appear , in order : an imperfect order - free loop , an imperfect ( perfect only in some special cases ) , linear or non - linear reduction loop , and a multilinear loop nest whose body is an independent basic block in which expressions are in normal form . we then characterize optimality for a finite element integration loop nest as follows . [ def : mln - optimality ] let @xmath57 be a generic loop nest , and let @xmath58 be a transformation function @xmath59 such that @xmath60 is semantically equivalent to @xmath57 ( possibly , @xmath61 ) . we say that @xmath62 is an optimal synthesis of @xmath57 if the total number of operations ( additions , products ) performed to evaluate the result is minimal .    the concept of local optimality , which relies on the particular class of _ flop - decreasing _ transformations , is also introduced .    a transformation which reduces the operation count is called flop - decreasing . [ def : mln - quasi - optimality ] given @xmath57 , @xmath60 and @xmath58 as in definition  [ def : mln - optimality ] , we say that @xmath62 is a locally optimal synthesis of @xmath57 if :    * the number of operations ( additions , products ) in the innermost loops performed to evaluate the result is minimal , and * @xmath58 is expressed as composition of flop - decreasing transformations . the restriction to flop - decreasing transformations aims to exclude those apparent optimizations that , to achieve flop - optimal innermost loops , would rearrange the computation at the level of the outer loops causing , in fact , a global increase in operation count . we also observe that definitions  [ def : mln - optimality ] and  [ def : mln - quasi - optimality ] do not take into account memory requirements . if the execution of loop nest were memory - bound  the ratio of operations to bytes transferred from memory to the cpu being too low  then optimizing the number of flops would be fruitless . henceforth we assume we operate in a cpu - bound regime , evaluating arithmetic - intensive expressions . in the context of finite elements , this is often true for more complex multilinear forms and/or higher order elements . achieving optimality in polynomial time is not generally feasible , since the @xmath56 sub - expressions can be arbitrarily unstructured . however , multilinearity results in a certain degree of regularity in @xmath54 and @xmath55 . in the following sections , we will elaborate on these observations and formulate an approach that achieves : ( i ) at least a local optimum in all cases ; ( ii ) global optimality whenever the monomials are `` sufficiently structured '' . to this purpose , we will construct :    * the space of legal transformations impacting the operation count ( sections  [ sec : sharing - elimination ]  [ sec : mem - const ] ) * an algorithm to select points in the transformation space ( section  [ sec : optimal - synthesis ] )      we start with introducing the fundamental notion of sharing . a statement within a loop nest @xmath57 presents sharing if at least one of the following conditions hold :    spatial sharing : :    there are at least two symbolically identical sub - expressions temporal sharing : :    there is at least one non - trivial sub - expression ( e.g. , an addition or    a product ) that is redundantly executed because it is independent of    @xmath63 .    to illustrate the definition , we show in figure  [ code : multi_loopnest ] how sharing evolves as factorization and code motion are applied to a trivial multilinear loop nest . in the original loop nest ( figure  [ code : multi_loopnest_a ] ) , spatial sharing is induced by the symbol @xmath64 . factorization eliminates spatial sharing and creates temporal sharing ( figure  [ code : multi_loopnest_b ] ) . finally , generalized code motion @xcite , which hoists sub - expressions that are redundantly executed by at least one loop in the nest , leads to optimality ( figure  [ code : multi_loopnest_c ] ) . = 13pt   = 13pt   = 4pt    in this section , we study _ sharing elimination _ , a transformation that aims to reduce the operation count by removing sharing through the application of expansion , factorization , and generalized code motion . if the objective were reaching optimality and the expressions lacked structure , a transformation of this sort would require solving a large combinatorial problem  for instance to evaluate the impact of all possible factorizations . our sharing elimination strategy , instead , exploits the structure inherent in finite element integration expressions to guarantee , after coordination with other transformations ( an aspect which we discuss in the following sections ) , local optimality . global optimality is achieved if stronger preconditions hold . setting local optimality , rather than optimality , as primary goal is essential to produce simple and computationally efficient algorithms  two necessary conditions for integration with a compiler . finite element expressions can be seen as composition of operations between tensors . often , the optimal implementation strategy for these operations is to be determined out of two alternatives . for instance , consider @xmath65 , with @xmath66 being the transposed inverse jacobian matrix for the change of ( two - dimensional ) coordinates , and @xmath67 a generic two - dimensional vector . the tensor operation will reduce to the scalar expression @xmath68 , in which @xmath69 and @xmath70 represent components of @xmath67 that depend on @xmath45 . to minimize the operation count for expressions of this kind , we have two options :    [ strategy : i ] eliminating temporal sharing through generalized code motion .    [ strategy : ii ] eliminating spatial sharing first  through product expansion and factorization  and temporal sharing afterwards , again through generalized code motion .    in the current example , we observe that , depending on the size of @xmath45 , applying strategy  [ strategy : ii ] could reduce the operation count since the expression would be recast as @xmath71 and some hoistable sub - expressions would be exposed . on the other hand , strategy  [ strategy : i ] would have no effect as @xmath67 only depends on a single loop , @xmath45 . in general , the choice between the two strategies depends on multiple factors : the loop sizes , the increase in operation count due to expansion ( in strategy  [ strategy : ii ] ) , and the gain due to code motion . a second application of strategy  [ strategy : ii ] was provided in figure  [ code : multi_loopnest ] . these examples motivate the introduction of a particular class of expressions , for which the two strategies assume notable importance . [ def : struct - expr ] we say that an expression is `` structured along a loop nest @xmath57 '' if and only if , for every symbol @xmath72 depending on at least one loop in @xmath57 , the spatial sharing of @xmath72 may be eliminated by factorizing all occurrences of @xmath72 in the expression .    [ prop : multi - struct ] an expression along a multilinear loop nest is structured . this follows directly from definition  [ def : linear - loop ] and definition  [ def : multi - linear - loop ] , which essentially restrict the number of occurrences of a symbol @xmath72 in a summand to at most 1 . if @xmath57 were an arbitrary loop nest , a given symbol @xmath72 could appear everywhere ( e.g. , @xmath37 times in a summand and @xmath53 times in another summand with @xmath73 , as argument of a higher level function , in the denominator of a division ) , thus posing the challenge of finding the factorization that maximizes temporal sharing . if @xmath57 is instead a finite element integration loop nest , thanks to proposition  [ prop : multi - struct ] the space of flop - decreasing transformations is constructed by `` composition '' of strategy  [ strategy : i ] and strategy  [ strategy : ii ] , as illustrated in algorithm  [ algo : sharing - elimination ] . finally , we observe that the @xmath56 sub - expressions can sometimes be considered `` weakly structured '' . this happens when a relaxed version of definition  [ def : struct - expr ] applies , in which the factorization of @xmath72 only `` minimizes '' ( rather than `` eliminates '' ) spatial sharing ( for instance , in the complex hyperelastic model analyzed in section  [ sec : perf - results ] ) . weak structure will be exploited by algorithm  [ algo : sharing - elimination ] in the attempt to achieve optimality . algorithm  [ algo : sharing - elimination ] describes sharing elimination assuming as input a tree representation of the loop nest . it makes use of the following notation and terminology :    * _ multilinear operand _ : any @xmath54 or @xmath55 in the input expression . * _ multilinear symbol _ : a symbol appearing within a multilinear operand depending on @xmath74 or @xmath51 ( e.g. , test functions , first order derivatives of test functions , etc . ) . examples will be provided in section  [ sec : se - examples ] .    ' '' ''    [ algo : sharing - elimination ] the input of the algorithm is a tree representation a finite element integration loop nest .    1 .   perform a depth - first visit of the loop tree to collect and partition multilinear operands into disjoint sets , @xmath75 . @xmath76 is such that all multilinear operands in each @xmath77 share the same set of multilinear symbols @xmath78 , whereas there is no sharing across different partitions . for all multilinear operands in @xmath77 such that @xmath79 , apply strategy  [ strategy : i ] . + _ note : as a consequence of proposition  [ prop : multi - struct ] , @xmath80 and @xmath81 represent the number of products in the innermost loop induced by @xmath82 if strategy  [ strategy : i ] or strategy  [ strategy : ii ] were applied _ 2 .   for each sub - expression @xmath83 depending on exactly one linear loop , collect the multilinear symbols and the temporaries produced at step ( 1 ) . partition them into disjoint sets , @xmath84 , such that @xmath85 includes all instances of a given symbol in @xmath83 . apply strategy  [ strategy : ii ] factorizing the symbols in each @xmath85 , provided that this leads to a reduction in operation count ; otherwise , apply strategy  [ strategy : i ] + _ note : the last check ensures the flop - decreasing nature of the transformation . in the cases in which expansion outweighs code motion , strategy  [ strategy : i ] is preferred . _ + _ note : the expansion cost is a function of the products wrapping a symbol ( how many of them and their arity ) , so it can be determined through tree visits . build the _ sharing graph _ @xmath86 . each @xmath87 represents a multilinear symbol or a temporary produced by the previous steps . an edge @xmath88 , @xmath89 indicates that a product @xmath90 would appear if the sub - expressions including @xmath91 and @xmath92 were expanded . + _ note : the following steps will only impact bilinear forms , since otherwise @xmath93 . _ 4 . partition @xmath47 into disjoint sets , @xmath94 , such that @xmath95 includes all instances of a given symbol @xmath96 in the expression . transform @xmath97 by merging @xmath98 into a unique vertex @xmath96 ( taking the union of the edges ) , provided that factorizing @xmath99 $ ] would not cause an increase in operation count . map @xmath97 to an integer linear programming ( ilp ) model for determining how to optimally apply strategy  [ strategy : ii ] . the solution is the set of symbols that will be factorized by strategy  [ strategy : ii ] . let @xmath100 ; the ilp model then is as follows : @xmath101 6 . perform a depth - first visit of the loop tree and , for each yet unhandled or hoisted expression , apply the most profitable between strategy  [ strategy : i ] and strategy  [ strategy : ii ] . + _ note : this pass speculatively assumes that expressions are ( weakly ) structured along the reduction loop . if the assumption does not hold , the operation count will generally be sub - optimal because only a subset of factorizations and code motion opportunities may eventually be considered . _    ' '' ''    although the primary goal of algorithm  [ algo : sharing - elimination ] is operation count minimization within the multilinear loop nest , the enforcement of flop - decreasing transformations ( steps ( 2 ) and ( 4 ) ) and the re - scheduling of sub - expressions within outer loops ( last step ) also attempt to optimize the loop nest globally . we will further elaborate this aspect in section  [ sec : proof ] . consider again figure  [ code : multi_loopnest_a ] . we have @xmath102 , with @xmath103 , @xmath104 , and @xmath105 . for all @xmath106 , we have @xmath107 , although applying strategy  [ strategy : i ] in step ( 1 ) has no effect . the sharing graph is @xmath108 , and @xmath109 . the ilp formulation leads to the code in figure  [ code : multi_loopnest_c ] .    in figure [ code : poisson ] , algorithm  [ algo : sharing - elimination ] is executed in a very simple realistic scenario , which originates from the bilinear form of a poisson equation in two dimensions . we observe that @xmath110 , with @xmath111 and @xmath112 . in addition , @xmath113 , so strategy  [ strategy : i ] is applied to both partitions ( step ( 1 ) ) . we then have ( step ( 3 ) ) @xmath114 . since there are no more factorization opportunities , the ilp formulation becomes irrelevant . = 7pt   = 7pt    for reasons of space , further examples , including the hyperelastic model evaluated in section  [ sec : perf - results ] and other non - trivial ilp instances , are made available online . sharing elimination uses three operators : expansion , factorization , and code motion . in this section , we discuss the role and legality of a fourth operator : reduction pre - evaluation . we will see that what makes this operator special is the fact that there exists a single point in the transformation space of a monomial ( i.e. , a specific factorization of test , trial , and coefficient functions ) ensuring its correctness . we start with an example . consider again the loop nest and the expression in figure  [ code : loopnest ] . we pose the following question : are we able to identify sub - expressions for which the reduction induced by @xmath45 can be pre - evaluated , thus obtaining a decrease in operation count proportional to the size of @xmath45 , @xmath43 ? the transformation we look for is exemplified in figure  [ code : loopnest_rednored ] with a simple loop nest . the reader may verify that a similar transformation is applicable to the example in figure  [ code : poisson_a ] . = 19pt   = 5pt    pre - evaluation can be seen as the generalization of tensor contraction ( section  [ sec : tc ] ) to a wider class of sub - expressions . we know that multilinear forms can be seen as sums of monomials , each monomial being an integral over the equation domain of products ( of derivatives ) of functions from discrete spaces . a monomial can always be reduced to the product between a `` reference '' and a `` geometry '' tensor . in our model , a reference tensor is simply represented by one or more sub - expressions independent of @xmath49 , exposed after particular transformations of the expression tree . this leads to the following algorithm .    ' '' '' [ algo : pre - evaluation ] consider a finite element integration loop nest @xmath115 $ ] . we dissect the normal form input expression into distinct sub - expressions , each of them representing a monomial . each sub - expression is then factorized so as to split constants from @xmath116$]-dependent terms . this transformation is feasible , as a consequence of the results in  . these @xmath116$]-dependent terms are hoisted outside of @xmath57 and stored into temporaries . as part of this process , the reduction induced by @xmath45 is computed by means of symbolic execution . finally , @xmath45 is removed from @xmath57 .    ' '' ''    the pre - evaluation of a monomial introduces some critical issues :    1 .   depending on the complexity of a monomial , a certain number , @xmath117 , of temporary variables is required if pre - evaluation is performed . such temporary variables are actually @xmath37-dimensional arrays of size @xmath47 , with @xmath37 and @xmath47 being , respectively , the arity and the extent ( iteration space size ) of the multilinear loop nest ( e.g. , @xmath118 and @xmath119 in the case of bilinear forms ) . for certain values of @xmath120 , pre - evaluation may dramatically increase the working set , which may be counter - productive for actual execution time the transformations exposing @xmath116$]-dependent terms increase the arithmetic complexity of the expression ( e.g. , expansion tends to increase the operation count ) . this could outweigh the gain due to pre - evaluation . 3 .   a strategy for coordinating sharing elimination and pre - evaluation is needed . we observe that sharing elimination inhibits pre - evaluation , whereas pre - evaluation could expose further sharing elimination opportunities . we expand on point ( 1 ) in the next section , while we address points ( 2 ) and ( 3 ) in section  [ sec : optimal - synthesis ] . we have just observed that the code motion induced by monomial pre - evaluation may dramatically increase the working set size . even more aggressive code motion strategies are theoretically conceivable . imagine @xmath57 is enclosed in a time stepping loop . one could think of exposing ( through some transformations ) and hoisting time - invariant sub - expressions for minimizing redundant computation at each time step . the working set size would then increase by a factor @xmath121 , and since @xmath122 , the gain in operation count would probably be outweighed , from a runtime viewpoint , by a much larger memory pressure .    since , for certain forms and discretizations , hoisting may cause the working set to exceed the size of some level of local memory ( e.g. the last level of private cache on a conventional cpu , the shared memory on a gpu ) , we introduce the following _ memory constraints_.    [ const : le ] the size of a temporary due to code motion must not be proportional to the size of @xmath49 . [ const : th ] the total amount of memory occupied by the temporaries due to code motion must not exceed a certain threshold , ` t_h ` .    constraint  [ const : le ] is a policy decision that the compiler should not silently consume memory on global data objects . it has the effect of shrinking the transformation space . constraint  [ const : th ] has both theoretical and practical implications , which will be carefully analyzed in the next sections . in this section , we build a transformation algorithm that , given a memory bound , systematically reaches a local optimum for finite element integration loop nests . we address the two following issues :    1 . _ coordination of pre - evaluation and sharing elimination . _ recall from section  [ sec : pre - evaluation ] that pre - evaluation could either increase or decrease the operation count in comparison with that achieved by sharing elimination . 2 .   _ optimizing over composite operations . _ consider a form comprising two monomials @xmath123 and @xmath124 . assume that pre - evaluation is profitable for @xmath123 but not for @xmath124 , and that @xmath123 and @xmath124 share at least one term ( for example some basis functions ) . if pre - evaluation were applied to @xmath123 , sharing between @xmath123 and @xmath124 would be lost . we then need a mechanism to understand which transformation pre - evaluation or sharing elimination  results in the highest operation count reduction when considering the whole set of monomials ( i.e. , the expression as a whole ) . let @xmath125 be a cost function that , given a monomial @xmath126 , returns the gain / loss achieved by pre - evaluation over sharing elimination . in particular , we define @xmath127 , where @xmath128 and @xmath129 represent the operation counts resulting from applying sharing elimination and pre - evaluation , respectively . thus pre - evaluation is profitable for @xmath53 if and only if @xmath130 . we return to the issue of deriving @xmath128 and @xmath129 in section  [ sec : op_count ] . having defined @xmath131 , we can now describe the transformation algorithm ( algorithm  [ algo : gamma ] ) .    ' '' ''    [ algo : gamma ] the algorithm has three main phases : initialization ( step 1 ) ; determination of the monomials preserving the memory constraints that should be pre - evaluated ( steps 2 - 4 ) ; application of pre - evaluation and sharing elimination ( step 5 ) .    1 .   perform a depth - first visit of the expression tree and determine the set of monomials @xmath132 . let @xmath47 be the subset of monomials @xmath53 such that @xmath133 . the set of monomials that will _ potentially _ be pre - evaluated is @xmath134 . + _ note : there are two fundamental reasons for not pre - evaluating @xmath135 straight away : 1 ) the potential presence of spatial sharing between @xmath123 and @xmath136 , which impacts the search for the global optimum ; 2 ) the risk of breaking constraint  [ const : th ] . build the set @xmath137 of all possible bipartitions of @xmath82 . let @xmath138 be the dictionary that will store the operation counts of different alternatives . 3 .   discard @xmath139 if the memory required after applying pre - evaluation to the monomials in @xmath140 exceeds @xmath141 ( see constraint  [ const : th ] ) ; otherwise , add @xmath142 = \mathrm{\theta}^{se}(s \cup b_s ) + \mathrm{\theta}^{pre}(b_p)$ ] . + _ note : @xmath143 is in practice very small , since even complex forms usually have only a few monomials . this pass can then be accomplished rapidly as long as the cost of calculating @xmath128 and @xmath129 is negligible . we elaborate on this aspect in section  [ sec : op_count ] . _ take @xmath144 $ ] . 5 .   apply pre - evaluation to all monomials in @xmath140 . apply sharing elimination to all resulting expressions . + _ note : because of the reuse of basis functions , pre - evaluation may produce some identical tables , which will be mapped to the same temporary variable . sharing elimination is therefore transparently applied to all expressions , including those resulting from pre - evaluation . _    ' '' ''    the output of the transformation algorithm is provided in figure  [ code : loopnest - opt ] , assuming as input the loop nest in figure  [ code : loopnest ] . we tie up the remaining loose end : the construction of the cost function @xmath131 . we recall that @xmath145 , with @xmath128 and @xmath129 representing the operation counts after applying sharing elimination and pre - evaluation . since @xmath131 is deployed in a working compiler , simplicity and efficiency are essential characteristics . in the following , we explain how to derive these two values . the most trivial way of evaluating @xmath128 and @xmath129 would consist of applying the actual transformations and simply count the number of operations . this would be tolerable for @xmath128 , as algorithm  [ algo : sharing - elimination ] tends to have negligible cost . however , the overhead would be unacceptable if we applied pre - evaluation  in particular , symbolic execution  to all bipartitions analyzed by algorithm  [ algo : gamma ] . we therefore seek an analytic way of determining @xmath129 . the first step consists of estimating the _ increase factor _ , @xmath146 . this number captures the increase in arithmetic complexity due to the transformations exposing pre - evaluation opportunities . for context , consider the example in figure  [ code : increase_factor ] . one can think of this as the ( simplified ) loop nest originating from the integration of the action of a mass matrix . the sub - expression ` f_0*b_{i0}+f_1*b_{i1}+f_2*b_{i2 } ` represents the coefficient @xmath147 over ( tabulated ) basis functions ( array @xmath137 ) . in order to apply pre - evaluation , the expression needs be transformed to separate @xmath147 from all @xmath116$]-dependent quantities ( see algorithm  [ algo : pre - evaluation ] ) . by product expansion , we observe an increase in the number of @xmath50$]-dependent terms of a factor @xmath148 .    in general , however , determining @xmath146 is not so straightforward since redundant tabulations may result from common sub - expressions . consider the previous example . one may add one coefficient in the same function space as @xmath147 , repeat the expansion , and observe that multiple sub - expressions ( e.g. , @xmath149 and @xmath150 ) will reduce to identical tables . to evaluate @xmath146 , we then use combinatorics . we calculate the @xmath151-combinations with repetitions of @xmath37 elements , where : ( i ) @xmath151 is the number of ( derivatives of ) coefficients appearing in a product ; ( ii ) @xmath37 is the number of unique basis functions involved in the expansion . in the original example , we had @xmath152 ( for @xmath153 , @xmath154 , and @xmath155 ) and @xmath156 , which confirms @xmath157 . in the modified example , there are two coefficients , so @xmath158 , which means @xmath159 . if @xmath160 ( the extent of the reduction loop ) , we already know that pre - evaluation will not be profitable . intuitively , this means that we are introducing more operations than we are saving from pre - evaluating @xmath45 . if @xmath161 , we still need to find the number of terms @xmath162 such that @xmath163 . the mass matrix monomial in figure  [ code : increase_factor ] is characterized by the dot product of test and trial functions , so trivially @xmath164 . in the example in figure  [ code : poisson ] , instead , we have @xmath165 after a suitable factorization of basis functions . in general , therefore , @xmath162 depends on both form and discretization . to determine this parameter , we look at the re - factorized expression ( as established by algorithm  [ algo : pre - evaluation ] ) , and simply count the terms amenable to pre - evaluation . we demonstrate that the orchestration of sharing elimination and pre - evaluation performed by the transformation algorithm guarantees local optimality ( definition  [ def : mln - quasi - optimality ] ) . the proof re - uses concepts and explanations provided throughout the paper , as well as the terminology introduced in section  [ sec : se - algo ] . [ prop : optimal - approach ] consider a multilinear form comprising a set of monomials @xmath132 , and let @xmath57 be the corresponding finite element integration loop nest . let @xmath58 be the transformation algorithm . let @xmath27 be the set of monomials that , according to @xmath58 , need to be pre - evaluated , and let @xmath166 . assume that the pre - evaluation of different monomials does not result in identical tables . then , @xmath167 is a local optimum in the sense of definition  [ def : mln - quasi - optimality ] and satisfies constraint  [ const : th ] . we first observe that the cost function @xmath131 predicts the _ exact _ gain / loss in monomial pre - evaluation , so @xmath27 and @xmath168 can actually be constructed . let @xmath169 denote the operation count for @xmath57 and let @xmath170 be the subset of innermost loops ( all @xmath51 loops in figure  [ code : loopnest - opt ] ) . we need to show that there is no other synthesis @xmath171 satisfying constraint  [ const : th ] such that @xmath172 . this holds if and only if    1 . _ the coordination of pre - evaluation with sharing elimination is optimal_. this boils down to prove that 1 . _ pre - evaluating any @xmath173 would result in @xmath174 _ 2 . _ not pre - evaluating any @xmath175 would result in @xmath174 _ 2 . _ sharing elimination leads to a ( at least ) local optimum . _ +    we discuss these points separately    1 .   1 . let @xmath176 represent the set of tables resulting from applying pre - evaluation to a monomial @xmath53 . consider two monomials @xmath177 and the respective sets of pre - evaluated tables , @xmath178 and @xmath179 . if @xmath180 , at least one table is assignable to the same temporary . @xmath58 , therefore , may not be optimal , since @xmath131 only distinguishes monomials in `` isolation '' . we neglect this scenario ( see assumptions ) because of its purely pathological nature and its  with high probability negligible impact on the operation count . 2 .   let @xmath181 and @xmath182 be two monomials sharing some generic multilinear symbols . if @xmath123 were carelessly pre - evaluated , there may be a potential gain in sharing elimination that is lost , potentially leading to a non - optimum . this situation is prevented by construction , because @xmath58 exhaustively searches all possible bipartitions on order to determine an optimum which satisfies constraint  [ const : th ] . recall that since the number of monomials is in practice very small , this pass can rapidly be accomplished . 2 .   consider algorithm  [ algo : sharing - elimination ] . proposition  [ prop : multi - struct ] ensures that there are only two ways of scheduling the multilinear operands in @xmath77 : through generalized code motion ( strategy  [ strategy : i ] ) or factorization of multilinear symbols ( via strategy  [ strategy : ii ] ) . if applied , these two strategies would lead , respectively , to performing @xmath80 and @xmath81 multiplications at every loop iteration . since strategy  [ strategy : i ] is applied if and only if @xmath183 and does not change the structure of the expression ( it requires neither expansion nor factorization ) , step ( 1 ) can not prune the optimum from the search space . + after structuring the sharing graph @xmath97 in such a way that only flop - decreasing transformations are possible , the ilp model is instantiated . at this point , proving optimality reduces to establishing the correctness of the model , which is relatively straightforward because of its simplicity . the model aims to minimize the operation count by selecting the most promising factorizations . the second set of constraints is to select all edges ( i.e. , all multiplications ) , exactly once . the first set of inequalities allows multiplications to be scheduled : once a vertex @xmath96 is selected ( i.e. , once a symbol is decided to be factorized ) , all multiplications involving @xmath96 can be grouped .    throughout the paper we have reiterated the claim that algorithm  [ algo : gamma ] achieves a globally optimal flop count if stronger preconditions on the input variational form are satisfied . we state here these preconditions , in increasing order of complexity .    1 . there is a single monomial and only a specific coefficient ( e.g. , the coordinates field ) . this is by far the simplest scenario , which requires no particular transformation at the level of the outer loops , so optimality naturally follows . 2 .   there is a single monomial , but multiple coefficients are present . optimality is achieved if and only if all sub - expressions depending on coefficients are structured ( see section  [ sec : se - rln ] ) . this avoids ambiguity in factorization , which in turn guarantees that the output of step ( 7 ) in algorithm  [ algo : sharing - elimination ] is optimal . 3 .   there are multiple monomials , but either at most one coefficient ( e.g. , the coordinates field ) or multiple coefficients not inducing sharing across different monomials are present . this reduces , respectively , to cases ( 1 ) and ( 2 ) above . there are multiple monomials , and coefficients are shared across monomials . optimality is reached if and only if the coefficient - dependent sub - expressions produced by algorithm  [ algo : sharing - elimination ]  that is , the by - product of factorizing test / trial functions from distinct monomials  preserve structure . sharing elimination and pre - evaluation , as well as the transformation algorithm , have been implemented in coffee , the compiler for finite element integration routines adopted in firedrake . in this section , we briefly discuss the aspects of the compiler that are relevant for this article . coffee implements sharing elimination and pre - evaluation by composing building block transformation operators , which we refer to as _ rewrite operators_. this has several advantages . the first is extensibility . new transformations , such as sum factorization in spectral methods , could be expressed by composing the existing operators , or with small effort building on what is already available . second , generality : coffee can be seen as a lightweight , low level computer algebra system , not necessarily tied to finite element integration . third , robustness : the same operators are exploited , and therefore tested , by different optimization pipelines . the rewrite operators , whose ( python ) implementation is based on manipulation of abstract syntax trees ( asts ) , comprise the coffee language . a non - exhaustive list of such operators includes expansion , factorization , re - association , generalized code motion . coffee aims to be independent of the high level form compiler . it provides an interface to build generic asts and only expects expressions to be in normal form ( or sufficiently close to it ) . for example , firedrake has transitioned from a version of the fenics form compiler @xcite modified to produce asts rather than strings , to a newly written compiler , while continuing to emply coffee . thus , coffee decouples the mathematical manipulation of a form from code optimization ; or , in other words , relieves form compiler developers of the task of fine scale loop optimization of generated code .      for several reasons , basis function tables may be block - sparse ( e.g. , containing zero - valued columns ) . for example , the fenics form compiler implements vector - valued functions by adding blocks of zero - valued columns to the corresponding tabulations ; this extremely simplifies code generation ( particularly , the construction of loop nests ) , but also affects the performance of the generated code due to the execution of `` useless '' flops ( e.g. , operations like ` a + 0 ` ) . in  , a technique to avoid iteration over zero - valued columns based on the use of indirection arrays ( e.g. ` a[b[i ] ] ` , in which ` a ` is a tabulated basis function and ` b ` a map from loop iterations to non - zero columns in a ) was proposed . this technique , however , produces non - contiguous memory loads and stores , which nullify the potential benefits of vectorization . coffee , instead , handles block - sparse basis function tables by restructuring loops in such a manner that low level optimization ( especially vectorization ) is only marginally affected . this is based on symbolic execution of the code , which enables a series of checks on array indices and loop bounds which determine the zero - valued blocks which can be skipped without affecting data alignment . experiments were run on a single core of an intel i7 - 2600 ( sandy bridge ) cpu , running at 3.4ghz , 32 kb l1 cache ( private ) , 256 kb l2 cache ( private ) and 8 mb l3 cache ( shared ) . the intel turbo boost and intel speed step technologies were disabled . the intel ` icc 15.2 ` compiler was used . the compilation flags used were ` -o3 , -xhost ` . the compilation flag ` xhost ` tells the intel compiler to generate efficient code for the underlying platform . the zenodo system was used to archive all packages used to perform the experiments : firedrake @xcite , petsc @xcite , petsc4py @xcite , fiat @xcite , ufl @xcite , ffc @xcite , pyop2 @xcite and coffee @xcite . the experiments can be reproduced using a publicly available benchmark suite  @xcite . we analyze the execution time of four real - world bilinear forms of increasing complexity , which comprise the differential operators that are most common in finite element methods . in particular , we study the mass matrix ( `` ` mass ` '' ) and the bilinear forms arising in a helmholtz equation ( `` ` helmholtz ` '' ) , in an elastic model ( `` ` elasticity ` '' ) , and in a hyperelastic model ( `` ` hyperelasticity ` '' ) . the complete specification of these forms is made publicly available . we evaluate the speed - ups achieved by a wide variety of transformation systems over the `` original '' code produced by the fenics form compiler ( i.e. , no optimizations applied ) . we analyze the following transformation systems :    quad : :    optimized quadrature mode . work presented in  , implemented in in the    fenics form compiler . tens : :    tensor contraction mode . work presented in  , implemented in the fenics    form compiler . auto : :    automatic choice between ` tens ` and ` quad ` driven by heuristic    ( detailed in   and summarized in section  [ sec : qualitative ] ) .    implemented in the fenics form compiler . ufls : :    uflacs , a novel back - end for the fenics form compiler whose primary    goals are improved code generation and execution times . cfo1 : :    generalized loop - invariant code motion . work presented in  ,    implemented in coffee . cfo2 : :    optimal loop nest synthesis with handling of block - sparse tables . work    presented in this article , implemented in coffee . the values that we report are the average of three runs with `` warm cache '' ; that is , with all kernels retrieved directly from the firedrake s cache , so code generation and compilation times are not counted . the timing includes however the cost of both local assembly and matrix insertion , with the latter minimized through the choice of a mesh ( details below ) small enough to fit the l3 cache of the cpu .    for a fair comparison , small patches were written to the make ` quad ` , ` tens ` , and ` ufls ` compatible with firedrake . by executing all simulations in firedrake , we guarantee that both matrix insertion and mesh iteration have a fixed cost , independent of the transformation system employed . the patches adjust the data storage layout to what firedrake expects ( e.g. , by generating an array of pointers instead of a pointer to pointers , by replacing flattened arrays with bi - dimensional ones ) .    for constraint  [ const : th ] , discussed in section  [ sec : mem - const ] , we set @xmath184 ; that is , the size of the processor l2 cache ( the last level of private cache ) . when the threshold had an impact on the transformation process , the experiments were repeated with @xmath185 . the results are documented later , individually for each problem .    following the methodology adopted in  , we vary the following parameters :    * the polynomial degree of test , trial , and coefficient ( or `` pre - multiplying '' ) functions , @xmath186 * the number of coefficient functions @xmath187    while constants of our study are    * the space of test , trial , and coefficient functions : lagrange * the mesh : tetrahedral with a total of 4374 elements * exact numerical quadrature ( we employ the same scheme used in  , based on the gauss - legendre - jacobi rule )          we report the results of our experiments in figures  [ fig : mass ] ,  [ fig : helmholtz ] ,  [ fig : elasticity ] , and  [ fig : hyperelasticity ] as three - dimensional plots . the axes represent @xmath188 , @xmath189 , and code transformation system . we show one subplot for each problem instance @xmath190 , with the code transformation system varying within each subplot . the best variant for each problem instance is given by the tallest bar , which indicates the maximum speed - up over non - transformed code . we note that if a bar or a subplot are missing , then the form compiler failed to generate code because it either exceeded the system memory limit or was otherwise unable to handle the form . the rest of the section is organized as follows : we first provide insights into the general outcome of the experimentation ; we then comment on the impact of a fundamental low - level optimization , namely autovectorization ; finally , we motivate , for each form , the performance results obtained . [ [ high - level - view ] ] high level view + + + + + + + + + + + + + + +    our transformation strategy does not always guarantee minimum execution time . in particular , about 5@xmath191 of the test cases ( 3 out of 56 , without counting marginal differences ) show that ` cfo2 ` was not optimal in terms of runtime . the most significant of such test cases is the elastic model with @xmath192 $ ] . there are two reasons for this . first , low level optimization can have a significant impact on the actual performance . for example , the aggressive loop unrolling in ` tens ` eliminates operations on zeros and reduces the working set size by not storing entire temporaries ; on the other hand , preserving the loop structure can maximize the chances of autovectorization . second , the transformation strategy adopted when @xmath141 is exceeded plays a key role , as we will later elaborate . [ [ autovectorization ] ] autovectorization + + + + + + + + + + + + + + + + +    we chose the mesh dimension and the function spaces such that the inner loop sizes would always be a multiple of the machine vector length . this ensured autovectorization in the majority of code variants . the biggest exception is ` quad ` , due to the presence of indirection arrays in the generated code . in ` tens ` , loop nests are fully unrolled , so the standard loop vectorization is not feasible ; the compiler reports suggest , however , that block vectorization @xcite is often triggered . in ` ufls ` , ` cfo1 ` , and ` cfo2 ` the iteration spaces have identical structure , with loop vectorization being regularly applied .    [ [ mass - matrix ] ] mass matrix + + + + + + + + + + +    we start with the simplest of the bilinear forms investigated , the mass matrix . results are in figure  [ fig : mass ] . we first notice that the lack of improvements when @xmath193 is due to the fact that matrix insertion outweighs local assembly . for @xmath194 , ` cfo2 ` generally shows the highest speed - ups . it is worth noting why ` auto ` does not always select the fastest implementation : ` auto ` always opts for ` tens ` , while for @xmath195 ` quad ` tends to be preferable . on the other hand , ` cfo2 ` always makes the optimal decision about whether to apply pre - evaluation or not . surprisingly , despite the simplicity of the form , the performance of the various code generation systems can differ significantly .    [ [ helmholtz ] ] helmholtz + + + + + + + + +    as in the case of mass matrix , when @xmath193 the matrix insertion phase is dominant . for @xmath194 , the general trend is that ` cfo2 ` outperforms the competitors . in particular :    @xmath196 : :    pre - evaluation makes ` cfo2 ` notably faster than ` cfo1 ` , especially for    high values of @xmath188 ; ` auto ` correctly selects ` tens ` , which    is comparable to ` cfo2 ` . @xmath197 : :    ` auto ` picks ` tens ` ; the choice is however sub - optimal when    @xmath198 and @xmath199 . this can indirectly be    inferred from the large gap between ` cfo2 ` and ` tens / auto ` : ` cfo2 `    applies sharing elimination , but it correctly avoids pre - evaluation    because of the excessive expansion cost . @xmath200 and @xmath201 : :    ` auto ` reverts to ` quad ` , which would theoretically be the right    choice ( the flop count is much lower than in ` tens ` ) ; however , the    generated code suffers from the presence of indirection arrays , which    break autovectorization and `` traditional '' code motion . the slow - downs ( or marginal improvements ) seen in a small number of cases exhibited by ` ufls ` can be attributed to the presence of sharing in the generated code . an interesting experiment we additionally performed was relaxing the memory threshold by setting @xmath185 . we found that this makes ` cfo2 ` generally slower for @xmath195 , with a maximum slow - down of 2.16@xmath202 with @xmath203 . this effect could be worse when running in parallel , since the l3 cache is shared and different threads would end up competing for the same resource .    [ [ elasticity ] ] elasticity + + + + + + + + + +    the results for the elastic model are displayed in figure  [ fig : elasticity ] . the main observation is that ` cfo2 ` never triggers pre - evaluation , although in some occasions it should . to clarify this , consider the test case @xmath204 , in which ` tens / auto ` show a considerable speed - up over ` cfo2 ` . ` cfo2 ` finds pre - evaluation profitable in terms of operation count , although it is eventually not applied to avoid exceeding @xmath141 . however , running the same experiments with @xmath185 resulted in a dramatic improvement , even higher than that obtained by ` tens ` . the reason is that , despite exceeding @xmath141 by roughly 40@xmath191 , the saving in operation count is so large ( 5@xmath202 in this specific problem ) that pre - evaluation would in practice be the winning choice . this suggests that our objective function should be improved to handle the cases in which there is a significant gap between potential cache misses and reduction in operation count . we also note that :    * the differences between ` cfo2 ` and ` cfo1 ` are due to the perfect sharing elimination and the zero - valued blocks avoidance technique presented in section  [ sec : zeros ] . * when @xmath197 , ` auto ` prefers ` tens ` over ` quad ` , which leads to sub - optimal operation counts and execution times . * ` ufls ` often results in better execution times than ` quad ` and ` tens ` . this is due to multiple factors , including avoidance of indirection arrays , preservation of loop structure , and a more effective code motion strategy . [ [ hyperelasticity ] ] hyperelasticity + + + + + + + + + + + + + + +    in the experiments on the hyperelastic model , shown in figure  [ fig : hyperelasticity ] , ` cfo2 ` exhibits the largest gains out of all problem instances considered in this paper . this is a positive result , since it indicates that our transformation algorithm scales well with form complexity . the fact that all code transformation systems ( apart from ` tens ` ) show quite significant speed - ups suggests two points . first , the baseline is highly inefficient . with forms as complex as in the hyperelastic model , a trivial translation of integration routines into code should always be avoided as even the best general - purpose compiler available ( the intel compiler on an intel platform at maximum optimization level ) fails to exploit the structure inherent in the expressions . second , the strategy for removing spatial and temporal sharing has a tremendous impact . sharing elimination as performed by ` cfo2 ` ensures a critical reduction in operation count , which becomes particularly pronounced for higher values of @xmath188 . we have developed a theory for the optimization of finite element integration loop nests . the article details the domain properties which are exploited by our approach ( e.g. , linearity ) and how these translate to transformations at the level of loop nests . all of the algorithms shown in this paper have been implemented in coffee , a compiler publicly available fully integrated with the firedrake framework . the correctness of the transformation algorithm was discussed . the performance results achieved suggest the effectiveness of our methodology . we have defined sharing elimination and pre - evaluation as high level transformations on top of a specific set of rewrite operators , such as code motion and factorization , and we have used them to construct the transformation space . there are three main limitations in this process . first , we do not have a systematic strategy to optimize sub - expressions which are independent of linear loops . although we have a mechanism to determine how much computation should be hoisted to the level of the integration ( reduction ) loop , it is not clear how to effectively improve the heuristics used at step ( 6 ) in algorithm  [ algo : sharing - elimination ] . second , lower operation counts may be found by exploiting domain - specific properties , such as redundancies in basis functions ; this aspect is completely neglected in this article . third , with constraint  [ const : le ] we have limited the applicability of code motion . this constraint was essential given the complexity of the problem tackled . another issue raised by the experimentation concerns selecting a proper threshold for constraint  [ const : th ] . to solve this problem would require a more sophisticated cost model , which is an interesting question deserving further research . we also identify two additional possible research directions : a complete classification of forms for which a global optimum is achieved ; and a generalization of the methodology to other classes of loop nests , for instance those arising in spectral element methods .
we present an algorithm for the optimization of a class of finite element integration loop nests . 
 this algorithm , which exploits fundamental mathematical properties of finite element operators , is proven to achieve a locally optimal operation count . in specified circumstances 
 the optimum achieved is global . 
 extensive numerical experiments demonstrate significant performance improvements over the state of the art in finite element code generation in almost all cases . 
 this validates the effectiveness of the algorithm presented here , and illustrates its limitations . 
 this work was supported by by the department of computing at imperial college london , the engineering and physical sciences research council [ grant number ep / l000407/1 ] , and the natural environment research council [ grant numbers ne / k008951/1 and ne / k006789/1 ] , and by a hipeac collaboration grant . 
 the authors would like to thank dr . 
 andrew t.t . 
 mcrae , dr . lawrence mitchell , and dr . 
 francis russell for their invaluable suggestions and their contribution to the firedrake project . 
 author s addresses : fabio luporini @xmath0 paul h. j. kelly , department of computing , imperial college london ; david a. ham , department of mathematics , imperial college london ;