MathCode C++

A Basic Introduction and Overview of MathCode

by Mats Jirstrand
1999 Mathematica Developer Conference,
21-23 October, Champaign

What is MathCode C++?

MathCode C++ is an application package for Mathematica that generates optimized C++ code for numerical computations.

Compiled code can easily be called from within Mathematica.

Stand-alone code can be generated (i.e. code that can be executed independent of Mathematica).

Increase in performance up to 1000 times for some computations.

External libraries of C, C++, and Fortran77 code can easily be connected and called from within Mathematica.

MathCode C++ turns Mathematica into a powerful environment for prototyping of advanced numerical algorithms for production code development.

Who can benefit from MathCode C++?

Users that are using Mathematica as an environment for developing applications and then want to automatically generate efficient C++ code for usage in large software systems: generation of embedded code.

Users that would like to use Mathematica as an interface to old function libraries or single programs written in Fortran, C, or C++.

Users that would like to write Mathematica packages for numerical computations and then compile them to executables for proprietary reasons and/or for easy distribution.

Users that would like to take advantage of the excellent prototyping environment in Mathematica but also have easy access to special purpose routines already written in C, C++, or Fortran.

Users that have performance problems e.g. due to very large data files or functions that are computationally expensive to compute.

Basic Usage

Load the package

Needs["MathCode`"]

Since a number of temporary files will be generated by MathCode C++ we create a temporary directory and set this to the current working directory. However, first we change the current directory to the same as the one that this notebook belongs to.

CreateDirectory["Tmp"];
SetDirectory[ToFileName[{".", "Tmp"}]];

Define a simple function but add type information about argument and return value (blue color)

f[Real x_]->Real:=Sin[x]/(1 + (x^2))

The function can be used a usual

Plot[f[x], {x, -π, π}];

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

BuildCode[]

The generated executable is connected to Mathematica

InstallCode[];

Now, the executable file is used in computations instead

Plot[f[x], {x, -π, π}];

We take a look at the generated C++ file

Run["Notepad Global.cc"];

Uninstall the executable, remove the temporary files, and remove the function f

UninstallCode[];
CleanMathCodeFiles[Confirm -> False, 
        CleanAllBut -> {}];
Remove[f];

Hence, at the very basic level the MathCode user only have to

Add type information (function argument and return value)

Build the C++ code using BuildCode[]

Install the executable file using InstallCode[]

The Compilable Subset

MathCode C++ generates efficient C++ code for a subset of Mathematica.

A necessary condition for using MathCode C++ to generate code for a Mathematica function is that the input and output are numerical entities such as vectors or matrices of integers or real numbers.

Code for intermediate symbolic computations can be generated if Mathematica can evaluate these computations at compile time.

The subset of the Mathematica language that MathCode C++ can handle is denotet the "compilable subset". The following list gives a brief overview and characterization of this subset

Statically typed functions where the types of function arguments and return values are given by the above types.

Scoping constructs:
	BeginPackage[], EndPackage[],
	Module[], Block[], With[]

Procedural constructs:
	For[], While[], If[], Which[], Do[]

Lists and tables:
	List[], Table[], Array[], Range[], IdentityMatrix[]

Size functions:
	Dimensions[], Length[]

Arithmetic and logic expressions e.g.:	
	+, -, *, /, ==, !=, >, !, &&, ||, ...

Elementary functions e.g.:
	Sin[], Exp[], ArcSin[], Sqrt[], ...
and others e.g.:
	Round[], Max[], Cross[], Transpose[], Dot[], ...

Constants:
	True, False, E, Pi

Assignments:
	:=, =

Functional commands:
	Map[], Apply[]

Some special commands:
	Sum[], Product[]

A complete listing of the compilable subset can be found in the MathCode C++ manual or at our homepage www.mathcore.com.

Functions not in the compilable subset can be used in external code by call-backs to Mathematica, see below for an example.

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

Integrate[], Solve[], FindRoot[], LinearSolve[], Expand[], Factor[]

