Peter Fritzson, Vadim Engelson, Krishnamurthy Sheshadri

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.

    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.

    Table 1. Approximate performance comparison for the calcPlot example.

    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.

    Table 2. Performance comparison for the Gauss example.

    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
    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
    kshesh@gmail.com