Download package here

ToExternal
A package for high-speed numerical computation

 Introduction

The ToExternal package provides the capability of converting numerical Maple procedures into C code, compiling them, and linking them into the currently running kernel. This can be seen as a faster evalhf, as it is primarily designed for numerical computations.

In addition it can automatically be used by dsolve/numeric  for IVP and BVP problems, and pdsolve/numeric  for 1-time, 1-space dimension problems.

 Requirements & Limitations

The ToExternal package is a tool that ties together a number of things, but has strict requirements for use:

At this time, the intersection between the classes of functions that ToExternal can handle, and those that evalhf can handle is quite large, but there are the following major differences:

These restrictions are largely a limitation resulting from the primary application area of interest, namely fast solution of ODE and PDE problems.

 Setup for Use & Platform Issues

On all platforms, the library and help files should be unpacked in a single directory, and that directory should be added to Maple's libname on startup.

The easiest way to do this is by adding the line:

        libname := "installed directory",libname:

to your ~/.mapleinit  file (on UNIX) or   your maple.ini  file (on Windows - if not present this needs to be put in your Users directory)

On UNIX, you are now done, but for Windows there are additional steps needed to prepare it for use (due to operating system limitations).

1) You need to tell ToExternal the path to your mingw bin directory. If you installed mingw in "C:\Program Files\mingw" this can be accomplished (in Maple) with the command:

         ToExternal:-SetOption(mingw,compiler_path,"C:\\Program Files\\mingw\\bin");  

2) You need to set a temp directory. If, for example, you want to use "C:\temp" this can be accomplished with:

         ToExternal:-SetOption(mingw,temp_directory,"C:\\Temp");  

And now you should be ready to run.

Important : Why the extra requirements on Windows?

The second requirement is by far the ugliest. The problem is that under Windows it is not possible to delete a shared library while it is in use (Windows locks it), and unfortunately it is in use by Maple until after Maple is closed (so code running in Maple cannot remove it).

What this means is that any shared libraries created under Windows will accumulate, and never be deleted.

In addition, the Maple system directory is the default current directory when launching Maple, so it would be used for the shared libraries.

Choosing a temporary directory to place these in is useful, as then they can be cleaned out easily.

Once no Maple sessions are running that use the shared libraries, simply delete them all via " del m2c*.dll " from the temporary directory.

The first requirement is due to DOS shell limitations and PATH limitations (as the compile is a command line process).

Interface-specific issues:

When running ToExternal calls in Windows in the GUI, a dos shell will pop up whenever a system command (read as "call to the compiler") is run. This is not very pretty, but it is another O/S limitation.

 Functions

There are currently only 3 top-level functions available in the package, and these are:

 Automatic use by dsolve/numeric and pdsolve/numeric

Use of ToExternal by dsolve/numeric and pdsolve/numeric is accomplished through environment variable settings.

Maple's default IVP and DAE solvers can be set to use ToExternal by setting:

>    _Env_ivpnumeric_ToExternal := true:

Maple's default BVP solver can be set to use ToExternal by setting:

>    _Env_bvpnumeric_ToExternal := true:

And Maple's default PDE solver can be set to use ToExternal by setting:

>    _Env_pdenumeric_ToExternal := true:

In all cases, problems which cannot use ToExternal (for example, and ODE system involving Bessel functions) will be detected, and the default method will be used instead.

Detection of the success of the use of ToExternal can be accomplished by setting the infolevel of the following to 1:

>    infolevel[ToExternal:-HotLink] := 0:

>   

 Examples

IVP Example, essentially a pair of second order ODE, solved for 72 sets of initial conditions simultaneously.

>    N := 72:
Sodes :=[seq(op([ diff(x[i](t),t,t)=-x[i](t)^2+y[i](t)^2,
                  diff(y[i](t),t,t)=2*x[i](t)*y[i](t)]),
             i=1..N)]:
Sics := [seq(op([ x[i](0)=cos(2*Pi*(i-1)/N), D(x[i])(0)=0,
                  y[i](0)=sin(2*Pi*(i-1)/N), D(y[i])(0)=0]),
             i=1..N)]:

First with ToExternal off:

>    _Env_ivpnumeric_ToExternal := false:
ds := dsolve({op(Sodes),op(Sics)},'numeric',abserr=1e-12,relerr=1e-12):
tt := time():
plots[odeplot](ds, [seq([x[i](t),y[i](t)],i=1..N)], 0..2.9,
               colour=black, view=[-3..3,-3..3],
               scaling=CONSTRAINED, labels=["",""]);
time()-tt;

[Maple Plot]

.971

Now with it on:

>    _Env_ivpnumeric_ToExternal := true:
ds := dsolve({op(Sodes),op(Sics)},'numeric',abserr=1e-12,relerr=1e-12):
tt := time():
plots[odeplot](ds, [seq([x[i](t),y[i](t)],i=1..N)], 0..2.9,
               colour=black, view=[-3..3,-3..3],
               scaling=CONSTRAINED, labels=["",""]);
time()-tt;

[Maple Plot]

.210

Not a massive difference, but noticable.

 

Stiff DAE Example, 4 weight pendulum in cartesian co-ordinates, until length strings attaching each pair of unit weights:

>    g := 98/10:
dsys := {
   diff(x[1](t),t,t)=-l[1](t)*(x[1](t)-0)      +l[2](t)*(x[2](t)-x[1](t)),
   diff(x[2](t),t,t)=-l[2](t)*(x[2](t)-x[1](t))+l[3](t)*(x[3](t)-x[2](t)),
   diff(x[3](t),t,t)=-l[3](t)*(x[3](t)-x[2](t))+l[4](t)*(x[4](t)-x[3](t)),
   diff(x[4](t),t,t)=-l[4](t)*(x[4](t)-x[3](t))+0,
   diff(y[1](t),t,t)=-l[1](t)*(y[1](t)-0)      +l[2](t)*(y[2](t)-y[1](t))-g,
   diff(y[2](t),t,t)=-l[2](t)*(y[2](t)-y[1](t))+l[3](t)*(y[3](t)-y[2](t))-g,
   diff(y[3](t),t,t)=-l[3](t)*(y[3](t)-y[2](t))+l[4](t)*(y[4](t)-y[3](t))-g,
   diff(y[4](t),t,t)=-l[4](t)*(y[4](t)-y[3](t))+0                        -g,
   (x[1](t)-0)^2      +(y[1](t)-0)^2      =1,
   (x[2](t)-x[1](t))^2+(y[2](t)-y[1](t))^2=1,
   (x[3](t)-x[2](t))^2+(y[3](t)-y[2](t))^2=1,
   (x[4](t)-x[3](t))^2+(y[4](t)-y[3](t))^2=1}:
ini := {x[1](0)=1,D(x[1])(0)=0,y[1](0)= 0,D(y[1])(0)=0,
        x[2](0)=1,D(x[2])(0)=0,y[2](0)=-1,D(y[2])(0)=0,
        x[3](0)=1,D(x[3])(0)=0,y[3](0)=-2,D(y[3])(0)=0,
        x[4](0)=1,D(x[4])(0)=0,y[4](0)=-3,D(y[4])(0)=0}:
        

>    _Env_ivpnumeric_ToExternal := false:
ds := dsolve(dsys union ini,numeric,stiff=true,implicit=true):
tt := time():
plots[odeplot](ds,[x[4](t),y[4](t)],0..10,numpoints=1000);
time()-tt;

[Maple Plot]

13.821

>    _Env_ivpnumeric_ToExternal := true:
ds := dsolve(dsys union ini,numeric,stiff=true,implicit=true):
tt := time():
plots[odeplot](ds,[x[4](t),y[4](t)],0..10,numpoints=1000);
time()-tt;

[Maple Plot]

.200

This difference is somewhat more noticable.

Jeffery-Hamel Flow, 4rth order BVP Problem with continuation and large Reynolds number (Re=1000)

>    dsys := {diff(Q(x),x$4)+2*Re*diff(Q(x),x)*diff(Q(x),x,x)+4*diff(Q(x),x,x),
         Q(0)=0, Q(Pi/2)=-2/3, (D@@2)(Q)(0)=0, D(Q)(Pi/2)=0};
csys := subs(Re=500+500*lambda,dsys):

dsys := {Q(0) = 0, Q(1/2*Pi) = -2/3, `@@`(D,2)(Q)(0) = 0, D(Q)(1/2*Pi) = 0, diff(Q(x),`$`(x,4))+2*Re*diff(Q(x),x)*diff(Q(x),`$`(x,2))+4*diff(Q(x),`$`(x,2))}

>    _Env_bvpnumeric_ToExternal := false:
tt := time():
ds := dsolve(csys,numeric,continuation=lambda,abserr=1e-10,maxmesh=256);
time()-tt;

ds := proc (x_bvp) local res, data, solnproc, ndsol, outpoint, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; _EnvDSNumericSaveDigits := Digits; Digits := 15; if _EnvInFsol...

.599

>    _Env_bvpnumeric_ToExternal := true:
tt := time():
ds := dsolve(csys,numeric,continuation=lambda,abserr=1e-10,maxmesh=256);
time()-tt;

ds := proc (x_bvp) local res, data, solnproc, ndsol, outpoint, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; _EnvDSNumericSaveDigits := Digits; Digits := 15; if _EnvInFsol...

.461

For this problem, the timing includes the compile time. The difference is less significant here because the function evaluations are a smaller expense in the solution process.

