*MathCode* is a package that translates a subset of *Mathematica* into a compiled language like Fortran or C++. The chief goal of the design of *MathCode* is to add extra performance and portability to the symbolic prototyping capabilities offered by *Mathematica*. This article discusses several important features of *MathCode,* such as adding type declarations, examples of functions that can be translated, ways to extend the compilable subset, and generating a stand-alone executable, and presents a few application examples.

### Introduction

*MathCode* is a *Mathematica* add-on that translates a *Mathematica* program into C++ or Fortran 90. The subset of *Mathematica* that *MathCode* is able to translate involves purely numerical operations, and no symbolic operations. In the following sections we provide a variety of examples that show precisely what we mean. The code that is generated can be called and run from within *Mathematica*, as if you were running a *Mathematica* function.

There are two important purposes that are served by *MathCode*. Firstly, the C++/Fortran 90 code runs faster, typically by a factor of about a few hundreds (or about 50 to 100) over interpreted (compiled) *Mathematica* code, resulting in considerable performance gains, while still requiring hardly any knowledge of C++/Fortran 90 on the part of the user. Secondly, the generated code can also be executed as a stand-alone program outside *Mathematica*, offering a portability otherwise not possible. You should note, however, that these advantages come at some loss of generality since integer and floating point overflow are not trapped and switched to arbitrary precision as in standard *Mathematica* code. Here the user is responsible for ensuring an appropriate choice of scaling and input data to avoid such problems. The measurements in this article were made using *Mathematica* 6.

There are situations in which having a system such as *MathCode* can be particularly helpful and effective, like when a certain calculation involves a symbolic phase followed by a numerical one. In such a hybrid situation, *Mathematica* can be employed for the symbolic part to give a set of expressions involving only numerical operations that can be made part of a *Mathematica* function, which can then be translated into C++/Fortran 90 using *MathCode*.

In this article, we describe some of the more important features of *MathCode*. For a more detailed discussion the reader is referred to [1]. For brevity, we simply say C++ when we actually mean C++ or Fortran 90: *MathCode* can generate code in both C++ and Fortran, although we illustrate C++ code generation in this article.

In Section 2, we show how to quickly get started with *MathCode* using a simple example of a function to add integers.

Section 3 presents many useful features of *MathCode*. In Section 3.1, we discuss the way the system works, the various auxiliary files generated and what to make of them, and how to build C++ code and install the executable. We then compare the execution times of the interpreted *Mathematica* code and the compiled C++ code. This section also illustrates how *MathCode* works with packages.

Section 3.2 briefly makes a few points about types and type declarations in *MathCode*. There are two ways to declare argument types and return types of a function mentioned in this section.

In Section 3.3, we show how to generate a stand-alone C++ executable. This executable can be run outside of *Mathematica.* We illustrate how to design a suitable main program that the executable runs.

It should be emphasized that *MathCode* can generate C++ for only that subset of *Mathematica* functions referred to as the *compilable subset*. Section 3.4 gives a sample of this subset, while Section 3.5 presents three ways to extend it with the already-available features of *MathCode*: Sections 3.5.1 through 3.5.3 discuss, respectively, symbolic expansion of function bodies, callbacks to *Mathematica*, and handling external functions. Each of these extensions has its own strengths and limitations.

Section 3.6 discusses common subexpression elimination, a feature that is aimed at enhancing the efficiency of generated code.

Section 3.7 presents some shortcuts available in *MathCode* to extract and make assignments to elements of matrices and submatrices, while Section 3.8 is about array declarations.

In Section 4, we present several examples of effectively using *MathCode*. Section 4.1 provides a summary of the examples.

Section 4.2 discusses an essentially simple example, that of computing the function over a grid in the - plane, but done in a somewhat roundabout manner so as to illustrate various features of *MathCode*.

Section 4.3 discusses an implementation of the Gaussian elimination algorithm [2] to solve matrix systems of the type , where is a square matrix of size and (the solution vector) and are vectors of size . In this section, we make a detailed performance study by computing the solution of a matrix system by turning on a few compilation options available in *MathCode*, and also make comparisons with `LinearSolve`.

In Section 4.4, we show how to call external libraries and object files from a C++ program that is automatically generated by *MathCode*. We take the example of a well-known matrix library called SuperLU [3], and demonstrate how to solve, using one of its object modules, a sparse matrix system arising from a partial differential equation.

The *MathCode* User Guide that is available online discusses more advanced aspects, like a detailed account of types and declarations, the numerous options available in *MathCode* with the aid of which the user can control code generation and compilation, and other features. We refer interested readers to [1].

