| Using the C Math Library | |
Example - Passing Functions As Arguments (ex4.c)
To illustrate function-functions, this example program uses the ordinary differential equation (ODE) solver mlfOde23() to compute the trajectory of the Lorenz equation. Given a function, F, and a set of initial conditions expressing an ODE, mlfOde23() integrates the system of differential equations, y' = F(t,y), over a given time interval. mlfOde23() integrates a system of ordinary differential equations using second and third order Runge-Kutta formulas. In this example, the name of the function being integrated is lorenz.
For convenience, this example has been divided into three sections; in a working program, all of the sections would be placed in a single file. The first code section specifies header files, declares global variables including the local function table, and defines the lorenz function.
/* ex4.c */
#include <stdlib.h>
#include "matlab.h" /* 1 */
double SIGMA, RHO, BETA; /* 2 */
static mlfFuncTabEnt MFuncTab[] = /* 3 */
{
{"lorenz", (mlfFuncp)lorenz, _lorenz_thunk_fcn_ },
{ 0, 0, 0}
};
mxArray *lorenz(mxArray *tm, mxArray *ym) /* 4 */
{
mxArray *ypm = NULL;
double *y, *yp;
mlfEnterNewContext(0, 2, tm, ym);
mlfAssign(&ypm, mlfDoubleMatrix(3, 1, NULL, NULL)); /* 5 */
y = mxGetPr(ym);
yp = mxGetPr(ypm);
yp[0] = -BETA*y[0] + y[1]*y[2]; /* 6 */
yp[1] = -SIGMA*y[1] + SIGMA*y[2];
yp[2] = -y[0]*y[1] + RHO*y[1] - y[2];
mlfRestorePreviousContext(0, 2, tm, ym);
return mlfReturnValue(ypm);
}
The numbered items in the list below correspond to the numbered comments in the code example:
matlab.h". This file contains the declaration of the mxArray data structure and the prototypes for all the functions in the library. stdlib.h contains the definition of EXIT_SUCCESS.
SIGMA, RHO, and BETA, which are the parameters for the Lorenz equations. The main program sets their values, and the lorenz function uses them.
MFuncTab[], of type mlfFuncTabEnt. This variable stores a function table entry that identifies the function that mlfOde23() calls. A table entry contains three parts: a string that names the function ("lorenz"), a pointer to the function itself ((mlfFuncp)lorenz), and a pointer to the thunk function that actually calls lorenz, (_lorenz_thunk_fcn_). The table is terminated with a {0, 0, 0} entry.
Before you call mlfOde23() in the main program, pass MFuncTab to the
function mlfFevalTableSetup(), which adds your entry to the built-in
function table maintained by the MATLAB C Math Library. Note that a
table can contain more than one entry.
tm, containing the value of t, and a 3-by-1 array, ym, containing the values of y. The result is a new 3-by-1 array, ypm, containing the values of the three derivatives of the equation at time = t.
lorenz function.
lorenz doesn't use the input time step, tm, which is provided by mlfOde23.) Store the values directly in the real part of the array that lorenz returns. yp points to the real part of ypm, the return value.
The next section of this example defines the thunk function that actually calls lorenz. You must write a thunk function whenever you want to pass a function that you've defined to one of MATLAB's function-functions.
static int _lorenz_thunk_fcn_(mlfFuncp pFunc, int nlhs, /* 1 */
mxArray **lhs, int nrhs,
mxArray **rhs )
{
typedef mxArray *(*PFCN_1_2)( mxArray * , mxArray *); /* 2 */
mxArray *Out;
if (nlhs > 1 || nrhs > 2) /* 3 */
{
return(0);
}
Out = (*((PFCN_1_2)pFunc))( /* 4 */
(nrhs > 0 ? rhs[0] : NULL),
(nrhs > 1 ? rhs[1] : NULL)
);
if (nlhs > 0) /* 5 */
lhs[0] = Out;
return(1); /* 6 */
}
The numbered items in the list below correspond to the numbered comments in the code example:
lorenz function. A thunk function acts as a translator between your function's interface and the interface needed by the MATLAB C Math Library.
The thunk function takes five arguments that describe any function with
two inputs and one output (in this example the function is always lorenz()):
an mlfFuncp pointer that points to lorenz(), an integer (nlhs) that
indicates the number of output arguments required by lorenz(), an array
of mxArray's (lhs) that stores the results from lorenz(), an integer (nrhs)
that indicates the number of input arguments required by lorenz(), and an
array of mxArray's (rhs) that stores the input values. The lhs (left-hand side)
and rhs (right-hand side) notation refers to the output arguments that
appear on the left-hand side of a MATLAB function call and the input
arguments that appear on the right-hand side.
lorenz function pointer. The pointer to lorenz comes into the thunk function with the type mlfFuncp, a generalized type that applies to any function.
mlfFuncp is defined as follows:
typedef void (*mlfFuncp)(void)
The function pointer type that you define here must precisely specify the
return type and argument types required by lorenz. The program casts
pFunc to the type you specify here.
The name PFCN_1_2 makes it easy to identify that the function has 1 output
argument (the return) and 2 input arguments. Use a similar naming scheme
when you write other thunk functions that require different numbers of
arguments. For example, use PFCN_2_3 to identify a function that has two
output arguments and three input arguments.
lorenz expects two input arguments and one output argument. (The return value counts as one output argument.) Exit the thunk function if too many input or output arguments have been provided. Note that the thunk function relies on the called function to do more precise checking of arguments.
lorenz, casting pFunc, which points to the lorenz function, to the type PFCN_1_2. Verify that the two expected arguments are provided. If at least one argument is passed, pass the first element from the array of input values (rhs[0]) as the first input argument; otherwise pass NULL. If at least two arguments are provided, pass the second element from the array of input values (rhs[1]) as the second argument; otherwise pass NULL as the second argument. The return from lorenz is stored temporarily in the local variable Out.
This general calling sequence handles optional arguments. It is technically
unnecessary in this example because lorenz has no optional arguments.
However, it is an essential part of a general purpose thunk function.
Note that you must cast the pointer to lorenz to the function pointer type
that you defined within the thunk function.
lorenz to the appropriate position in the array of output values. The return value is always stored at the first position, lhs[0]. If there were additional output arguments, values would be returned in lhs[1], lhs[2], and so on.
The next section of this example contains the main program. Keep in mind that in a working program, all parts appear in the same file.
int main( )
{
mxArray *tm = NULL, *ym = NULL; /* 1 */
mxArray *tsm = NULL, *ysm = NULL;
double tspan[] = { 0.0, 10.0 };
double y0[] = { 10.0, 10.0, 10.0 };
double *t, *y1, *y2, *y3;
int k, n;
mlfEnterNewContext(0, 0);
mlfFevalTableSetup ( MFuncTab ); /* 2 */
SIGMA = 10.0; /* 3 */
RHO = 28.0;
BETA = 8.0/3.0;
mlfAssign(&tsm, mlfDoubleMatrix(2, 1, tspan, NULL));
mlfAssign(&ysm, mlfDoubleMatrix(1, 3, y0, NULL));
mlfAssign(&tm, mlfOde23(&ym, mlfVarargout(NULL), /* 4 */
mxCreateString("lorenz"), tsm, ysm, NULL, NULL));
n = mxGetM(tm); /* 5 */
t = mxGetPr(tm);
y1 = mxGetPr(ym);
y2 = y1 + n;
y3 = y2 + n;
mlfPrintf(" t y1 y2 y3\n"); /* 6 */
for (k = 0; k < n; k++) {
mlfPrintf("%9.3f %9.3f %9.3f %9.3f\n",
t[k], y1[k], y2[k], y3[k]);
}
mxDestroyArray(tsm); /* 7 */
mxDestroyArray(ysm);
mxDestroyArray(tm);
mxDestroyArray(ym);
mlfRestorePreviousContext(0,0);
return(EXIT_SUCCESS);
}
The numbered items in the list below correspond to the numbered comments in the code example:
tspan stores the start and end times. y0 is the initial value for the lorenz iteration and contains the vector 10.0, 10.0, 10.0.
feval function table by calling mlfFevalTableSetup(). The argument, MFuncTab, associates the string "lorenz" with a pointer to the lorenz function and a pointer to the lorenz thunk function. When mlfOde23() calls mlfFeval(), mlfFeval() accesses the library's built-in function table to locate the function pointers that are associated with a given function name, in this example, the string "lorenz".
SIGMA, RHO, and BETA. These parameters are shared between the main program and the lorenz function. The lorenz function uses the parameters in its computation of the values of the Lorenz equations.
tsm and ysm, which are passed as input arguments to the mlfOde23 function. Initialize tsm to the values stored in tspan. Initialize ysm to the values stored in y0.
mlfOde23(). The return value and the first argument store results. mlfOde23() is a varargout function; mlfVarargout(NULL) indicates that you are not interested in supplying any varargout arguments. Pass the name of the function, two required input arguments, and NULL values for the two optional input arguments.
mlfOde23() calls mlfFeval() to evaluate the lorenz function. mlfFeval()
searches the function table for a given function name. When it finds a match,
it composes a call to the thunk function that it finds in the table, passing the
thunk function the pointer to the function to be executed, also found in the
table. In addition, mlfFeval() passes the thunk function arrays of input and
output arguments. The thunk function actually executes the target function.
n time steps, the values in column 1 occupy positions 0 through n-1 in the result, the values in column 2, positions n through 2n-1, and so on.
tm returned from mlfOde23. The function mxGetM returned the number of rows in its mxArray argument.
Output
The output from this program is several pages long. Here are the last lines of the output.
t y1 y2 y3
9.390 41.218 12.984 2.951
9.405 39.828 11.318 0.498
9.418 38.530 9.995 -0.946
9.430 37.135 8.678 -2.043
9.442 35.717 7.404 -2.836
9.455 34.229 6.117 -3.409
9.469 32.711 4.852 -3.778
9.484 31.185 3.632 -3.972
9.500 29.657 2.477 -4.029
9.518 28.123 1.402 -3.989
9.539 26.563 0.415 -3.899
9.552 25.635 -0.116 -3.845
9.565 24.764 -0.576 -3.807
9.580 23.861 -1.014 -3.796
9.598 22.818 -1.478 -3.833
9.620 21.682 -1.948 -3.964
9.645 20.488 -2.429 -4.245
9.674 19.280 -2.960 -4.761
9.709 18.143 -3.618 -5.642
9.750 17.275 -4.545 -7.097
9.798 17.162 -6.000 -9.461
9.843 18.378 -7.762 -12.143
9.873 20.156 -9.147 -13.971
9.903 22.821 -10.611 -15.464
9.931 26.021 -11.902 -16.150
9.960 29.676 -12.943 -15.721
9.988 32.932 -13.430 -14.014
10.000 34.012 -13.439 -12.993
| | Extending the mlfFeval() Table | Replacing Argument Lists with a Cell Array | |