The first example shows how the package can be used to create a Fortran subroutine to compute a vector valued function.
|
>
|
fg := makeproc(v,parameters=[t,x,y]);
|
|
subroutine fg(t,x,y,crea_par)
doubleprecision t
doubleprecision x
doubleprecision y
doubleprecision crea_par(2)
doubleprecision t1
doubleprecision t5
doubleprecision v(2)
t1 = exp(-t)
v(1) = -y*x*t1+1
t5 = x**2
v(2) = -t5*t1*y+x
crea_par(1) = v(1)
crea_par(2) = v(2)
return
return
end
| |
The function f below computes the square of the distance from a point (x,y) to a circle of radius r with centre (h,k) under translation by (C,D) followed by rotation of theta radians. We use automatic differentiation to compute the gradient of the function f with respect to C,D,theta and generate C code for the result. Finally we output the computational cost of the optimized code.
|
>
|
f := proc(C,D,theta,x,y,h,k,r)
local s,c,xbar,ybar,d1,d2,d;
s := sin(theta);
c := cos(theta);
xbar := (x+C)*c + (y+D)*s;
ybar := (y+D)*c - (x+C)*s;
d1 := (h-xbar)^2;
d2 := (k-ybar)^2;
d := sqrt( d1+d2 ) - r;
d^2
end proc:
|
|
#include <math.h>
void G(C,D,theta,x,y,h,k,r,grd)
float C;
float D;
float theta;
float x;
float y;
float h;
float k;
float r;
float grd[3];
{
float c;
float d1;
float d2;
float dfr0[7];
float s;
float t1;
float t10;
float t12;
float t20;
float t22;
float t3;
float t7;
float t8;
float t9;
{
s = sin(theta);
c = cos(theta);
t1 = x+C;
t3 = y+D;
t7 = -c*t1-s*t3+h;
d1 = t7*t7;
t8 = -c*t3+s*t1+k;
d2 = t8*t8;
t9 = d1+d2;
t10 = sqrt(t9);
dfr0[6] = 2.0*t10-2.0*r;
t12 = sqrt(t9);
dfr0[5] = 1/t12*dfr0[6]/2.0;
dfr0[4] = dfr0[5];
dfr0[3] = -2.0*t8*dfr0[4];
dfr0[2] = -2.0*t7*dfr0[4];
t20 = dfr0[3];
t22 = dfr0[2];
dfr0[1] = t1*t22+t3*t20;
dfr0[0] = -t1*t20+t3*t22;
grd[0] = 0.0;
grd[1] = c*t20+s*t22;
grd[2] = c*dfr0[0]-s*dfr0[1];
return;
}
}
| |
In this example, we compute 1+x+x^2/2+x^3/6+...+x^n/n! We do this in two ways. Firstly, we use a symbolic sum. To obtain C code we explicitly convert the symbolic sum to a for loop.
|
>
|
F := makeproc(f,parameters=[n,x],locals=[i]);
|
|
double F(n,x)
int n;
double x;
{
double i;
int i1;
double s1;
double t1;
{
if( 0.0 < -n )
s1 = 0.0;
else
{
t1 = 1.0;
s1 = 1.0;
for(i1 = 1;i1 <= n;i1++)
{
t1 = x/i1*t1;
s1 += t1;
}
}
return(s1);
}
}
| |
In the second approach, our program uses a loop to compute the finite sum. We use the intermediate representation for programs provided by the codegen package. This representation is an expression tree. It can be converted into a Maple procedure using the intrep2maple command which can then be converted into Fortran or C code if desired.
In this example, we compute the gradient and hessian of a simple function, manipulate the procedures then join the two together.
|
>
|
F := makeproc(f,[x,y,t]);
|
This final example shows the functionality of the maple2intrep and intrep2maple commands. The maple2intrep command converts a Maple procedure into an intermediate representation which is suitable for manipulation. One must be careful about how one evaluates this representation because the functions in the code will evaluate.
|
>
|
f := proc(x,n) local A,i;
A := Array(0..n);
A[0] := 1;
for i to n do A[i] := x*A[i-1] end do;
A
end proc:
|