In Section 5, we summarize the salient aspects of *MathCode* and discuss the kinds of applications for which *MathCode* is particularly useful. We conclude the article with a brief summary of various points made. The first version of *MathCode*, released in 1998, was partly developed from the code generator in the ObjectMath environment [4, 5]. The current version is almost completely rewritten and very much improved.

### 2. Getting Started with *MathCode*

#### 2.1. An Example Function

In this section we take the reader on a quick tour of *MathCode* using the simple example of a function to add integers.

The following command loads *MathCode*.

*MathCode* works by generating a set of files in the current directory (see Section 3.1). We can set the directory in the standard way as follows; here, `$MCRoot` is the *MathCode* root directory. The user can, however, use any other directory to store the files.

Let us now define a *Mathematica* function `sumint` to add the first natural numbers.

Note that the body of this function has purely numerical operations, like incrementing the loop index , adding two numbers, and assigning the result to a variable.

#### 2.2. Declaration of Types

We must now declare the data types of the parameter and the local variables and ; we must also specify the return type of the function. We do this using the function `Declare` that *MathCode* provides.

Note that *does not* mean ; the function `Declare` creates an environment in which this is interpreted as a type declaration, that is, an integer variable is being declared in the example. The type `Integer` is translated to a native C int type, and the type `Real` to a native C double type.

#### 2.3. C++ Code

To generate and compile the C++ code, we execute the following command.

Since we have not specified the context of `sumint`, its default context is `Global`. We could, therefore, have simply executed the following command instead.

With the following command, we seamlessly integrate an external program with *Mathematica*.

We can now run the external program in the same way that we would execute a *Mathematica* command.

If we want to run the *Mathematica* code (and not the generated C++ code) for `sumint`, we must first uninstall the C++ executable.

Now the *Mathematica* code for `sumint` will run.

### 3. A Tour of *MathCode*

#### 3.1. How the *MathCode* System Works

*MathCode* works by generating a set of files in the home directory. In the example of `sumint`, the default context is `Global` and the files generated by *MathCode* are: Global.cc (the C++ source file), Global.h and Global.mh (the header files), Globaltm.c, Global.tm and Globalif.cc (the -related files that enable transparently calling C++ versions of the function `sumint` from *Mathematica*), and Globalmain.cc, which contains the function main( ) needed when building a stand-alone executable.

We can also create a package (let us call it `foo`) that defines its own context `foo` instead of the default context `Global`. See Figure 1 for a block diagram of the way the overall system works. The *MathCode* code generator translates the *Mathematica* package to a corresponding C++ source file foo.cc. Additional files are automatically generated: the header file foo.h, the *MathCode* header file foo.mh, the *MathLink-*related files footm.c, foo.tm, foo.icc, and fooif.cc, which enable calling the C++ versions from *Mathematica*, and foomain.cc, which contains the function main that is needed when building a stand-alone executable for `foo` (see Section 3.3). The generated file foo.cc created from the package `foo`, the header file foo.h, and additional files are compiled and linked into two executables. In the case of *MathCode* F90, Fortran 90 is generated and a file foo.f90 is created. No header file is generated in that case since Fortran 90 provides directives for the use of module. External numerical libraries may be included in the linking process by specifying their inclusion (Sections 3.5.3 and 4.5). The executable produced, foo.exe, can be used for stand-alone execution, whereas fooml.exe is used when calling on the compiled C++ functions from *Mathematica* via *MathLink.*

**Figure 1. **Generating C++ code with *MathCode* for a package called `foo`.

Let us see how to work with a package again using the same `sumint` example.

If we are compiling the package `foo` using *MathCode*, we also need to mention `MathCodeContexts` within the path of the package.

We define the function `sumint`

and close the context `foo`.

We next declare the types, and then build and install as before.

Again, since the package `foo` has been defined, it is the default context, and so we could simply have executed the following command.

To run the executable from the notebook, we must install it.

Now the following command runs the C++ executable fooml.exe. The call to `sumint` via *MathLink* is executed 1000 times. The timing measurement includes *MathLink* overhead, which typically for small functions is much more than the execution time for the compiled function. This can be avoided if the loop is executed within the external function itself, as in the example in Section 4.2.5.

Here is the C++ code that was generated.

