MathCode

an Mathematica Application Package for Code Generation

by Peter Fritzson and Johan Gunnarsson, Linköping University
petfr@ida.liu.se, johan@isy.liu.se

Background

Why Code Generation?

Why using Mathematica?

MathCode Features

Generates a subset of Mathematica to C++ including:

Polynomials

Elementary mathematical functions

Numeric Arrays and connection to linear algebra operation library

Procedural styled programs.

Adds to Mathematica:

Convenient typsetting for variables and functions

Global variables

Local variables

Function signatures

Connection to external software written in C++ and Fortran.

Special keywords giving computations of a mixture of generated C++ code and symbolic expression in Mathematica.

Convenient Matlab/Fortran90 like Part[] for matrices.

The columns 3-5 from A

MathCode Developing Environment

Automatic generation of code, compilation and integration into Mathematica using Mathlink.

On the fly interchanging external code with the internal Mathematica code in computations.

Increasing the compatible subset of MathCode

Simple Example

Load the MathCode system into Mathematica

Needs["MathCode`"]

Set Working Directory

SetDirectory["C:/temp"];

Define a simple function

Typed definition

(Integer -> Real) foo[n_] := Sin[1/n]

Evaluation in Mathematica

foo[2]

foo is still a normal Mathematica function

??foo

the type information is stored somewhere else and is only used for the code generation

getDefinition[foo]

A Global Matrix

Define a matrix

Var[
	Real[3,3] xymatrix 
	]

This matrix definition will also create/initialize the matrix. This to make, e.g, xymatrix[[2,3]] valid directly.

xymatrix

Matrix Function

Initializer

(Real -> {}) xyinit[x_] := xymatrix = Table[x,{3},{3}];

Updates the matrix

(Real-> Real[_,_]) fee[x_] := Module[ { Real[3,3] m },
	    m = Table[foo[i+j],{i,3},{j,3}]; 
		xymatrix = xymatrix + m;
		xymatrix
]

xyinit[0. ]

fee[1]

fee[1]

fee[1]

Generate C++ and Build Executables

BuildCode["Global"]

Install the Mathlink Executable

InstallCode["Global"]

Run some examples

xyinit[1.]

fee[2.]

fee[2.]

Switch between Mathematica and Externals

The different definitions (DownValues) are stored.

?? foo

To use the Mathematica source definition instead

DeactivateCode[foo]

?? foo

Uninstall the code

UninstallCode["Global"]

which will clean up things and terminate the MathLink process.

?? foo

A glance at the Generated Code

!! Global.cc

This code is used to produce two executables.
One executeable for the MathLink connection

FileNames["*ml.exe"]

and another for standalone computations

FileNames["Global.exe"]

The main function for the standalone executable can be given to BuildCode as an option.

BuildCode["Global",
	MainFileAndFunction ->
		"int main() {
			Global_Txyinit(1.0);
			cout << Global_Tfee(1.0);
			return (0);
		}"]

Run the stand alone executable

! "Global > tmp"

!! tmp

Subexpression Elimination

For

A Large (Realistic) Example

Compute a rigidity matrix and second member for a finite element model a geometrically nonlinear 3D beam element. (Example from paubert@mathinsa.insa)

The Energy Expression

