Critcl NAG

The following examples show how the critcl package can be used to write wrappers for NAG [1] (Numerical Algorithms Group) library routines.

Here is a simple example which calls s07aaf, a tan(x) function [2].

 #
 # Test  tcl to nag via critcl
 # This was tested on a Solaris.

 package require critcl

 # This list of libraries worked OK for the Solaris installation of NAG.
 critcl::clibraries -lnag -lpdef -lbndy -lfunct1 -lfunct2 -llsfun1 -llsfun2 -lhess2 -llshes2

 critcl::ccode {
     #include <stdio.h> 
     #include <math.h>
 double  s07aaf_(double *, int *);
 }

 critcl::cproc Tan { double x } double {
   
    /* tan(x) function */
     double y;
     double xx;
     int ifail;
     xx = x;
     ifail = 0;

     y = s07aaf_(&xx,&ifail);
    
    return y;
 }

Save the above in a file tan.tcl. Test by running interactive tclsh shell (x is in Radians):-

 $ tclsh
 % source tan.tcl
 % Tan -.5
 -0.546302489844
 %

cproc is OK for simple routines with simple data types. The following code calls the NAG routine e02acf [3] to perform a minimax curve fit by polynomials

 #
 # Test  tcl to nag
 #

 package require critcl

 # Different linkage for Linux and Solaris.
 if { $tcl_platform(os) == "Linux" } {
     critcl::clibraries -lnag -lm  -lpgc -lpgftnrtl
 } else {
     critcl::clibraries -lnag -lpdef -lbndy -lfunct1 -lfunct2 -llsfun1 -llsfun2 -lhess2 -llshes2
 }

 critcl::ccode {
    #include <stdio.h> 
    #include <math.h>

 double  e02acf_(double *, double*, int *, double *, int *, double *);
 }


 critcl::ccommand nagcfit { dummy ip objc objv } {

    /* objv[1] list of x values */
    /* objv[2] list of y values */
    /* objv[3] int Degree of fit */

    int i, result ;
    int xlen, ylen ; 
    unsigned char *data;
    char errmsg[80];
    Tcl_Obj *obj = Tcl_GetObjResult(ip);
    Tcl_Obj *listElemPtr;
    Tcl_Obj **objCoeff;
    Tcl_Obj *tobj;
    int degree ;
    double *xdat;
    double *ydat;
    int npoints;
    double *coeff;
    double ref;

    /* Check number of args */
    if (objc <  4|| objc > 4) {
        Tcl_WrongNumArgs(ip, 1, objv, "{list of X} {list of Y} degree+1");
        return TCL_ERROR ;
    }
    /* get list of x value */
    result = Tcl_ListObjLength(ip, objv[1], &xlen);
    if (result != TCL_OK) {
        return TCL_ERROR ;
    }

    /* get list of y values */
    result = Tcl_ListObjLength(ip, objv[2], &ylen);
    if (result != TCL_OK) {
        return TCL_ERROR ;
    }

    /* Find the min number of xy points */
    if (xlen < ylen ) {
        npoints = xlen;
    } else {
        npoints = ylen;
    }

    /* Get int degree of the Polynomial fit */
    if (TCL_OK != Tcl_GetIntFromObj(ip, objv[3],&degree)) {
      return TCL_ERROR ; 
    }

    /* Range check degree > 1 ||  < npoints-1 or < 100 */ 
    if (degree < 1 || degree >= npoints || degree > 100) {
        Tcl_AppendStringsToObj(obj, " Degree must > 0, < number of data points+2 and less than 100", (char *) NULL);
        return TCL_ERROR ;
    }
    degree++; /* Bump degree to degree +1 */

    /* Allocate Storage for data and results */
    xdat = (double*) ckalloc(npoints *  sizeof(double));
    ydat = (double*) ckalloc(npoints *  sizeof(double));
    coeff = (double*) ckalloc(degree *  sizeof(double));
    objCoeff = (Tcl_Obj**) ckalloc(degree * sizeof(Tcl_Obj));
    
    /* Transfer the x and y value lists into arrays */

    for (i=0; i< npoints; i++) {
      /* Get pointer to next data x value */
      Tcl_ListObjIndex(ip, objv[1], i, &listElemPtr);
      /* Extract the double value from it. */
      result=Tcl_GetDoubleFromObj(ip,listElemPtr,&xdat[i]);

      /* Error if data item cannot be converted into  double */
      if (result != TCL_OK) {
          Tcl_AppendStringsToObj(obj, " in x data.", (char *) NULL);
          return TCL_ERROR ;
      }

      /* Check that x data value a accending order */
      if ( i > 0 ) {
          if ( xdat[i] < xdat[i-1] ) {
              sprintf(errmsg," X data must be accending order, see items %d and %d ",i-1,i); 
              Tcl_AppendStringsToObj(obj, errmsg, (char *) NULL);
              return TCL_ERROR ;
          }
      }
      /* Get pointer to next data y value */
      Tcl_ListObjIndex(ip, objv[2], i, &listElemPtr);
      /* Extract the double value from it. */
      result=Tcl_GetDoubleFromObj(ip,listElemPtr,&ydat[i]);
      /* Error if data item cannot be converted into  double */
      if (result != TCL_OK) {
          Tcl_AppendStringsToObj(obj, " in y data.", (char *) NULL);
          return TCL_ERROR ;
      }
    }
    #if 0
    for (i=0; i< npoints; i++) {
        printf("x[i]= %e, y[i]= %e\n",xdat[i],ydat[i]);
    }
    #endif
    /* The nag routine */
    e02acf_(xdat,ydat,&npoints,coeff,&degree,&ref);
    
    /* Copy resulting Coefficient array to new double objects */
    for (i=0; i < degree; i++) {
        objCoeff[i] = Tcl_NewDoubleObj(coeff[i]);
    } 

    /* This proc result object to be a list of double objects */
    Tcl_SetListObj(obj, degree, objCoeff);

    /* Free malloc space */
    ckfree(xdat);
    ckfree(ydat);
    ckfree(coeff);
    ckfree(objCoeff);

    return TCL_OK;
 }

Save the above in a file called nagcfit.tcl Again we can test interactively with tclsh:

 $ tclsh
 % source nagcfit.tcl
 % set x [ list 0.0 1.0 2.0 3.0 3.5 4.0 5.0 6.0 7.0 ]
 0.0 1.0 2.0 3.0 3.5 4.0 5.0 6.0 7.0
 % set y [ list 10.0 8.0 6.0 5.0 4.5 5.0 6.0 8.0 10.0 ]
 10.0 8.0 6.0 5.0 4.5 5.0 6.0 8.0 10.0
 % nagcfit $x $y 2
 10.3469387755 -3.14285714286 0.448979591837

The parameters passed to nagcfit are a list if x values, a list of y values and a integer degree of the polynomial to be fitted. It returns a list coefficients (cal(0) cal(1) ...cal(degree) for a polynomial.

Category Critcl

Category Mathematics