# An Algebraic Programming Style for Numerical Software and its Optimization

## 1 Introduction

The purpose of the Sophus approach to writing partial differential equation (PDE) solvers originally proposed in  is to close the gap between the underlying coordinate-free mathematical theory and the way actual solvers are written. The main ingredients of Sophus are:

1. A library of abstract datatypes corresponding to manifolds, scalar fields, tensors, and the like, figuring in the abstract mathematical theory.
2. Expressions involving these datatypes written in a side-effect free algebraic style similar to the expressions in the underlying mathematical theory.

Because of the emphasis on abstract datatypes, Sophus is most naturally combined with object-oriented languages or other languages supporting abstract datatypes. Hence, we will be discussing high-performance computing (HPC) optimization issues within an object-oriented or abstract datatype context, using abstractions that are suitable for PDEs.

Sophus is not simply object-oriented scientific programming, but a much more structured approach dictated by the underlying mathematics. The object-oriented numerics paradigm proposed in  is related to Sophus in that it uses abstractions corresponding to familiar mathematical constructs such as tensors and vectors, but these do not include continuous structures such as manifolds and scalar fields. The Sophus approach is more properly called coordinate-free numerics . A fully worked example of conventional vs. coordinate-free programming of a computational fluid dynamics problem (wire coating for Newtonian and non-Newtonian flows) is given in .

Programs in a domain-specific programming style like Sophus may need additional optimization in view of their increased use of expensive constructs. On the other hand the restrictions imposed by the style may lead to new high-level optimization opportunities that can be exploited by dedicated tools. Automatic selection of high-level HPC transformations (especially loop transformations) has been incorporated in the IBM XL Fortran compiler, yielding a performance improvement for entire programs of typically less than 2$`\times `$ \[14, p. 239\]. We hope Sophus style programming will allow high-level transformations to become more effective than this.

In the context of Sophus and object-oriented programming this article focuses on the following example. Object-oriented languages encourage the use of self-mutating (self-updating, mutative) objects rather than a side-effect free algebraic expression style as advocated by Sophus. The benefits of the algebraic style are considerable. We obtained a reduction in source code size using algebraic notation vs. an object-oriented style of up to 30% in selected procedures of a seismic simulation code, with a correspondingly large increase in programmer productivity and maintainability of the code as measured by the Cocomo technique , for instance. On the negative side, the algebraic style requires lots of temporary data space for (often very large) intermediate results to be allocated and subsequently recovered. Using self-mutating objects, on the other hand, places some of the burden of variable management on the programmer and makes the source code much more difficult to write, read, and maintain. It may yield much better efficiency, however. Now, by including certain restrictions as part of the style, a precise relationship between self-mutating notation and algebraic notation may be achieved. Going one step further, we see that the natural way of building a program from high-level abstractions may be in direct conflict with the way current compilers optimize program code. We propose a source-to-source optimization tool, called CodeBoost, as a solution to many of these problems. Some further promising optimization opportunities we have experimented with but not yet included in CodeBoost are also mentioned. The general approach may be useful for other styles and other application domains as well.

This paper is organized as follows. After a brief overview of tensor based abstractions for numerical programming and their realization as a software library (Section 2), we discuss the relationship between algebraic and self-mutating expression notation, and how the former may be transformed into the latter (Section 3). We then discuss the implementation of the CodeBoost source-to-source optimization tool (Section 4), and give some further examples of how software construction using class abstractions may conflict with efficiency issues as well as lead to new opportunities for optimization (Section 5). Finally, we present conclusions and future work (Section 6).

## 2 A Tensor Based Library for Solving PDEs

Historically, the mathematics of PDEs has been approached in two different ways. The solution-oriented approach uses concrete representations of vectors and matrices, discretisation techniques, and numerical algorithms, while the abstract approach develops the theory in terms of manifolds, scalar fields, tensors, and the like, focusing more on the structure of the underlying concepts than on how to calculate with them (see  for a good introduction).