\(Needs["LinearAlgebra`MatrixManipulation`"]; 

    
 (*\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ *) 
 (*\ 
      shape\ function\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ *) 
 (*
      \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ *) 
Ni\  = \ 1 - x/l; 
Nj\  = x/l; 

    
 (*\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ *) 
 (*\ interpolations\ and\ derivatives\ \ \ \ \ \ \ \ \ *) 
 (*
      \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ *) 
u1x\  = \ Ni*u1i\  + \ Nj*u1j; 
u2x\  = \ Ni*u2i\  + \ Nj*u2j; 
    
u3x\  = \ Ni*u3i\  + \ Nj*u3j; 
t1x\  = \ Ni*t1i\  + \ Nj*t1j; 

    t2x\  = \ Ni*t2i\  + \ Nj*t2j; 
t3x\  = \ Ni*t3i\  + \ Nj*t3j; 


    du1 = D[u1x, x]; 
du2 = D[u2x, x]; 
du3 = D[u3x, x]; 
dt1 = D[t1x, x]; 
    
dt2 = D[t2x, x]; 
dt3 = D[t3x, x]; 

    
 (*\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ *) 
 (*\ compute\ a\ closed\ form\ of\ the\ expential\ matrix\ *) 
 (*
      \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ *) 
al = Sqrt[t1x^2 + t2x^2 + t3x^2]; 
a\ \  = Sin[al]/al; 
    
a1 = 1 - a; 
a2 = 1/2*\((Sin[al/2]/\((al/2)\))\)^2	; 
x1 = 1/2*l; 

    theta = \ {{0, \(-t3x\), t2x}, {t3x, 0, \(-t1x\)}, {\(-t2x\), t1x, 0}}; 

    vec\ \ \ \  = \ {t1x, t2x, t3x}; 
uni\  = \ {1, 0, 0}; 

    du = {du1, du2, du3}; 
dt = {dt1, dt2, dt3}; 

    
 (*\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ *) 
      
 (*\ suite\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ *) 
 (*
      \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ *) 

    rexp = IdentityMatrix[3] + a\ Transpose[theta] + a2\ theta . theta; 

    tenseurmat = \((rexp - IdentityMatrix[3])\) . uni + rexp . du; 

    texp = a\ IdentityMatrix[3] + Outer[Times, vec, vec] + 
        a2\ Transpose[theta]; 
tenseurrot\  = texp . dt; 

    
 (*\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ *) 
 (*\ 
      Elasticity\ tensor\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ *) 
 (*
      \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ *) 
ElastMat = DiagonalMatrix[{EA, GAy, GAz}]; 

    ElastRot = DiagonalMatrix[{EIx, EIy, EIz}]; 


    Energie = 
      1/2\ tenseurmat . ElastMat . tenseurmat\  + 
        1/2\ \ tenseurrot . \ ElastRot . tenseurrot; \)

Energie  //  Short

The second memeber element

\(\( (*
      \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ *) 
 (*\ 
      vector\ of\ ddl\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ *) 
 (*
      \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ *) 
\(ddl = {u1i, u2i, u3i, t1i, t2i, t3i, u1j, u2j, u3j, t1j, t2j, 
        t3j}; \)
\)\)

\(\( (*
      \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ *) 
 (*\ second\ member\ F\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ *) 
      
 (*\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ *) 
\(F = Table[D[Energie, ddl[\([i]\)]], {i, 1, Length[ddl]}]; 
    \)\)\)

To print this is not feasible. Compute the size.

ByteCount[F]

Generate Code

Define Ffunction as the first 4 elements of F

({Real,Real,Real,Real,Real,Real,Real,Real,Real,Real,
Real,Real,Real,Real,Real,Real,Real,Real,Real,Real} -> Real[_]
) Ffunction[EA_,EIx_,EIy_,EIz_,GAy_,GAz_,l_,t1i_,t1j_,t2i_,t2j_,
	   t3i_,t3j_,u1i_,u1j_,u2i_,u2j_,u3i_,u3j_,x_] = Take[F,4]

To be able to messarue the perfomance I need function looping over Ffunction

\(\(\((Integer\  -> \ Real[_])\)\ FfunLoop[n_]\  := \ 
      Module[{Real[_]\ \ res}, 
		Print[n]; 
		
        Do[res\  = 
            Ffunction[0.65, \ 0.65, \ 0.65, \ 0.65, \ 0.65, \ 0.65, \ 0.65, \ 
              0.65, \ 0.65, \ 0.65, \ 0.65, \ 0.65, \ 0.65, \ 0.65, \ 0.65, \ 
              0.65, \ 0.65, \ 0.65\ , 0.65, \ 0.65], {n}]; 
		res
		]; 
    \)\)