Note that the function `sumint` appears as foo_Tsumint in the generated code. This is because the full name of the function is in fact `foo`sumint`, and *MathCode* replaces the backquote “” by “`_T`” in the C++ code.

To run the *Mathematica* function (and not its C++ equivalent) `sumint`, we must use the following command to uninstall the C++ code.

Now it is the *Mathematica* code that runs when you execute `sumint`.

You can see that the C++ executable together with the *MathLink* overhead runs about 15 times faster than the *Mathematica* code. The factor by which the performance is enhanced is problem dependent, however. The performance of the *Mathematica* code could also have been improved by using the built-in `Compile` function. In Section 4 we will see many more examples, some quite involved, where we get a range of performance enhancements, also including usage of the `Compile` function.

We clean up the current directory by removing the files automatically generated by *MathCode*.

#### 3.2. Types and Declarations

To be able to generate efficient code, the types of function arguments and return values must be specified, as we have seen in the preceding examples. The basic types used by *MathCode* are

Arrays (vectors and matrices) of these types can also be declared.

Type declarations can be given in two different ways:

The latter construction can be useful if you want to separate already existing *Mathematica* code with the type information needed to be able to generate C++ code using *MathCode*.

#### 3.3. Generating a Stand-Alone Program

So far we have only seen examples in which the installed C++ code can be run within *Mathematica*. However, we can also produce a stand-alone executable. This offers a degree of portability that can be useful in practice.

To illustrate, we take the same example function `sumint` that we discussed in the previous sections. The sequence of commands is very much as in the previous section, except for the option for the *MathCode* function `MakeBinary`, and an appropriate option `MainFileAndFunction` for the function `SetCompilationOptions` immediately after `BeginPackage`. Figure 2 illustrates the process of building the two kinds of executable, namely fooml.exe and foo.exe (on some systems foomain.exe) from a package called `foo`.

**Figure 2. **Building two executables from the package `foo`, possibly including numerical libraries.

The option `MainFileAndFunction` is used to specify the main file. The functions defined in *Mathematica* must have the prefix Global_T (*packagename*_T in general) to be recognized in the main file.

Now we are ready to generate and compile the C++ code for the package `foo`. We can do this in two ways: we can either employ the *MathCode* function `BuildCode`, as in the previous examples, or first execute `CompilePackage` (which generates the C++ source and header files) and then the function `MakeBinary` (which creates the executable).

The last command generates the stand-alone executable foo.exe that can be executed from a command line, or, alternatively, by using the *Mathematica* function `Run`.

If you desire, you can, in addition to the stand-alone executable foo.exe, also generate fooml.exe that can be run from within *Mathematica*, just like before.

Now the following command runs the C++ program foo.cc.

##### 3.3.1. Generating a DLL

Here we briefly mention the possibility of generating a DLL, without giving a full example. To generate a DLL from a package, you have to write a file containing one simple wrapper function in order to make a generated function visible outside the DLL. You write a wrapper function for each generated function. The flags used are as follows:

Here “ext.cpp” is a C++ file with wrapper functions, and “” is a flag for the Visual C++ linker. For other C++ compilers this procedure is not automatic and requires several operating system commands, but the wrapper functions are not needed.

#### 3.4. The Compilable Subset

*MathCode* generates C++ code for a subset of *Mathematica* functions, called the *compilable subset*. The following items give a sample of the compilable subset. For a complete list of *Mathematica* functions in the compilable subset, see [1].

- Statically typed functions, where the types of function arguments and return values are given by the types discussed in Section 3.2
- Scoping constructs: , , , ,
- Procedural constructs: , , , ,
- Lists and tables: , , , ,
- Size functions: ,
- Arithmetic and logical expressions, for example:
`+`,`-`,`*`,`/`,`==`,`!=`,`>`,`!`,`&&`,`||`, and so forth - Elementary functions and some others, for example: , , , , , , , ,
- Constants:
`True`,`False`,`E`,`Pi` - Assignments:
`:=`,`=` - Functional commands: ,
- Some special commands: ,

Functions not in the compilable subset can be used in external code by callbacks to *Mathematica* (see Section 3.5.2 for an example).

Examples of functions that are not a part of the compilable subset include: , , , , , .

These functions can be used if *Mathematica* can evaluate them at compile time to expressions that belong to the compilable subset. In general, *Mathematica* functions that perform symbolic operations are not in the compilable subset. Also, many functions in the subset are implemented with limitations, that is, more difficult cases are not always supported. However, *MathCode* currently provides several ways to extend the compilable subset, as we discuss in the next section.

#### 3.5 Ways to Extend the Compilable Subset

##### 3.5.1. Symbolic Expansion of Function Bodies

Functions not entirely written using *Mathematica* code in the compilable subset, but whose definitions can be evaluated symbolically to expressions that belong to the compilable subset, can be handled by *MathCode*.

Generate C++ code and compile it to an executable file.

The option `EvaluateFunctions` tells *MathCode* to let *Mathematica* expand the function body as much as possible. Everything works fine because the result belongs to the compilable subset.

The generated executable is connected to *Mathematica*:

##### 3.5.2. Callbacks to *Mathematica*

Consider the following function whose definition includes the `Zeta` function, which does not belong to the compilable subset.

Let us plot the function:

We now make the declarations:

These declare statements do not change the way *Mathematica* computes the function.

Let us now generate C++ code and compile it to an executable file. The option `CallBackFunctions` tells *MathCode *which functions have to be evaluated by *Mathematica*. As a result, although the function `Zeta` is not in the compilable subset, an executable is still generated and communicates with the kernel to evaluate `Zeta`.

The generated executable is connected to *Mathematica*:

Now it is the external code that is used to compute the function:

In this case the external code calls *Mathematica* when the `Zeta` function has to be evaluated. After the evaluation the computation proceeds in the external code.

Note that it is the installed code for the function `f` that is executed above, and not the original *Mathematica* function. In the installed code, the argument of `f` must be real, according to our declaration. As a result, , in which we pass a rational number as an argument, is left unevaluated.

We again plot the function, but this time using the external code to evaluate it:

##### 3.5.3. External Functions

We can have references to external objects in C++ code generated by *MathCode*. Let us consider three very simple external functions that compute , , and to illustrate the idea. These must be defined as follows in an external source file that must be in the working directory.

Observe here that each function definition, which is in C language syntax, is followed by a “wrapper” that enables *MathCode* to recognize the object as external. We can then create an object file corresponding to these functions and link the object as follows.

We define a function to create a list of numbers using the external functions.

We now compile the package. Since this is a very small example, we do not bother to create a special package for the code.

Let us now create the *MathLink* binary; to do this when there are external functions, we must specify the option `NeedsExternalObjectModule` as follows.

Here, as we noted above, `external1` and `external2` represent the external object modules external1.o and external2.o. Install the *MathCode*-compiled code so it is called using *MathLink*.

When we make the following plot, it is the external code for `extsqr`, `extexp`, and `extoscillation` that is used.

#### 3.6. Common Subexpression Elimination

Consider the following function whose definition contains a number of common subexpressions (e.g., and ).

There are very efficient algorithms to evaluate functions containing common subexpressions. The basic idea is to evaluate common subexpressions only once and put the results in temporary variables.

Now we generate C++ code using *MathCode* and run it.

*MathCode* does common subexpression elimination (CSE) when the option `EvaluateFunctions` is given to or . This basic strategy could be further improved for special cases in future versions of *MathCode*. Moreover, since mathematical expressions are intrinsically free of side effects and do not have a specific evaluation order, the CSE optimization may change the order of computing subexpressions if this improves performance. Changing the order can sometimes have a small influence on the result when floating-point arithmetic is used.

We take a look at the generated C++ file.

Note how the computation of the function has been divided into small subexpressions that are evaluated only once and then stored in temporary variables for future use. This gives very efficient code for large functions. The speed enhancement of roughly 150% brought about by CSE in this example is not appreciable because the example itself is rather small.

#### 3.7. Extended Matrix Operations

When dealing with matrices, it is very convenient to have a short notation for part extraction. *MathCode* extends the functionality of or to achieve this.

Consider the following matrix:

We can extract rows 2 to 4 as follows, with the shorthand available in *MathCode*.

We can extract the elements in all rows that belong to column 3 and higher:

We can assign values to a submatrix of `A`.

All these operations belong to the compilable subset and can result in compact code. Note: denotes the same *Mathematica* computation as , and is equivalent to .

#### 3.8. Array Declaration and Dimension

In this subsection, we give a few examples of array declarations. There are two main cases to consider.

- Arrays that are passed as function parameters or returned as function values, where the actual array size has been previously allocated
- Declaration of array variables, usually specifying both the type and the allocation of the declared array

There are five allowed ways to specify array dimension sizes in array types for function arguments and results.

- Integer constant dimension sizes, for example:
- Symbolic-constant dimension sizes, for example:
- Unknown dimension sizes with unnamed placeholders, for example:
- Unknown dimension sizes with named placeholders, for example:
- Unknown dimension sizes with variables as dimension sizes, for example:

The dimension sizes can be constant, in which case the size information is part of the type. Alternatively, the sizes are unknown and thus fixed later at runtime when the array is allocated. Such unknown dimension sizes are specified through named (e.g., `n_`) or unnamed (`_`) placeholders.

All arrays that are passed as arguments to functions have already been allocated at runtime. Thus, their sizes are already determined. These sizes might, however, be different for different calls. Therefore it is not allowed to specify conflicting dimension sizes through integer variables (e.g., ) in array types of function parameters or results, as can be done for ordinary declared variables. Only constants and named, or unnamed, placeholders are allowed.

We now give examples of the five different ways of specifying array dimension information in variable declarations. The examples show a global variable declaration using `Declare`, but the same kinds of declarations can also be used for local declarations in functions.

The fifth case is where sizes are specified through integer variables. This is needed to handle declaration and allocation of arrays for which the sizes are not determined until runtime.

`arr`:

`n`,

`m`(e.g., function parameters, local or global variables that are visible from the declaration):

Integer variables, such as `n` and `m`, are assumed to be assigned once; that is, their values are not changed after the initial assignment, so that the declared sizes of allocated arrays are kept consistent with the values of those variables. This single-assignment property is not checked by the current version of the system, however. Thus, the user is responsible for maintaining such consistency.

### 4. Application Examples

#### 4.1. Summary of Examples

In the following we present a few complete application examples using *MathCode*. The first example application is a small *Mathematica* program called `SinSurface` (Section 4.2), which has been designed to illustrate two basic modes of the code generator: compiling without symbolic evaluation (the default mode, in which the function body is translated into C++ as it is), and compilation preceded by symbolic expansion, which is indicated by setting the option (the function body is expanded using symbolic operations, simplified, and then translated).

The second example, presented in Section 4.3, is an implementation of the Gaussian elimination procedure to solve a linear algebraic system of equations (see any standard text on numerical techniques for a discussion of the procedure, e.g., [2]). Here we compile generated C++ code with various options and do a detailed performance analysis.

In Section 4.4, we discuss the example of SuperLU, an external library [3] that performs efficient sparse matrix operations. We give an example of a program useful in solving partial differential equations that calls the SuperLU library and some of its object modules to solve a matrix equation of the type , where is a very sparse square matrix.

#### 4.2. The SinSurface Application Example

Here we describe the `SinSurface` program example. The actual computation is performed by the functions `calcPlot`, `sinFun2`, and their helper functions. The two functions `calcPlot` and `sinFun2` in the `SinSurface` package will be translated into C++ and are declared together with a global array `xyMatrix`.

The array `xyMatrix` represents a grid on which the numerical function `sinFun2` will be computed. The function `calcPlot` accepts five arguments: four of these are coordinates describing a square in the - plane and one is a counter (`iter`) to make the function repeat the computation as many times as necessary in order to measure execution time. For each point on a grid, the numeric function `sinFun2` is called to compute a value that is stored as an element in the matrix representing the grid.

##### 4.2.1. Introduction

The `SinSurface` example application computes a function (here `sinFun2`) over a two-dimensional grid. The function values are stored in the matrix `xyMatrix`. The execution of compiled C++ code for the function `sinFun2` is over 500 times faster than evaluating the same function interpretively within *Mathematica*.

The function `sinFun2` computes essentially the same values as , but in a more complicated way, using a rather large expression obtained through converting the arguments into polar coordinates (through `ArcTan`) and then using a series expansion of both `Sin` and `Cos`, up to 10 terms. The resulting large symbolic expression (more than a page) becomes the body of `sinFun2`, and is then used as input to `CompileEvaluateFunction` to generate efficient C++ code. The symbolic expression and the call to `CompileEvaluateFunction` is initiated by using the `EvaluateFunctions` option.

##### 4.2.2. Initialization

We first set the directory in which *MathCode* will store the auxiliary files, the C++ code, and executable, and then load *MathCode*.

The `SinSurface` package starts in the usual way with a `BeginPackage` declaration that references other packages. `MathCodeContexts` is needed in order to call the code generation related functions.

Next we define possibly exported symbols. Even though it is not necessary here, we enclose these names within … as a kind of context bracket, since this can be put into a cell, which can be conveniently re-evaluated by itself if new names are added to the list.

Now we set compilation options as follows. This defines how the functions and variables in the package should be compiled to C++. By default, all typed variables and functions are compiled. However, the compilation process can be controlled in a more detailed way by giving compilation options to `CompilePackage` or via `SetCompilationOptions`. For example, in this package the function `sinFun2` should be symbolically evaluated before being translated to code, because it contains symbolic operations; the functions `sin`, `cos`, and `arcTan` should not be compiled at all, because they are expanded within the body of `sinFun2`. The remaining typed function, `calcPlot`, will be compiled in the normal way.

##### 4.2.3. The Body of the SinSurface Package

We begin the implementation section of the `SinSurface` package, where functions are defined. This is usually private, to avoid accidental name shadowing due to identical local variables in several packages.

Declare public global variables and private package-global variables:

Taylor-expanded `sin` and `cos` functions called by `sinFun2` are now defined, just for the sake of the example, even though such a series gives lower relative accuracy close to zero. A substitution of the symbol `z` for the actual parameter `x` is necessary to force the series expansion before replacing with the actual parameter.

Define `arcTan`, which converts a grid point to an angle, called by `sinFun2`:

`sinFun2` is the function to be computed and plotted, called by `calcPlot`. It provides a computationally heavy (series expansion) and complicated way of calculating an approximation to . This gives an example of a combination of symbolic and numeric operations as well as a rather standard mix of arithmetic operations. The expanded symbolic expression, which comprises the body of `sinFun2`, is about two pages long when printed.

Note that the types of local variables to `sinFun2` need not be declared, since setting the `EvaluateFunctions` option will make the whole function body be symbolically expanded before translation.

Note also that a function should be without side effects in order to be symbolically expanded before final code generation. For example, there should be no assignments to global variables or input/output, since the relative order of these actions when executing the code often changes when the symbolic expression is created and later rearranged and optimized by the code generator.

The function `calcPlot` calculates data for a plot of `sinFun2` over a grid, which is returned as a array.

##### 4.2.4. Execution

We first execute the application interpretively within *Mathematica*, and then use `Compile` on the key function and execute the application again. Then we compile the application to C++, build an executable, and call the same functions from *Mathematica* via *MathLink*.

Let us first do the *Mathematica* evaluation and plot.

Next, we redefine `sinFun2` to become a compiled version, using *Mathematica*’s standard `Compile`.

##### 4.2.5. Using the MathCode Code Generator

Compile the `SinSurface` package.

The warnings concern local variables in `sinFun2` that have no type information. This is not important because those variables disappear upon symbolic expansion.

The command `MakeBinary` compiles the generated code using a compiler (g++ in the present case). The object code is by default linked into the executable SinSurfaceml.exe for calling the compiled code via *MathLink*.

If any problems are encountered during code compilation, then warning and error messages are shown. Otherwise no messages are shown. When `MakeBinary` is called without arguments, the call applies to the current package.

The command `InstallCode` installs and connects the external process containing the compiled and linked `SinSurface` code.

Execute the generated C++ code for `calcPlot`.

Since the external computation was performed 3000 times, the time needed for one external computation is

Check that the result appears graphically the same.

##### 4.2.6. Performance Comparison

Let us now compare the running times for the three cases, the standard *Mathematica*, compiled *Mathematica*, and the generated C++ code.

The performance between the three forms of execution are compared in Table 1. The generated C++ code for this example is roughly 100 times faster than standard interpreted *Mathematica* code, and 50 times faster than code compiled by the internal *Mathematica* `Compile` command. This is on a Toshiba Satellite-2100, 400 Mhz AMD-K6, running Windows XP Pro SP2 and *Mathematica* 6, without inline and norange optimization. If the inline is specified, the inline directive is passed to the C++ compiler for all functions to be compiled. If norange is specified, array element index range checking is turned off in the code generated by the C++ compiler, resulting in faster but less safe code.

We should emphasize that the comparisons in Table 1 are rather crude for several reasons. From a separate measurement, the loop part of `calcPlot` excluding the call to `sinFun2` comprises 25% of the total `calcPlot` time executed in interpreted *Mathematica*. The `calcPlot` function itself cannot be compiled using `Compile`, since it contains an assignment to a global matrix variable that cannot currently be handled by `Compile`. This might be regarded as unfair to `Compile`. On the other hand, a *MathLink* overhead (divided by 500) in returning the matrix is embedded in the figure for *MathCode*, which can be regarded as unfair to *MathCode*. A better comparison for another small application example is available in Section 4.3.6.

#### 4.3. Gauss Application Example

##### 4.3.1. Introduction

In this section, we present a textbook algorithm, Gaussian elimination (e.g., [2]), to solve a linear equation system. The given linear system, represented by a matrix equation of the type , is subjected to a sequence of transformations involving a pivot, resulting in the solution to the system, contained in the matrix .

The following subsections illustrate the various aspects of the application.

##### 4.3.2. Initialization

Define exported symbols:

##### 4.3.3. Body of the Package

We now define the function `GaussSolveArraySlice`, based on the Gaussian elimination algorithm.

This function accepts three arguments in an attempt to solve a matrix equation of the form . The first two arguments are essentially the matrices and . The third argument specifies the number of times the body of the function must run; this is useful for an accurate measurement of the running time. The function output () has the same shape as the second argument ().

##### 4.3.4. Mathematica Execution

Let us create two random matrices.

In the following, specifies the number of times the body of `GaussSolveArraySlice` runs. The appropriate value of `loops` for reliable estimates of running time is system dependent. A reasonable value of `factor` for a 1.5 GHz computer is about 10. The output checks that the solution obtained is correct.

##### 4.3.5. Generating and Running the C++ Code

The command `BuildCode` translates the package and produces an executable.

Interpreted versions are removed, and compiled ones are used instead.

We now make two runs of the C++ code for the package `Gauss`. The first run evaluates the body of `GaussSolveArraySlice` `loops` times, and returns the solution only once. The second run evaluates the body of `GaussSolveArraySlice` only once, but does this inside a `Do`-loop for `loops` times, returning the solution `loops` times as a result. Clearly, there is overhead in the second run, and the time taken is expected to be higher, as can be seen from the following.

We now also solve the same system using `LinearSolve`.

##### 4.3.6. Performance Comparison

We present the performance analysis in Table 2. As we observe from the table, a performance enhancement by a factor of approximately 500 can be obtained for the compiled C++ code over interpreted *Mathematica*. More importantly, we are able to get a performance close to `LinearSolve`, although we have implemented a simple version of a Gaussian elimination algorithm directly from a textbook as straight-line code without any attempts at tuning or optimization. Also note that `LinearSolve` is more general in that it can also handle sparse arrays efficiently using the *SuperLU* package linked into the *Mathematica* kernel. To achieve similar generality with the *MathCode* package, the example would need to be extended and a *SuperLU* routine, for example, called as an external function from the generated code. In general, it is better to use already implemented robust and reliable routines from packages like *LAPACK* and *SuperLU*, which also can be called as external functions from *MathCode-*generated code. The Gauss example in this article is not intended to replace such routines but to be a simple example of using *MathCode*.

#### 4.4. External Libraries and Functions

We now demonstrate how to call external functions and libraries using *MathCode*. We have already presented an example of how to do this for three very simple functions, , , and , in Section 3.5.3. In this section we present a more realistic application example that illustrates how to employ an external library for handling sparse matrix systems that arise in the solution of partial differential equations [6].

We take as our example the problem of solving the one-dimensional diffusion equation using the method of finite differences.

In this method, the continuous domain is approximated by a set of discrete points called a *grid*, and each derivative is replaced by a certain linear function of values of the dependent variables, called a *finite difference*. For the previous equation, a variant of this method gives

where now and are assumed to take integer values, and and are step sizes along and directions, respectively. The algebraic equation must be solved at each grid point, thus resulting in a simultaneous system of equations, which is essentially a matrix system of the form . Since the matrix system in this case is very sparse, we solve it using the sparse matrix library called SuperLU [3].

The rest of this section assumes that the SuperLU library has been compiled. We now explain how to call the external objects based on this library using *MathCode*.

Here is the *Mathematica* code to solve the one-dimensional diffusion equation.

Note that this source file is somewhat different from the one in Section 3.5.3, mainly because arrays are involved here. This wrapper function makes a reference to a C function linsolve( ) that is defined in the following source file.

It is the function linsolve( ) that solves the matrix equation by calling other object modules of the SuperLU library; from these two C/C++ source codes, object files must be generated using suitable makefiles.

The matrices are expected to be in a special format called “column-compressed storage format,” so as to minimize storage space. Thus, the matrix elements of need not all be specified, since only a small number, , of them are nonzero; here is the number of spatial grid points. The matrix is specified through three row matrices `amat` and `asubmat` (that have a length ), and `xamat` (that has a length ). Our function takes these integers , , and as parameters; in addition, we must pass as parameters the time step and the solution vector of the PDE at time ; the function then returns the solution vector at time .

The function `linsolvepp` must now be defined as an external procedure using the following command.

Note the keyword `InOut` preceding the last argument of `ExternalProcedure`: in the calling function `SolveDiffusion1D`, the array `rhsmat` is passed to `linsolvepp` as input, but `linsolvepp` also returns the solution vector by destroying `rhsmat` and using it to store the solution vector. As a result, the array `rhsmat` is both an input and an output. The way to declare this is by using the keyword `InOut`.

We next declare the types, and then build and install.

We now create the executable. Note that an additional option `NeedsExternalLibrary` must also be specified in this example, since the external objects depend on other objects of the SuperLU library.

We take the following initial conditions.

Now the following command runs the C++ executable fooml.exe. We evolve from to with .

### 5. Summary and Conclusions

*MathCode* is an application package that generates optimized Fortran/C++ code for numerical computations. The code can be either compiled and run from within a notebook environment, or ported, and typically runs several hundred times faster than original *Mathematica* code.

*MathCode* is easy to use, since only the following three simple steps are involved for most applications:

- Add type declarations.
- Execute to generate C++ code and an executable program.
- Execute to connect the executable program to
*Mathematica*.

It must be remembered that only a subset of *Mathematica *functions and operations are translated into C++ by *MathCode*. However, *MathCode* also provides these ways to extend the subset:

- Symbolic evaluation
- Callbacks to
*Mathematica* - Use of external code

To conclude, we remark that *MathCode* can turn *Mathematica* into a powerful environment for prototyping advanced numerical algorithms and production code development. Since it can generate stand-alone code, applications that use *Mathematica* as an environment for development and need to automatically generate efficient C++ code as embedded code in large software systems can greatly benefit.

*MathCode* is a product available both for purchase and free trial (see the website of MathCore Engineering, [1]). Currently, both the C++ and Fortran 90 versions of the code generator are available.

### References

[1] | P. Fritzson, MathCode C++, Linköping, Sweden: MathCore Engineering AB, 1998 www.mathcore.com. |

[2] | C. F. Gerald and P. O. Wheatley, Applied Numerical Analysis, 5th ed., Reading, MA: Addison-Wesley, 1994. |

[3] | The SuperLU package is available for download at crd.lbl.gov/~xiaoye/SuperLU. |

[4] | L. Viklund, J. Herber, and P. Fritzson. “The implementation of ObjectMath—A High-Level Programming Environment for Scientific Computing,” in Compiler Construction, Proceedings of the Fourth International Workshop on Compiler Construction (CC 1992), Paderborn, Germany (U. Kastens and P. Pfahler, eds.), Lecture Notes in Computer Science, 641, London: Springer-Verlag, 1992 pp. 312-318. Also see the ObjectMath home page, www.ida.liu.se/~pelab/omath. |

[5] | P. Fritzson, V. Engelson, and L. Viklund, “Variant Handling, Inheritance and Composition in the ObjectMath Computer Algebra Environment,” in Proceedings of the International Symposium on Design and Implementation of Symbolic Computaton Systems (DISCO 1993), Gmunden, Austria (A. Miola, ed.), Lecture Notes in Computer Science, 722, London: Springer-Verlag, 1993 pp. 145-160. Also see the ObjectMath home page, www.ida.liu.se/~pelab/omath. |

[6] | K. Sheshadri and P. Fritzson, “MathPDE: A Package to Solve PDEs,” submitted to The Mathematica Journal, 2005. |

P. Fritzson, V. Engelson, and K. Sheshadri, “MathCode: A System for C++ or Fortran Code Generation from Mathematica,” The Mathematica Journal, 2011.dx.doi.org/10.3888/tmj.10.4-7. |

### About the Authors

Peter Fritzson, Ph.D., is Director of research at MathCore Engineering AB. He is also a full professor at Linköping University and Director of the Programming Environment Laboratory (PELAB) at the Department of Computer and Information Science, Linköping University, Sweden. He initiated the development of *MathCode* and the ObjectMath environment, and developed parts of the first version of *MathCode*. Professor Fritzson is regarded as one of the leading authorities on the subject of object-oriented mathematical modeling languages, and has recently published a book, *Principles of Object-Oriented Modeling and Simulation with Modelica 2.1*, Wiley-IEEE Press. He currently holds positions as Chairman of the Scandinavian Simulation Society, Secretary of EuroSim, and Vice Chairman of the Modelica Association, an organization he helped to establish. Professor Fritzson has published 10 books/proceedings and approximately 170 scientific papers.

Vadim Engelson, Ph.D., is a development engineer at MathCore Engineering AB. He was previously an assistant professor at Linköping University. Dr. Engelson received his Ph.D. in computer science at Linköping University in 2000 for his research in tools for design, interactive simulation, and visualization of object-oriented models in scientific computing. He has participated in several Swedish and European projects. Dr. Engelson is Technical Coordinator of the Scandinavian Simulation Society. He developed most of the *MathCode* *Mathematica* to C++ translator and several parts of *MathModelica*. He is the author of 21 publications.

**Krishnamurthy**** **Sheshadri, Ph.D., is a development engineer at *Connexios Life Sciences Private Limited*. He was previously a post-doctoral student at PELAB, Linköping University, where he used the *MathCode* system for advanced applications. Work on the project was done when Sheshadri was associated with Modeling and Simulation Research, Bangalore, India and Linköping University.

**Peter Fritzson
Vadim Engelson
**

*Mathcore Engineering AB*

Teknikringen 1B

SE 583 30 Linköping, Sweden

Teknikringen 1B

SE 583 30 Linköping, Sweden

*peter.fritzson@mathcore.com*

*vadim.engelson@mathcore.com*

**Krishnamurthy Sheshadri
**

*Connexios Life Sciences Private Limited*

49, Shilpa Vidya

First Main Road, J P Nagar 3rd Phase

Bangalor-560 078, India

49, Shilpa Vidya

First Main Road, J P Nagar 3rd Phase

Bangalor-560 078, India

*kshesh@gmail.com*