The former approach is the basis for most of the PDE software in existence today. The latter has very promising potential for the structuring of complicated PDE software when combined with template class based programming languages or other languages supporting abstract datatypes. As far as notation is concerned, the abstract mathematics makes heavy use of overloaded infix operators. Hence, user-definable operators and operator overloading are further desirable language features in this application domain. C++  comes closest to meeting these desiderata, but, with modules and user-definable operators, Fortran 90/95  can also be used. In its current form Java  is less suitable. It has neither templates nor user-definable operators. Also, Java’s automatic memory management is not necessarily an advantage in an HPC setting \[16, Section 4\]. Some of these problems may be alleviated in Generic Java . The examples in this article use C++.

### 2.1 The Sophus Library

The Sophus library provides the abstract mathematical concepts from PDE theory as programming entities. Its components are based on the notions of manifold, scalar field and tensor field, while the implementations are based on the conventional numerical algorithms and discretisations. Sophus is currently structured around the following concepts:

* Basic $`n`$-dimensional mesh structures. These are like rank $`n`$ arrays (i.e., with $`n`$ independent indices), but with operations like $`+`$, $``$ and $``$ mapped over all elements (much like Fortran 90/95 array operators) as well as the ability to add, subtract or multiply all elements of the mesh by a scalar in a single operation. There are also operations for shifting meshes in one or more dimensions. Parallel and sequential implementations of mesh structures can be used interchangeably, allowing easy porting between architectures of any program built on top of the mesh abstraction.
* Manifolds. These represent the physical space where the problem to be solved takes place. Currently Sophus only implements subsets of $`R^n`$.
* Scalar fields. These may be treated formally as functions from manifolds to reals, or as arrays indexed by the points of the manifold with reals as data elements. Scalar fields describe the measurable quantities of the physical problem to be solved. As the basic layer of “continuous mathematics” in the library, they provide the partial derivation operations. Also, two scalar fields on the same manifold may be pointwise added, subtracted or multiplied. The different discretisation methods provide different designs for the implementation of scalar fields. A typical implementation would use an appropriate mesh as underlying discrete data structure, use interpolation techniques to give a continuous interface, and map the $`+`$, $``$, and $``$ operations directly to the corresponding mesh operations. In a finite difference implementation partial derivatives are implemented using shifts and arithmetic operations on the mesh.
* Tensors. These are generalizations of vectors and matrices and have scalar fields as components. Tensors define the general differentiation operations based on the partial derivatives of the scalar fields, and also provide operations such as componentwise addition, subtraction and multiplication, as well as tensor composition and application (matrix multiplication and matrix-vector multiplication). A special class are the metric tensors. These satisfy certain mathematical properties, but their greatest importance in this context is that they can be used to define properties of coordinate systems, whether Cartesian, axiosymmetric or curvilinear, allowing partial differential equations to be formulated in a coordinate-free way. The implementation of tensors relies heavily on the arithmetic operations of the scalar field classes.

A partial differential equation in general provides a relationship between spatial derivatives of tensor fields representing physical quantities in a system and their time derivatives. Given constraints in the form of the values of the tensor fields at a specific instance in time together with boundary conditions, the aim of a PDE solver is to show how the physical system will evolve over time, or what state it will converge to if left by itself. Using Sophus, the solvers are formulated on top of the coordinate-free layer, forming an abstract, high level program for the solution of the problem. 

### 2.2 Sophus Style Examples

The algebraic style for function declarations can be seen in Figure 1, which shows specifications of some operations for multidimensional meshes, the lowest level in the Sophus library. The mesh class is parameterized by a class T, so all operations on meshes likewise are parameterized by T. Typical parameters would be a float or scalar field class. The operations declared are defined to behave like pure functions, i.e., they do not update any internal state or modify any of their arguments. Such operations are generally nice to work with and reason about, as their application will not cause any hidden interactions with the environment.