These function can be used if Mathematica can evaluate them at compile time to expressions that belong to the "compilable subset".

Enlargement of the compilable subset

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 belongs to the compilable subset, can be handled by MathCode C++.

f[Real a_,  Real b_] -> Real := 
      Integrate[x Sin[x], {x, a, b}]

f[1. , 2. ]

?f

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

BuildCode[EvaluateFunctions -> {f}]

The option EvaluateFunctions tells MathCode C++ to let Mathematica expand the function body as much as possible. Everything works fine since the result belongs to the compilable subset

Integrate[x Sin[x], {x, a, b}]

The generated executable is connected to Mathematica

InstallCode[];

f[1. , 2. ]

Uninstall the executable, remove the temporary files, and remove the function f

UninstallCode[];
CleanMathCodeFiles[Confirm -> False, 
        CleanAllBut -> {}];
Remove[f];

Call-backs to Mathematica

Consider the following function whose definition includes the Zeta function. The Zeta function does not belong to the compilable subset.

f[x_] := ((Sin[x] Cos[x])/(1 + 
              Tan[x]^2)) (E^(-(x^2))) Zeta[x]

A plot of the function

Plot[f[x], {x, 2, 4}, PlotRange -> All];

To generate C++ code we have to declare the function arguments and return values.

Declare[f[Real x_] -> Real]

Declare[Zeta[Real x_] -> Real]

These declare statements does not change the way Mathematica computes the function

{f[2.5], f[(5/2)]}

Generate C++ code and compile it to an executable file. The option CallBackFunctions tells MathCode C++ which functions that has to be evaluated by Mathematica.

BuildCode[CallBackFunctions -> {Zeta}]

The generated executable is connected to Mathematica

InstallCode[];

Now the external code is used to compute the function

{f[2.5], f[(5/2)]}

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

A plot using the external code

Plot[f[x], {x, 2, 4}, PlotRange -> All];

Uninstall the executable, remove the temporary files, and remove the function f

UninstallCode[];
CleanMathCodeFiles[Confirm -> False, 
      CleanAllBut -> {}]
Remove[f]

Other Features

In this section we briefly describe some more features of MathCode.

Common subexpession elimination

Consider the following function whose definition contains a number of common subexpressions (e.g. 1 + (x^2) and (Sqrt[1 + (x^2)])).

g[Real x_] -> Real:=((x Cos[x] Cos[(Sqrt[1 + (x^2)])])/((((1 + (x^2)))^(3/2)) (1 + Cos[x]^2)))-((2 x Cos[x] Sin[(Sqrt[1 + (x^2)])])/((((1 + (x^2)))^2) (1 + Cos[x]^2)))+((2 Cos[x]^2 Sin[x] Sin[(Sqrt[1 + (x^2)])])/((1 + (x^2)) (((1 + Cos[x]^2))^2)))-((Sin[x] Sin[(Sqrt[1 + (x^2)])])/((1 + (x^2)) (1 + Cos[x]^2)))

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

We plot the function

pl1 = Plot[g[x], {x, -5, 5}];

We can tell MathCode C++ to do Common Subexpression Elimination (CSE) by the option EvaluateFunctions to CompileCode[ ] or BuildCode[ ].

BuildCode[EvaluateFunctions -> {g}]

The generated executable is connected to Mathematica

InstallCode[];

and we plot the externally evaluated function

pl2 = Plot[g[x], {x, -5, 5}];

We take a look at the generated C++ file

Run["Notepad Global.cc"];

Note how the computation of the function has been divided into small subexpressions that are only evaluated once and then stored in temporary variables for future use. This gives very efficient code for large functions!

Uninstall the executable, remove the temporary files, and remove the function g

UninstallCode[];
CleanMathCodeFiles[Confirm -> False, 
        CleanAllBut -> {}];
Remove[g];

Stand-Alone Code

MathCode C++ can be used to generate C++ programs that can be executed independently of Mathematica, that is to produce stand-alone programs. We define a simple function that prints something as result (hence, return type Null)

f[Real x_] -> Null := Print[(x^2)]

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.

