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:

Directly in the function definition

In a separate command

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.

Integer constant dimension sizes using the array arr:

Symbolic constant dimension sizes:

Unknown dimension sizes with unnamed placeholders:

Unknown dimension sizes with named placeholders:

Unknown dimension sizes that are specified and fixed to the values of integer variables, for example, 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.