Selected parts of the implementation of a continuous scalar field class are shown in Figure 2. This scalar field represents a multi-dimensional torus, and is implemented using a mesh class as the main data structure. The operations of the class have been implemented as self-mutating operations (Section 3), but are used in an algebraic way for clarity. It is easy to see that the partial derivation operation is implemented by shifting the mesh longer and longer distances, and gradually scaling down the impact these shifts have on the derivative, yielding what is known as a four-point, finite difference, partial derivation algorithm. The addition and multiplication operations are implemented using the element-wise mapped operations of the mesh.

The meshes used in a scalar field tend to be very large. A TorusScalarField may typically contain between 0.2 and 2MB of data, perhaps even more, and a program may contain many such variables. The standard translation technique for a C++ compiler is to generate temporary variables containing intermediate results from subexpressions, adding a considerable run-time overhead to the algebraic style of programming. An implementation in terms of self-mutating operators might yield noticeable efficiency gains. For the addition, subtraction and multiplication algorithms of Figure 2 a self-mutating style is easily obtained. The derivation algorithm will require extensive modification, such as shown in Figure 5, with a marked deterioration in readability and maintainability as a result.

## 3 Algebraic Notation and Self-Mutating Implementation

### 3.1 Self-Mutating Operations

Let a, b and c be meshes with operators as defined in Figure 1. The assignment

```
c = a * 4.0 + b + a

```

is basically evaluated as

```
temp1 = a * 4.0;
temp2 = temp1 + b;
c     = temp2 + a.

```

This involves the creation of the meshes temp1, temp2, c, the first two of which are temporary. Obviously, since all three meshes have the same size and the operations in question are sufficiently simple, repeated use of a single mesh would have been possible in this case. In fact, for predefined types like integers and floats an optimizing C or C++ compiler would translate the expression to a sequence of compound assignments<sup>1</sup><sup>1</sup>1Not to be confused with the C notion of compound statement, which is a sequence of statements enclosed by a pair of braces.

```
c = a; c *= 4.0; c += b; c += a,

```

which repeatedly uses variable c to store intermediate results.

We would like to be able to do a similar optimization of the mesh expression above as well as other expressions involving $`n`$-ary operators or functions of a suitable nature for user-defined types as proposed in . In an object-oriented language, it would be natural to define self-mutating methods (i.e., methods mutating this) for the mesh operations in the above expression. These would be closely similar to the compound assignments for predefined types in C and C++, which return a pointer to the modified data structure. Sophus demands a side-effect free expression style close to the underlying mathematics, however, and forbids direct use of self-mutating operations in expressions. Note that with a self-mutating += operator returning the modified value of its first argument, the expression (a += b) += a would yield $`2(a+b)`$ rather than $`(2a)+b`$.

By allowing the user to define self-mutating operations and providing a way to use them in a purely functional manner, their direct use can be avoided. There are basically two ways to do this, namely, by means of wrapper functions or by program transformation. These will be discussed in the following sections.

### 3.2 Wrapper Functions

Self-mutating implementations can be made available to the programmer in non-self-mutating form by generating appropriate wrapper functions. We developed a C++ preprocessor SCC doing this. It scans the source text for declarations of a standard form and automatically creates wrapper functions for the self-mutating ones. This allows the use of an algebraic style in the program, and relieves the programmer of the burden of having to code the wrappers manually.

A self-mutating operator op= is related to its algebraic analog op by the basic rule

or, if the second argument is the one being updated,<sup>2</sup><sup>2</sup>2This does not apply to built-in compound assignments in C or C++, but user-defined compound assignments in C++ may behave in this way. by the rule

where $``$ denotes equivalence of the left- and right-hand sides, x, y, z are C++ variables, and copy makes a copy of the entire data structure. Now, the Sophus style does not allow aliasing or sharing of objects, and the (overloaded) assignment operator x = y is always given the semantics of x = copy(y) as used in (1) and (2). Hence, in the context of Sophus (1) can be simplified to

