|
|
Evaluating a multi-dimensional integral (d01wcc) |
In this example we will demonstrate two ways of providing a user-defined function to a NAG routine in Scilab. The function we shall use is d01wcc, which evaluates a multidimensional integral of the function f(x) where x is a multidimensional vector. The user-supplied function will evaluate this function at any x.
The two ways of providing this function to the NAG routine are either to write the function in C, or to write it as a Scilab function. Both approaches have their advantages and disadvantages. The C code is a lot easier to write and requires fewer lines, but in order to change it, each time you would need to rebuild the interface. Writing the function in Scilab code is slightly harder, as you need to write a second interface for the Scilab function in C, but it does give you the flexibility to change your function without having to rebuild the Scilab-NAG interface.
According to the C Library Manual, the prototype for function d01wcc looks like this:
#include <nag.h>
#include <nagd02.h>
void nag_multid_quad_adapt_1(Integer ndim, double (*f)(Integer ndim, const double x[], Nag_User *comm),
const double a[], const double b[], Integer *minpts, Integer maxpts, double eps,
double *finval, double *acc, Nag_User *comm, NagError *fail)
Note that the user-supplied function is named f, and is
of type double.
For more information about the argument list of nag_multid_quad_adapt_1,
refer to the NAG C Library manual.
Here is the source code of our interface written in C nag_intext4a.c with the user-supplied function also written in C:
/* Example 4a
=========
Shows how to implement a Scilab wrapper to a NAG function using a User
defined function written in C */
#include "stack-c.h"
#undef Complex
#define Complex NagComplex
#include <nag.h>
#include <nag_stdlib.h>
#include <nagd01.h>
#define SciComplex doublecomplex
#include "stdio.h"
static double f(Integer n, double z[], Nag_User *comm);
int nag_intext4a(char *fname)
{
// to call this function in scilab use:
// [finval,acc,minpts,fail] = nag_multid_quad_funa(ndim,a,b,minpts,maxpts,eps)
int m1,n1,l1;
int m2,n2,l2;
int m3,n3,l3;
int m4,n4,l4;
int m5,n5,l5;
int m6,n6,l6;
int l7;
int l8;
int l9;
int l10;
int n, i, min, max;
Integer ndim, minpts, maxpts;
double eps, finval, acc;
double *a=0, *b=0;
Nag_User comm;
NagError ifail;
// define minimum and maximum left and right hand arguments
int minlhs=1, minrhs=6, maxlhs=4, maxrhs=6;
Nbvars = 0;
CheckRhs(minrhs, maxrhs);
CheckLhs(minlhs,maxlhs);
// Initialise ifail
INIT_FAIL(ifail);
GetRhsVar(1, "i", &m1, &n1, &l1); // input ndim
GetRhsVar(2, "d", &m2, &n2, &l2); // input a
GetRhsVar(3, "d", &m3, &n3, &l3); // input b
GetRhsVar(4, "i", &m4, &n4, &l4); // input minpts
GetRhsVar(5, "i", &m5, &n5, &l5); // input maxpts
GetRhsVar(6, "d", &m6, &n6, &l6); // input eps
//******************************************************************
// Read in input variables from Scilab and convert to NAG variables
//******************************************************************
// Read in ndim
//*============*
if (m1!=1 || n1!=1)
{
sciprint("%s: Dimension should be 1x1 character for arg 1\r\n",fname);
Error(999); goto END;
}
else
{
n = *istk(l1);
ndim = (Integer) n;
}
// Allocate memory for a(ndim) and b(ndim)
//========================================
if ( !(a = NAG_ALLOC(n,double) ) ||
!(b = NAG_ALLOC(n,double)))
{
sciprint("Allocation failure\n");
Error(999); goto END;
}
// Read in a
//===========
if (m2!= n || n2!=1)
{
sciprint("%s: Dimension should be ndim by 1 for arg 2\r\n",fname);
sciprint("n = %i", n);
Error(999); goto END;
}
else
{
for (i=0; i<n; i++)
{
a[i] = *stk(l2+i);
}
}
// Read in b
//==============
if (m3!=n || n3!=1)
{
sciprint("%s: Dimension should be ndim by 1 for arg 3\r\n",fname);
Error(999); goto END;
}
else
{
for(i=0; i<n; i++)
{
b[i] = *stk(l3+i);
}
}
// Read in minpts
//============
if (m4!=1 || n4!=1)
{
sciprint("%s: Dimension should be 1x1 for arg 4\r\n",fname);
Error(999); goto END;
}
else
{
min = *istk(l4);
minpts = (Integer)min;
}
// Read in maxpts
//============
if (m5!=1 || n5!=1)
{
sciprint("%s: Dimension should be 1x1 for arg 5\r\n",fname);
Error(999); goto END;
}
else
{
max = *istk(l5);
maxpts = (Integer) max;
}
// Read in eps
//==============
if (m6!=1 || n6!=1)
{
sciprint("%s: Dimension should be 1x1 for arg 6\r\n",fname);
Error(999); goto END;
}
else
{
eps = *stk(l6);
}
//***********************************
// End of reading in input variables
//***********************************
// nag_multid_quad_adapt_1 (d01wcc).
//* Multi-dimensional adaptive quadrature, thread-safe
//*====================================================
nag_multid_quad_adapt_1(ndim, f, a, b, &minpts, maxpts, eps, &finval, &acc,
&comm, &ifail);
// Print NAG error message if routine fails
//=========================================
if (ifail.code != NE_NOERROR)
sciprint("Ifail message is %s", ifail.message);
// Create output variables for Scilab
//====================================
CreateVar(7, "d", &n1, &m1, &l7); // finval
CreateVar(8, "d", &n1, &m1, &l8); // acc
CreateVar(9, "i", &n1, &m1, &l9); // minpts
CreateVar(10, "i", &n1, &m1, &l10); // ifail
// Assign scilab output variables
//================================
*stk(l7) = finval;
*stk(l8) = acc;
*istk(l9) = (int) minpts;
*istk(l10) = (int) ifail.code;
// Assign order for output variables
//===================================
LhsVar(1) = 7;
LhsVar(2) = 8;
LhsVar(3) = 9;
LhsVar(4) = 10;
END:
// Deallocate memory
//===================
if (a) NAG_FREE(a);
if (b) NAG_FREE(b);
return(0);
}
static double f(Integer n, double z[], Nag_User *comm)
{
// Function describing the integrand for the NAG routine used above
//==================================================================
double tmp_pwr;
tmp_pwr = z[1]+1.0+z[3];
return z[0]*4.0*z[2]*z[2]*exp(z[0]*2.0*z[2])/(tmp_pwr*tmp_pwr);
}
Points to note about this code:
To compile and link this function, you can run the build script nag_builder4a.sce which links this interface to the function name fname that will be seen in Scilab. This must be run in the same directory as the interface file, and when run successfully, generates a loader script loader.sce that needs to be executed in order to dynamically link the newly created library into Scilab.
exec nag_builder4a.sce
exec loader.sce
If the function built and the loader.sce file loaded the function without problem, then we can use the function within Scilab, either on the command-line or in a script. Here we run the example script, nag_example4a.sce.
--> exec nag_example4a.sce
finval = 0.5753623
Tip: If you get any error messages in Scilab check the troubleshooting section.
The only thing that is different between the two source codes nag_intext4a.c and nag_intext4b.c is that in the second one, instead of the user-supplied function being written in C, we have an interface to the user-supplied function written in Scilab. This is so that any variable types can be converted and so that the interface to the NAG routine is correct. Below we just show the part of the file containing the new interface, where we have named the function fSci.
static double fSci(Integer n, double z[], Nag_User *comm);
static double fSci(Integer n, double z[], Nag_User *comm)
{
int mlhs, mrhs, ibegin;
int m_u, n_u, l_u;
int m_z, n_z, l_z;
int i;
double retval;
char name[] = "fUser";
// Put input variables for scilab function onto stack in the right order
// On input u (stored in space 7) is a dummy input to be overwritten with
// the output from fUser
// The inputs for fUser must be in the correct order in the last places of
// the stack. They are then scrapped and the outputs are created starting
// at the same place as the first input variable.
m_u = 1;
n_u = 1;
m_z = (int) n;
n_z = 1;
CreateVar(7, "d", &m_u, &n_u, &l_u); // first input
CreateVar(8, "d", &m_z, &n_z, &l_z); // second input
// Write the input vector z onto the stack in the second input
for (i=0; i<m_z; i++)
{
*stk(l_z+i) = z[i];
}
ibegin = 7; // the first input variable on the stack
mlhs = 1; // number of left hand side outputs to write onto the stack
mrhs = 2; // number of right hand side inputs to read off the stack
SciString(&ibegin,name,&mlhs,&mrhs); // Call scilab function 'name'
// output from scilab function is now on the stack so return the result
retval = *stk(l_u);
return retval;
}
Points to note about this code (other than in the previous code):
deff('u=fUser(uu,z)',
['tmp_pwr = z(2)+1.0+z(4)';
'u = z(1)*4.0*z(3)*z(3)*exp(z(1)*2.0*z(3))/(tmp_pwr*tmp_pwr)']);
To compile and link this function, you need to run the build script nag_builder4b.sce which links this interface to the function name fname that will be seen in Scilab. This must be run in the same directory as the interface file, and when run successfully, generates a loader script loader.sce that needs to be executed in order to dynamically link the newly created library into Scilab.
exec nag_builder4b.sce
exec loader.sce
If the function built and the loader.sce file loaded the function correctly, then we can use the function within Scilab, either on the command-line or in a script. Here we run the example script, nag_example4b.sce.
--> exec nag_example4b.sce
finval = 0.5753623
Tip: If you get any error messages in Scilab check the troubleshooting section.