CompilePackage[MainFileAndFunction -> "
	void main()
	{
		Global_Tf(1.2);
	}
"]

The stand-alone executable is generated

MakeBinary[StandAloneExecutable -> True]

We run the program "Global" and put the output in the file named result.txt

Run["Global>result.txt"];

The content of result.txt

Run["Notepad result.txt"]

Uninstall the executable, remove the temporary files, and remove the function f

UninstallCode[];
CleanMathCodeFiles[Confirm -> False, 
        CleanAllBut -> {}];
Remove[f];

Extended matrix operations

When dealing with matrices it is very convenient to have short notation for part extraction. The functionality of Part[] or [[ ]] has been extended to obtain this.

Consider the following 5×5 matrix:

A = Table[a[i, j], {i, 4}, {j, 5}]; A // MatrixForm

We extract row 2 to 4

A[[2 | 4]]  //  
      MatrixForm

We extract the elements in all rows which also belongs to column 3 and higher

A[[_, 3 | _]] // 
      MatrixForm

We can assign values to a submatrix of A.

A[[2 | 3, 
          2 | 3]]  =  {{1, 2}, {3, 4}};

A  //  MatrixForm

All these operations belongs to the compilable subset.

Code Generation from a Package

In this section we generate code for a package called MySurface. The example is similar to what is shown in the SinSurface.nb demo notebook.

The Package

Beginning of Package

BeginPackage["MySurface`", {MathCodeContexts}];
Off[General::"spell1"]
Clear["MySurface`*"];

Exported Symbols

Begin["MySurface`"];
\(	\(xyMatrix;\)\)
\(	\(sin;\)\)
\(	\(cos;\)\)
\(	\(arcTan;\)\)
\(	\(sinxy;\)\)
\(	\(computeGridValues;\)\)
End[];

Compilation Options

\(\(SetCompilationOptions[
	
        EvaluateFunctions -> {sin, cos, arcTan}, 
	
        CallBackFunctions -> {FresnelS}, 
	
        StandAloneExecutable -> True, 
	
        MainFileAndFunction -> "
		void main(){
		int i, j;
		i=-5;
		while(i<=5){
			i++; j=-5;
			while(j<=5){
				j++;
				printf("%f  ",MySurface_Tsinxy(0.4*i,0.4*j)); }
				printf("\n"); } }"
];\)\)

Private Context

Begin["MySurface`Private`"];

A global variable.

Declare[ Real[21, 21] xyMatrix ];

An function that have to be computed by a call-back to Mathematica.

Declare[FresnelS[Real x_] -> Real]

A number of functions with right-hand-sides (rhs) not in the compilable subset. However, the rhs can be evaluated by Mathematica to expressions that belong to the compilable subset.

sin[Real x_ ] -> Real := 
      Normal[Series[
            Sin[z],  {z, 0, 
              10}]]  /.  z -> x 
cos[Real x_ ] -> Real := 
      1 - Integrate[sin[z], {z, 0, x}]
arcTan[Real x_, Real y_] -> Real :=  
      If[x < 0, Pi, 0] + 
        If[x == 0, Sign[y]*Pi/2,  ArcTan[y/x] ]

The main function

\(\(sinxy[Real\ x_, Real\ y_] -> Real := \ 
	
        Module[
		{\ Real\ {r, xx, yy}\ }, 
		r = 
            Sqrt[x^2 + y^2]; 
		xx = \ r*cos[arcTan[x, y]]; 
		yy = 
            r*sin[arcTan[x, y]]; 
		sin[xx + yy]
\ \ ];\)\)

A function that computes function values of sinxy[ ] over a grid. A scaling using the FresnelS[ ] function is also performed.

\(computeGridValues[Integer\ n_]\  -> \ 
        Real[21, \ 
          21]\  := \ 
\ \ Module[{Integer\ {i, \ j, \ k}}, 
\ \ For[k = 1, 
          k ≤ n, 
          xyMatrix = 
            Table[sinxy[\(-2. \) + 0.2*i, \(-\ 2. \) + 0.2*j], {i, 1, 
                21}, {j, 1, 21}], \(k++\)]; 
		FresnelS[1. ]*
          xyMatrix]\)

End of Private Context

End[];

End of Package

EndPackage[];

Generate code for all functions in the package

Generate C++ code according to the compilation options given in the package

CompilePackage["MySurface`"];

The generated code

!! MySurface.cc

Use function expansion to avoid the auxiliary functions

Generate C++ code where the body of sinxy[ ] is symbolically expanded before code generation. This makes the auxiliary functions sin, cos, arcTan unnecessary and we choose not to compile them.

\(\(CompilePackage["MySurface`", 
	
        EvaluateFunctions -> {sinxy}, 
	
        UnCompiledFunctions -> {sin, cos, arcTan}];\)\)

The generated code

!! MySurface.cc

Generate Stand-Alone Code

Generate C++ code without code for FresnelS[ ] and computeGridValues[ ] since the main function does not use these (FresnalS[ ] cannot be evaluated without call-backs to Mathematica).

\(\(CompilePackage["MySurface`", 
	
        CallBackFunctions -> {}, 
	
        UnCompiledFunctions -> {FresnelS, 
            computeGridValues}];\)\)

MakeBinary[]

We run the external program and put the result in a file: res.txt

Run["MySurface.exe > res.txt"];

The contents of res.txt

!! res.txt

and a plot of the result

ListPlot3D[
        ReadList["res.txt", Number, 
          RecordLists -> True]];

The generated code

!! MySurface.cc

Macro rules

Instead of using call-backs to Mathematica to evaluate functions not in the compilable subset so called macro rules can be used to replace these functions with compilable expressions. We illustrate this by substituting the FresbelS[ ] function above with a polynomial.

Generate some data

data = Table[{x, FresnelS[x]}, {x, -2, 2,  .1}];

Define the function ownFresnelS[ ] whose right hand side is a polynomial fit of 10th degree to the data with "small" terms removed

ownFresnelS[x_] = Fit[data, (x^Range)[0, 10], x] // Chop

We check the result

Plot[{FresnelS[x], ownFresnelS[x]}, {x, -2.5, 2.5}];

Generate C++ code with

\(\(CompilePackage["MySurface`", 
	
        CallBackFunctions -> {}, 
	
        MacroRules -> {FresnelS[x_] -> 
              ownFresnelS[x]}];\)\)

The generated code

!! MySurface.cc

Clean-up

Uninstall the executable, remove the temporary files, and remove the functions in the MySurface` package.

CleanMathCodeFiles[Confirm -> False, 
        CleanAllBut -> {}];
Remove["MySurface`*"];

Delete the temporary directory and all files it contains.

SetDirectory[ToFileName[".."]];
DeleteDirectory["Tmp",  
        DeleteContents -> True];

Conclusions

Simple usage: the most basic usage of MathCode C++ can be summarized by three steps
	1)  Add type declarations
	2)  BuildCode[ ] generates C++ code and an executable program
	3)  InstallCode[ ] connects the executable program to Mathematica

MathCode C++ provides an easy way of calling external programs without knowledge about MathLink.

The combination of Mathematica and MathCode C++ is an excellent tool for prototyping of numerical algorithms.

Speeds up large numerical computations that can be performed by compiled code outside Mathematica.

Boost in performance for large numerical computations.

Transparent calling of compiled code from within Mathematica

Mathematica independent stand-alone code and programs can be generated

Keep track of the compilable subset of Mathematica
- Procedural programming constructs (For, While, Do, If, ...)
- Elementary functions (Sin, Exp, Log, ...)
- Arithmetic and logic expressions (+,-,*,/,||,<,!=,...) 
- Map, Apply, List, Table, Array, Module, Sum, Product, Dot, Max, Transpose, DiagonalMatrix, ...

How to use functions which does not belong to the compilable subset?
- Symbolic evaluation
- Call-backs to Mathematica
- Use external code
- Reimplement the function using functions in the compilable subset