| x = y op z;$``$x = y; x op= z; (3) |
| --- |
|  |

and similarly for (2). We note the special case

| x = x op z;$``$x op= z; (4) |
| --- |
|  |

and the obvious generalizations

| x = x op e;$``$x op= e; (5) |
| --- |
|  |
| x = e1 op e2;$``$x = e1; x op= e2; (6) |
|  |

where e, e1, and e2 are expressions. SCC uses rules such as (6) to obtain purely functional behavior from the self-mutating definitions in a straightforward way. Figure 4 shows the wrappers created by SCC for the self-mutating mesh operations of Figure 3. The case of $`n`$-ary operators and functions is similar ($`n1`$). We note that, unlike C and C++ compound assignments, Sophus style self-mutating operators do not return a reference to the variable being updated and cannot be used in expressions. This simpler behavior facilitates their definition in Fortran 90/95 and other languages of interest to Sophus.

The wrapper approach is superficial in that it does not minimize the number of temporaries introduced for expression evaluation as illustrated in Section 3.1. We therefore turn to a more elaborate transformation scheme.

### 3.3 Program Transformation

Transformation of algebraic expressions to self-mutating form with simultaneous minimization of temporaries requires a parse of the program, the collection of declarations of self-mutating operators and functions, and matching them with the types of the operators and functions actually used after any overloading has been resolved. Also, declarations of temporaries have to be added with the proper type. Such a preprocessor would be in a good position to perform other source-to-source optimizations as well. In fact, this second approach is the one implemented in CodeBoost with promising results.

Figure 5 shows an optimized version of the partial derivation operator of class TorusScalarField (Figure 2) that might be obtained in this way. In addition to the transformation to self-mutating form, an obvious rule for ushift was used to incrementalize shifting of the mesh.

Assuming the first argument is the one being updated, some further rules for binary operators used in this stage are

| x op1= e1 op2 e2;$``$ |
| --- |
| {T t = e1; t op2= e2; x op1= t;} (7) |
|  |
| {T t1 = e1; s1;}{T t2 = e2; s2;}$``$ |
| {T t = e1; s1; t = e2; s2;}. (8) |

Here x, t, t1, t2 are variables of type T; e1, e2 are expressions; and self-mutating operators op=, op1=, op2= correspond to operators op, op1, op2, respectively. Recall that Sophus does not allow aliasing. Rule (7) introduces a temporary variable t in a local environment and rule (8) reduces the number of temporary variables by merging two local environments declaring a temporary into a single one.

### 3.4 Benchmarks

#### 3.4.1 Two Kernels

Consider C++ procedures F and P shown in Figure 6. F computes $`x^2+2x`$ using algebraic notation while P computes the same expression in self-mutating form using a single temporary variable temp1. Both were run with meshes of different sizes. The corresponding timing results are shown in Figures 78, and 9.

The mesh size is given in the leftmost column. Mesh elements are single precision reals of 4B each. The second column indicates the benchmark procedure (F or P) or the ratio of the corresponding timings (F/P). The columns NC, NS, OC, and OS give the time in seconds of several iterations over each mesh so that a total of 16 777 216 elements were updated in each case. This corresponds to $`\mathrm{32\hspace{0.17em}768}`$ iterations for mesh size $`8^3`$, $`64`$ iterations for mesh size $`64^3`$, $`1`$ iteration for mesh size $`256^3`$, and so forth. In columns \_C (conventional style) the procedure calls are performed for each element of the mesh, while in columns \_S (Sophus style) they are performed as operations on the entire mesh variables. 

Columns N\_ give the time for unoptimized code (no compiler options), while columns O\_ give the time for code optimized for speed (compiler option -fast for the SUN CC compiler and -Ofast for the Silicon Graphics/Cray CC compiler). The timings represent the median of 5 test runs. These turned out to be relatively stable measurements, except in columns NS and OS, rows $`256^3`$ F and P of Figure 7, where the time for an experiment could vary by more than 100%. This is probably due to paging activity on disk dominating the actual CPU time. Also note that the transformations done by the optimizer are counterproductive in the P case, yielding an NS/OS ratio of 0.8.

