Calling NAG Library Routines from Scilab

Start of report   Skip to Example 1   Skip to Example 2   Skip to Example 3   Skip to Example 5   Reference Page

3.4. Example 4

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.

Contents

1. Function prototype from the NAG C Library Manual
2. Writing the User-supplied Function in C
   1. C Interface
   2. Compiling with the builder script
   3. Example of calling the function
3. Writing the User-Supplied Function in Scilab
   1. C Interface for the Scilab Code
   2. User-supplied functions in Scilab
   3. Compiling with the builder script
   4. Calling the function

1. Function prototype from the NAG C Library Manual

   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)
   

   For more information about the argument list of this function, refer to the NAG C Library manual.

2. Writing the User-supplied Functions in C

   1. C Interface

      Here is the source code of our interface written in C nag_intext4a.c with the user-supplied function 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 m7,n7,l7;
            int m8,n8,l8;
            int m9,n9,l9;
            int m10, n10, 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:

      • The input function for the NAG routine has a specific interface that has to be adhered to, which is documented in the NAG C Library manual [5]
      • The interface doesn't take in the input arguments of the function as its arguments, but only takes fname the name of the function in Scilab. In this example we have chosen fname to be the same as the NAG C library name.
      • The Scilab function GetRhsVar is used to retrieve the input arguments of fname
      • CreateVar creates a variable into which the NAG output is placed
      • The interface doesn't actually return the outputs of the function fname - instead the Scilab array LhsVar stores the variable numbers that are stored on Scilab's stack stk()

   2. Compiling with the Builder Script

      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
      

   3. Calling the function

      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.

3. Writing the user-supplied functions in Scilab

   1. C Interface for the Scilab code

      Here is the source code of our second interface written in C nag_intext4b.c:

      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 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.

          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):

      • The User-defined function must still have the same interface as specified in the C Library manual and therefore we must write an interface to the Scilab function.
      • The Scilab function fUser is not called directly, but through a Scilab function called SciString. The arguments for SciString are
        • ibegin: the variable number of the first input argument. Note that the input arguments must be next to each other and in the order they are to passed into fUser with the first input argument being at ibegin.
        • name: the name of the Scilab function.
        • mlhs and mrhs: the number of input and output variables of the Scilab function. The input arguments are destroyed after use, and the output arguments are then put on the stack in order starting at ibegin.

   2. User-supplied functions in Scilab

      Source code of our Scilab function is in the example script nag_example4b.sce but could equally have been in a separate script. The source also demonstrates two different ways of writing the function.

          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)']);
      

   3. Compiling with the Builder Script

      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
      

   4. Calling the function

      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.

Start of report   Skip to Example 1   Skip to Example 2   Skip to Example 3   Skip to Example 5   Reference Page