NAG Technical Report 5/2009

Calling NAG Library Routines from Scilab


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

3.5. Example 5

pseudorand

Pseudorandom Numbers taken from the Normal Distribution

We show how to call a routine from the NAG Fortran Library: G05LAF Although we are calling a NAG Fortran Library routine, we must use C code for the glue between Scilab and NAG. In this example we will demonstrate calling a random number generator which returns a vector of pseudorandom numbers generated from the Normal distribution. We will include an option to generate either a repeatable or non-repeatable sequence of random numbers.

Contents

  1. Function prototype from the NAG Fortran Library Manual
  2. C interface
  3. Compiling with the Builder Script
  4. Calling the function
  1. Function prototype from the NAG Fortran Library Manual

    According to the Fortran Library Manual, the prototype for function G05LAF looks like this:

      SUBROUTINE G05LAF(XMU, VAR, N, X, IGEN, ISEED, IFAIL)
      integer N, IGEN, ISEED(4), IFAIL
      double precision XMU, VAR, X(N)
    
    See the NAG FORTRAN Library Manual [6] for detailed argument descriptions.
  2. C interface

    Here is the source code of our C interface nag_intext5.c:
    /* Example 5
       =========
    
       Shows how to implement a Scilab wrapper calling a NAG routine from the
       Fortran library */
    
    #include "stack-c.h"
    #undef Complex
    #define Complex NagComplex
    #define SciComplex doublecomplex
    
    #include "stdio.h"
    
    ... (other header file definitions missing here which are in the source code) ...
    
    int nag_intext5(char *fname)
    {
      // to call this function in scilab use:
      // [x,ifail] = nag_psdrnd_norm_fun(n,xmu,var,repeatable)
      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 n, i, min, max;
      int repeatable;
    
      Integer igen, ifail, ndim;
      Integer *iseed=0;
    
      double xmu, var;
      double *x;
    
      // define minimum and maximum left and right hand arguments
      int minlhs=1, minrhs=4, maxlhs=2, maxrhs=4;
      Nbvars = 0;
      CheckRhs(minrhs, maxrhs);
      CheckLhs(minlhs,maxlhs);
    
      GetRhsVar(1, "i", &m1, &n1, &l1); // input n
      GetRhsVar(2, "d", &m2, &n2, &l2); // input xmu
      GetRhsVar(3, "d", &m3, &n3, &l3); // input var
      GetRhsVar(4, MATRIX_OF_BOOLEAN_DATATYPE, &m4, &n4, &l4); // repeatable
    
      // there are also other datatypes for integers, real, string etc.
      // MATRIX_OF_INTEGER_DATATYPE on the istk
      // MATRIX_OF_RATIONAL_DATATYPE on the sstk
      // MATRIX OF DOUBLE DATATYPE on the stk
      // MATRIX_OF_COMPLEX_DATATYPE on the zstk
      // MATRIX_OF_BOOLEAN_DATATYPE on the istk
      // STRING_DATATYPE on the cstk
    
    
      //******************************************************************
      // Read in input variables from Scilab and convert to NAG variables
      //******************************************************************
    
      // Read in n
      //*============*
      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 output array x
      // ==================================
      if ( !(x = malloc(n*sizeof(double))))
        {
          sciprint("Allocation failure\n");
          Error(999); goto END;
    
        }
    
    
      // Read in xmu
      //===========
      if (m2!= 1 || n2!=1)
        {
          sciprint("%s: Dimension should be 1x1 for arg 2\r\n",fname);
          sciprint("n = %i", n);
          Error(999); goto END;
        }
      else
        {
          xmu = *stk(l2);
        }
    
      // Read in var
      //==============
      if (m3!=1 || n3!=1)
        {
          sciprint("%s: Dimension should be 1x1 for arg 3\r\n",fname);
          Error(999); goto END;
        }
      else
        {
          var = *stk(l3);
        }
    
      // Read in repeatable
      //====================
      if (m4!=1 || n4!=1)
        {
          sciprint("%s: Dimension should be 1x1 for arg 4\r\n",fname);
          Error(999); goto END;
        }
      else
        {
          repeatable = *istk(l4);
        }
    
    
      //***********************************
      // End of reading in input variables
      //***********************************
    
    
      if ( !(iseed = malloc(4*sizeof(Integer))))
        {
          sciprint("Allocation failure\n");
          Error(999); goto END;
    
        }
    
      // Call initialisation routine for this rng
      // ========================================
      igen = 0; // Basic generator
    
      iseed[0] = 23;
      iseed[1] = 4;
      iseed[2] = 10;
      iseed[3] = 1985;
    
      // Initialise the rng either repeatably or not repeatably
      if (repeatable)
        {
          g05kbf_(&igen, iseed);
        }
      else
        {
          g05kcf_(&igen, iseed);
        }
    
      // Call random number generator based on Normal PDF
      //* ================================================
      ifail = -1;
    
      g05laf_(&xmu, &var, &ndim, x, &igen, iseed, &ifail);
    
      // Print NAG error message if routine fails
      //=========================================
      if (ifail != 0)
        sciprint("g05laf failed");
    
      // Create output variables for Scilab
      //====================================
      CreateVar(5, "d", &n, &m1, &l5); // x
      CreateVar(6, "i", &n1, &m1, &l6); // ifail
    
      // Assign scilab output variables
      //================================
      for (i=0;i<n;i++)
        {
          *stk(l5+i) = x[i];
        }
    
      *istk(l6) = (int) ifail;
    
      // Assign order for output variables
      //===================================
      LhsVar(1) = 5;
      LhsVar(2) = 6;
    
     END:
    
      // Deallocate memory
      //===================
      if (x) free(x);
      if (iseed) free(iseed);
    
      return(0);
    }
    

    Points to note about this code:

  3. Compiling with the Builder Script

    To compile and link this function, you need to run the build script nag_builder5.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.

    This builder script looks slightly different from the previous ones, as we are linking to the NAG Fortran library.

        exec nag_builder5.sce
        exec loader.sce
    
  4. 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_example5.sce.

        --> exec nag_example5.sce
        ifail  =
               0.
    
          x  =
    
               0.9074485
               1.1109788
               0.8802789
               1.4147883
               2.1821836
               3.5388501
               1.3572214
               3.0131266
               2.5570866
               2.1172563
    

    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 4   Reference Page