Cell["Generate, compile and link", "Text"

BuildCode["Global", EvaluateFunctions -> {Ffunction}]

Take Look at the Code

Do this with an editor in this case

Compare Performance

Restore some memory

Share[]

Normal Mathematica Evaluation

Ffunction[Sequence@@ Table[0.65, {20}]]  //  Timing

mathtime  =  %[[1]];

External Computation

Load the externaly compiled code

InstallCode["Global"]

Create a new timing function

\(SetAttributes[AbsTime, HoldFirst]; 


    AbsTime[x_]\  := \ 
      Module[{
\ \ start, res
\ }, 
\ \ \ start\  = \ AbsoluteTime[]; 

        \ \ \ res = x; 

        \ \ \ {\((AbsoluteTime[] - start)\)\ Second, res}
]; \)

To remove the influence the MathLink duration we compute the

AbsTime[FfunLoop[200000]]

exttime  =  %[[1]] / 200000

How much faster?

mathtime / exttime

Mathematica Compile

Switch to the non external implementation

DeactivateCode[Ffunction,  FfunLoop]

Use Compile on the Ffunction

\(\(fc\  = \ 
      Compile[{{EA, \ _Real}, 
			{EIx, _Real}, \ {EIy, \ _Real}, 
          \ {EIz, \ _Real}, \ {GAy, \ _Real}, \ {GAz, \ _Real}, 
          \ {l, \ _Real}, \ {t1i, \ _Real}, \ {t1j, \ _Real}, 
          \ {t2i, \ _Real}, 
          \ {t2j, \ _Real}, {\ t3i, \ _Real}, {\ t3j, \ _Real}, 
          \ {u1i, \ _Real}, {\ u1j, \ _Real}, \ {u2i, \ _Real}, 
          \ {u2j, \ _Real}, {\ u3i, \ _Real}, \ {u3j, \ _Real}, 
          \ {x, \ _Real}}, \ 
        Evaluate[
          Ffunction[EA, \ EIx, \ EIy, \ EIz, \ GAy, \ GAz, \ l, \ t1i, \ t1j, 
            \ t2i, \ t2j, \ t3i, \ t3j, \ u1i, \ u1j, \ u2i, \ u2j, \ u3i, \ 
            u3j, \ x]]]; \)\)

fc[Sequence@@ Table[0.65, {20}]]  //  Timing

comptime  =  %[[1]];

mathtime / comptimecomptime / exttime

UninstallCode["Global"];

Multiple Package Code Generation

Create Two Packages

BeginPackage["Foo1`", {MathCodeContexts}]\(	fee1\)Begin["`Private`"]\(	\((Real\  -> \ Real)\)\ fee1[x_] := \ Sin[1/x]\)End[]EndPackage[]

BeginPackage["Foo2`", {"Foo1`", MathCodeContexts}]\(	fee2; 
	fee3; 
Begin["`Private`"]\)\(	\((Real\  -> \ Real)\)\ fee2[x_] := \ fee3[1/x] + fee1[x]\)\(	\((Real\  -> \ Real)\)\ fee3[x_]\  := \ Cos[x]\)End[]EndPackage[]

Code Generation

Try to generate code for Foo2

CompilePackage["Foo2"]

The file Foo1.mh is called the Math Header file, which stores necessary information of a compiled package such that code generation of packages depending on this package can be performed.

CompilePackage["Foo1"]

!! Foo1.mh

Try now to compiled the Foo2 package

CompilePackage["Foo2"]

!! Foo2.mh

The generated code

!! Foo2.cc

MakeBinary["Foo2"]

MathCode Utilities

Not supported Mathematica function. What to do?

Connecting to External Code

The Inverse example

BeginPackage["MyNumeric`",{MathCodeContexts}]

myinv;

Begin["`Private`"]

(Real[_, _] -> Real[_, _]) myinv[mat_] := Inverse[mat + 1];

End[];
EndPackage[];

BuildCode["MyNumeric"]

!! MyNumeric.mh

InstallCode["MyNumeric"]

m = Table[Random[], {10}, {10}];

m1  =  Inverse[m];

m2  =  myinv[m - 1];

(m1 - m2)[[1]]

CallBack function

Some functions should be left to Mathematica to compute