...it will do this by asking for a pointer to a function. You will usually
be writing this function. You don't need to write it in Fortran; you
may write it in C.
However, be aware of the Fortran conventions; for example, in Fortran
a function taking double and returning double is
actually a function taking a double pointer and returning
a double. So, if you wanted to pass the sine function to your Slatec
function you would need to define a function like
Consider the Slatec zero finding function dfzero which finds
zeros of function F taking doubles and returning double:
SUBROUTINE DFZERO (F, B, C, R, RE, AE, IFLAG)
The function F above is said to be, reference the Slatec
source code documentation:
C F :EXT - Name of the DOUBLE PRECISION external function. This
C name must be in an EXTERNAL statement in the calling
C program. F must be a function of one DOUBLE
C PRECISION argument.
You could write your functions in Fortran. For example, the function below is acceptable for dfzero:
DOUBLE PRECISION FUNCTION F(X)
DOUBLE PRECISION X
F=X*X+2.0*X+3.0
RETURN
END
However, for the novice in Fortran, writing your functions in C is preferable
in order to avoid Fortran programming and to keep the compiling simpler.
To imitate the type of Fortran function above we declare a type of function
in C below:
typedef double ddfortran(double*);
Now, some examples of functions in C are:
ddfortran fortranSin;
ddfortran fortranCos;
double fortranSin(double* x) { /* sin() in fortran format */
double y;
y=(*x);
return(sin(y));
}
double fortranCos(double* x) { /* cos() in fortran format */
return(cos(*y));
}
Click here for a commentated example testing dfzero.
Click here for an extension of dfzero, written in C, that will search for all the zeros in an interval.
Another example of Slatec functions requiring the use of user specified functions is in the collection of ode-solvers, (DEPAC). Here is an example of the use of one these solvers, DDEABM (Adams Bashford Moulton) solver.
To describe a system of first order differential equations, we need a function,
C DF -- This is the name of a subroutine which you provide to C define the differential equations. C DF -- Provide a subroutine of the form C DF(X,U,UPRIME,RPAR,IPAR) C to define the system of first order differential equations C which is to be solved. For the given values of X and the C vector U(*)=(U(1),U(2),...,U(NEQ)) , the subroutine must C evaluate the NEQ components of the system of differential C equations DU/DX=DF(X,U) and store the derivatives in the C array UPRIME(*), that is, UPRIME(I) = * DU(I)/DX * for C equations I=1,...,NEQ. C C Subroutine DF must NOT alter X or U(*). You must declare C the name df in an external statement in your program that C calls DDEABM. You must dimension U and UPRIME in DF. C C RPAR and IPAR are DOUBLE PRECISION and INTEGER parameter C arrays which you can use for communication between your C calling program and subroutine DF. They are not used or C altered by DDEABM. If you do not need RPAR or IPAR, C ignore these parameters by treating them as dummy C arguments. If you do choose to use them, dimension them in C your calling program and in DF as arrays of appropriate C length.
The correct type of function, in C, for the above DF is:
typedef void (defunc)(double*,double*,double*,double*,int*);For the damped driven pendulum (one dim'l second order=two dim'l first order DE) we would define a function like this:
defunc df;
void df(double* t, double* X, double* V, double* dpar, int* ipar) {
double om, al, f,x,v;
om=dpar[0];
al=dpar[1];
f=dpar[2];
x=X[0];
v=X[1];
V[0]=X[1]; /* X[0] is position X[1] is velocity. */
V[1]=-W0*W0*sin(x)-al*v+cos(om*(*t));
} /* damped driven pendulum */