When run on the SUN the tests where the only really active processes, while the Cray was run in its normal multi-user mode but at a relatively quiet time of the day (Figure 10). As can be seen the load was moderate (around 58) and although fully utilized, resources where not overloaded.

In the current context, only columns NS and OS are relevant, the other ones are explained in Section 5.1. As expected, the self-mutative form P is a better performer than the algebraic form F when the Sophus style is used. Disregarding the cases with disk paging mentioned above, we see that the self-mutating mesh operations are 1.8–2.4 times faster than their algebraic counterparts, i.e., the CodeBoost transformation roughly doubles the speed of these benchmarks.

#### 3.4.2 Full Application: SeisMod

We also obtained preliminary results on the Silicon Graphics/Cray Origin 2000 for a full application, the seismic simulation code SeisMod, which is written in C++ using the Sophus style. It is a collection of applications using the finite difference method for seismic simulation. Specific versions of SeisMod have been tailored to handle simulations with simple or very complex geophysical properties.<sup>3</sup><sup>3</sup>3SeisMod is used and licensed by the geophysical modelling company UniGEO A.S. (Bergen, Norway). We compared a version of SeisMod implemented using SCC generated wrapper functions and a self-mutating version produced by the CodeBoost source-to-source optimizer:

* The algebraic expression style version turned out to give a 10–30% reduction in source code size and greatly enhanced readability for complicated parts of the code. This implies a significant programmer productivity gain as well as a significant reduction in maintenance cost as measured by the Cocomo technique , for instance
* A 30% speed increase was obtained after 10 selected procedures out of 150 procedures with speedup potential had been brought in self-mutating form. This speedup turned out to be independent of C++ compiler optimization flag settings.

This shows that although a more user-friendly style may give a performance penalty compared to a conventional style, it is possible to regain much of the efficiency loss by using appropriate optimization tools. Also, a more abstract style may yield more cost-effective software, even without these optimizations, if the resulting development and maintenance productivity improvement is taken into account.

## 4 Implementation of CodeBoost

CodeBoost is a dedicated C++ source-to-source transformation tool for Sophus style programs. It has been implemented using the ASF+SDF language prototyping system . ASF+SDF allows the required transformations to be entered directly as conditional rewrite rules whose right- and left-hand sides consist of language (in our case C++) patterns with variables and auxiliary transformation functions. The required language specific parsing, rewriting, and prettyprinting machinery is generated automatically by the system from the high-level specification. Program transformation tools for Prolog and the functional language Clean implemented in ASF+SDF are described in .

An alternative implementation tool would have been the TAMPR program transformation system , which has been used successfully in various HPC applications. We preferred ASF+SDF mainly because of its strong syntactic capabilities enabling us to generate a C++ environment fairly quickly given the complexity of the language.

Another alternative would have been the use of template metaprogramming and/or expression templates . This approach is highly C++ specific, however, and cannot be adapted to Fortran 90/95.

Basically, the ASF+SDF implementation of CodeBoost involves the following two steps:

1. Specify the C++ syntax in SDF, the syntax definition formalism of the system.
2. Specify the required transformation rules as conditional rewrite rules using the C++ syntax, variables, and auxiliary transformation functions.

As far as the first step is concerned, specification of the large C++ syntax in SDF would involve a considerable effort, but fortunately a BNF-like version is available from the ANSI C++ standards committee. We obtained a machine-readable preliminary version  and translated it largely automatically into SDF format. The ASF+SDF language prototyping system then generated a C++ parser from it. The fact that the system accepts general context-free syntax rather then only LALR or other restricted forms of syntax also saved a lot of work in this phase even though the size of the C++ syntax taxed its capabilities.