The following is a fairly difficult Hyperbolic PDE problem, originally given as a complex valued problem, but converted to a real one by splitting the real and imaginary parts of the solution.

>    theta[B] := 185456/10^6:
lambda := 708/10^13:
X[h] := -1901/10^9+I*1589/10^11:
Xmax := 180/10^6:
Zmin := 0:
Zmax := 800/10^6:
pde1 := -sin(theta[B])*diff(d0(x,z),x) + cos(theta[B])*diff(d0(x,z),z)-I/2*2*Pi/lambda*X[h]*dh(x,z):
pde2 := sin(theta[B])*diff(dh(x,z),x) + cos(theta[B])*diff(dh(x,z),z)-I/2*2*Pi/lambda*X[h]*d0(x,z):
fd0 := (x,z)->d0r(x,z)+I*d0i(x,z):
fdh := (x,z)->dhr(x,z)+I*dhi(x,z):
psys1 := eval(expand(eval({pde1,pde2},{d0=fd0,dh=fdh})),I=II):
psys2 := map(coeff,psys1,II,0) union map(coeff,psys1,II,1);
f(x) := piecewise(x < -5/10^6, 0,
                  -5/10^6 <= x and x <= 5/10^6,
                      10*(x^2-25/10^12)^8/((5.0/10^6)^16),
                   x > 5/10^6, 0):
IBC := {d0r(x,0) = f(x),  d0i(x,0)=0,
        d0r(Xmax,z) = 0,  d0i(Xmax,z) = 0,
        dhr(x,0) = 0,     dhi(x,0) = 0,
        dhr(-Xmax,z) = 0, dhi(-Xmax,z) = 0};

psys2 := {-sin(11591/62500)*diff(d0r(x,z),x)+cos(11591/62500)*diff(d0r(x,z),z)+39725/177*Pi*dhr(x,z)-4752500/177*Pi*dhi(x,z), sin(11591/62500)*diff(dhr(x,z),x)+cos(11591/62500)*diff(dhr(x,z),z)+39725/1...
psys2 := {-sin(11591/62500)*diff(d0r(x,z),x)+cos(11591/62500)*diff(d0r(x,z),z)+39725/177*Pi*dhr(x,z)-4752500/177*Pi*dhi(x,z), sin(11591/62500)*diff(dhr(x,z),x)+cos(11591/62500)*diff(dhr(x,z),z)+39725/1...
psys2 := {-sin(11591/62500)*diff(d0r(x,z),x)+cos(11591/62500)*diff(d0r(x,z),z)+39725/177*Pi*dhr(x,z)-4752500/177*Pi*dhi(x,z), sin(11591/62500)*diff(dhr(x,z),x)+cos(11591/62500)*diff(dhr(x,z),z)+39725/1...
psys2 := {-sin(11591/62500)*diff(d0r(x,z),x)+cos(11591/62500)*diff(d0r(x,z),z)+39725/177*Pi*dhr(x,z)-4752500/177*Pi*dhi(x,z), sin(11591/62500)*diff(dhr(x,z),x)+cos(11591/62500)*diff(dhr(x,z),z)+39725/1...

IBC := {d0r(x,0) = PIECEWISE([0, x < -1/200000],[.6553600001e86*(x^2-1/40000000000)^8, -1/200000 <= x and x <= 1/200000],[0, 1/200000 < x]), d0r(9/50000,z) = 0, d0i(9/50000,z) = 0, dhr(-9/50000,z) = 0,...
IBC := {d0r(x,0) = PIECEWISE([0, x < -1/200000],[.6553600001e86*(x^2-1/40000000000)^8, -1/200000 <= x and x <= 1/200000],[0, 1/200000 < x]), d0r(9/50000,z) = 0, d0i(9/50000,z) = 0, dhr(-9/50000,z) = 0,...

First without ToExternal:

>    _Env_pdenumeric_ToExternal := false:
sol := pdsolve(psys2, IBC, numeric, spacestep=Xmax/500, timestep=Xmax/100):
tt:=time():
sol:-plot([[d0r,color=red],[d0i,color=blue],
           [dhr,color=green],[dhi,color=cyan]],
            z=Xmax,view=[-5e-5..5e-5,DEFAULT]);
time()-tt;

[Maple Plot]

2.101

And with:

>    _Env_pdenumeric_ToExternal := true:
sol := pdsolve(psys2, IBC, numeric, spacestep=Xmax/500, timestep=Xmax/100):
tt:=time():
sol:-plot([[d0r,color=red],[d0i,color=blue],
           [dhr,color=green],[dhi,color=cyan]],
            z=Xmax,view=[-5e-5..5e-5,DEFAULT]);
time()-tt;

[Maple Plot]

.329

>