With the C++ parser in place, the required program transformation rules were entered as conditional rewrite rules. In general, a program transformer has to traverse the syntax tree of the program to collect the context-specific information used by the actual transformations. In our case, the transformer needs to collect the declaration information indicating which of the operations have a self-mutating implementation. Also, in Sophus the self-mutating implementation of an operator (if any) need not update this but can indicate which of the arguments is updated using the upd flag. The transformer therefore needs to collect not only which of the operations have a self-mutating implementation but also which argument is being mutated in each case. As a consequence, CodeBoost has to traverse the program twice: once to collect the declaration information and a second time to perform the actual transformations. Two other issues have to be taken into account:

* C++ programs cannot be parsed before their macros are expanded. Some Sophus-specific language elements are implemented as macros, but are more easily recognized before expansion than after. An example is the upd flag indicating which argument of an operator or function is the one to be updated.
* Compared to the total number of constructs in C++, there are relatively few constructs of interest. Since all constructs have to be traversed, this leads to a plethora of trivial tree traversal rules. As a result, the specification gets cluttered up by traversal rules, making it a lot of work to write as well as hard to understand. One would like to minimize or automatically generate the part of the specification concerned with straightforward program traversal.

Our approach to the above problems is to give the specification a two-phase structure as shown in Figure 11. Under the reasonable assumption that the declarations are not spoiled by macros, the first phase processes the declarations of interest prior to macro expansion using a stripped version of the C++ grammar that captures the declaration syntax only. We actually used a Perl script for this, but it could have been done in ASF+SDF as well. It yields an ASF+SDF specification that is added to the specification of the second phase. The effect of this is that the second phase is specialized for the program at hand in the sense that the transformation rules in the second phase can assume the availability of the declaration information and thus can be specified in a more algebraic, i.e., context independent manner. As a consequence, they are easy to read, consisting simply of the rules for the constructs that may need transformation and using the ASF+SDF system’s built-in innermost tree traversal mechanism. In this way, we circumvented the last-mentioned problem.

As CodeBoost is developed further, it will have to duplicate more and more functions already performed by any C++ preprocessor/compiler. Not only will it have to do parsing (which it is already doing now), but also template expansion, overloading resolution, and dependence analysis. It would be helpful if CodeBoost could tap into an existing compiler at appropriate points rather than redo everything itself. One of the candidates we are considering is the IBM Montana C++ compiler/programming environment , which provides an open architecture with APIs giving access to various compiler intermediate representations with pointers back to the source text.

## 5 Software Structure vs. Efficiency

As noted in Section 1, programs in a domain-specific programming style like Sophus may need additional optimization in view of their increased use of expensive constructs. On the other hand, the restrictions imposed by the style may lead to new high-level optimization opportunities that can be exploited by a CodeBoost-like optimization tool. We give some further examples of both phenomena.

### 5.1 Inefficiencies Caused by the Use of an Abstract Style

We consider an example. As explained in Section 2.1, scalar field operations like $`+`$ and $``$ are implemented on top of mesh operations $`+`$ and $``$. The latter will typically be implemented as iterations over all array elements, performing the appropriate operations pairwise on the elements. For scalar fields, expressions like

$`X_1=A_{1,1}V_1+A_{1,2}V_2`$,   
$`X_2=A_{2,1}V_1+A_{2,2}V_2`$

will force 8 traversals over the mesh data structure. If the underlying meshes are large, this may cause many cache misses for each traversal. Now each of the scalar fields $`A_{i,j}`$, $`V_j`$, and $`X_j`$ are actually implemented using a mesh, i.e., an array of $`n`$ elements, and are represented in the machine by A\[i,j,k\], V\[j,k\] and X\[j,k\] for $`𝚔=1,\mathrm{},𝙺`$, where K is the number of mesh points of the discretisation. In a conventional implementation this would be explicit in the code more or less as follows:

```
for k := 1,K
    for j := 1,2
        X[j,k] := 0
        for i := 1,2
            X[j,k] += A[i,j,k]*V[j,k]
        endfor
    endfor
endfor

```

It would be easy for an optimizer to partition the loops in such a way that the number of cache misses is reduced by a factor of 8.

In the abstract case aggressive in-lining is necessary to expose the actual loop nesting to the optimizer. Even though most existing C++ compilers do in-lining of abstractions, this would be non-trivial since many abstraction layers are involved from the programmer’s notation on top of the library of abstractions down to the actual traversals being performed.

Consider once again the timing results shown in Figure 7, Figure 8, and Figure 9. As was explained in Section 3.4, the procedure calls in columns \_C (conventional style) are performed for each element of the mesh, while they are performed as operations on the entire mesh variables in columns \_S (Sophus style). Columns OS/OC for row P give the relevant figures for the performance loss of optimized Sophus style code relative to optimized conventional style code as a result of Sophus operating at the mesh level rather than at the element level. The benchmarks show a penalty of 1.1–5.3, except for data structures of less than 128kB on the SUN, where a speedup of up to 1.4 (penalty of 0.7) can be seen in Figure 9. As is to be expected, for large data structures the procedure calls in column OC are more efficient than those in column OS, as the optimizer is geared towards improving the conventional kind of code consisting of large loops with procedure calls on small components of data structures. Also, cache and memory misses become very costly when large data structures have to be traversed many times.

The figures for P in column OS of Figure 9 are somewhat unexpected. In these cases OS is the fastest alternative up to a mesh size somewhere between $`32^3`$ and $`64^3`$. This may be due to the smaller number of procedure calls in the OS case than in the OC case. In the latter case F and P are called once per element, i.e., 16 777 216 times, while in the OS case they are called only once and the self-mutating operations are called only 4 times.

Another interesting phenomenon can be seen in column NC of Figure 7 and Figure 8. Here the self-mutating version takes longer than the algebraic version, probably because the compiler automatically puts small temporaries in registers for algebraic expressions, but cannot do so for self-mutating forms. The OC column shows that the optimizer eliminates the difference.

### 5.2 New Opportunities for Optimization

The same abstractions that were a source of worry in the previous section at the same time provide the specificity and typing making the use of high-level optimizations possible. Before they are removed by inlining, the information the abstractions provide can be used to select and apply appropriate datatype specific optimization rules. Sophus allows application of such rules at very high levels of abstraction. Apart from the expression transformation rules (1)–(8) (Section 3), which are applicable to a wide range of operators and functions, further examples at various levels of abstraction are:

* The laws of tensor algebra. In Sophus the tensors contain the continuous scalar fields as elements (Section 2.1), thus making the abstract tensor operations explicit in appropriate modules.
* Specialization of general tensor code for specific coordinate systems. A Cartesian coordinate system gives excellent simplification and axiosymmetric ones also give good simplification compared to general curvilinear code.
* Optimization of operations on matrices with many symmetries. Such symmetries offer opportunities for optimization in many cases, including the seismic modelling application mentioned in Section 3.4.2.

## 6 Conclusions and Future Work

* The Sophus class library in conjunction with the CodeBoost expression transformation tool shows the feasibility of a style of programming PDE solvers that attempts to stay close to the abstract mathematical theory in terms of both the datatypes and the algebraic style of expressions used.
* Our preliminary findings for a full application, the Sophus style seismic simulation code SeisMod, indicate significant programmer productivity gains as a result of adopting the Sophus style.
* There are numerous further opportunities for optimization by CodeBoost in addition to replacement of appropriate operators and functions by their self-mutating versions. Sophus allows datatype specific rules to be applied at very high levels of abstraction.

## Acknowledgments

Hans Munthe-Kaas, André Friis, Kristin Frøysa, Steinar Søreide, and Helge Gunnarsli have contributed to Sophus in various ways.