Calling NAG Library Routines from Scilab - Example 3
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 4   Skip to Example 5   Reference Page

3.3. Example 3

eigvals

Complex Eigenvalalues and Eigenvectors of a Complex General Matrix (f02gcc)

Here we show how to write an interface to a NAG C Library function with complex input and output matrix variables. The function we will demonstrate with is f02gcc which computes the eigenvalues and eigenvectors of a complex system of equations.

Contents

  1. Function prototype from the NAG C Library Manual
  2. C Interface
  3. Compiling with the builder script
  4. Example of calling the function

  1. Function prototype from the NAG C Library Manual

    According to the C Library Manual, the prototype for function f02gcc looks like this:

         #include <nag.h>

     #include <nagf02.h>

     void nag_complex_eigensystem_sel(Nag_Select_Eigenvalues crit, Integer n, Complex a[], Integer tda, double wl,
                                      double wu, Integer mest, Integer *m, Complex w[], Complex v[], Integer tdv,
                                      NagError *fail)
    This function takes an input matrix of type Complex amongst other inputs and computes the selected eigenvalues and their corresponding eigenvectors.

  2. C Interface

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

        /* Example 3
       =========

      Shows how to implement a Scilab wrapper calling in Complex matrix
      arguments and characters */

    #include "stack-c.h"
    #undef Complex
    #define Complex NagComplex
    #include <nag.h>
    #include <nag_stdlib.h>
    #include <nagf02.h>

    #define SciComplex doublecomplex

    #include "stdio.h"

    // Macros to allow matrix notation in the code.
    // Here the matrices are accessed in row order due to working in C
    #define A_NAG(I,J) a_nag[(I)*n + J]
    #define AMAT(I,J) amat[I*n + J]
    #define V_NAG(I,J) v_nag[(I)*n + J]
    #define VMAT(I,J) vmat[(I)*n + J]

    int nag_intext3(char *fname)
    {
      // to call this function in scilab use:
      // [w,v,m,ifail] = nag_cmplx_eisys_sel_fun(crit,n,a,wl,wu,mest)
      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, tda, mest,i,j, tdv;
      Integer m;
      double wl, wu;
      char crit;
      SciComplex *amat=0, *vmat=0, *wmat=0;
      NagComplex *a_nag=0, *v_nag=0, *w_nag=0;
      Nag_Select_Eigenvalues crit_nag;
      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, "c", &m1, &n1, &l1); // input crit
      GetRhsVar(2, "i", &m2, &n2, &l2); // input n
      GetRhsVar(3, "z", &m3, &n3, &l3); // input a
      GetRhsVar(4, "d", &m4, &n4, &l4); // input wl
      GetRhsVar(5, "d", &m5, &n5, &l5); // input wu
      GetRhsVar(6, "i", &m6, &n6, &l6); // input mest

      // Allocate memory for a_nag and amat
      if ( !(a_nag = NAG_ALLOC(m3*n3,NagComplex) ) ||
           !(amat = NAG_ALLOC(m3*n3,SciComplex)))
        {
          sciprint("Allocation failure\n");
          Error(999); goto END;
        }

      //******************************************************************
      // Read in input variables from Scilab and convert to NAG variables
      //******************************************************************

      // Read in crit
      //*============*
      if (m1!=1 || n1!=1)
        {
          sciprint("%s: Dimension should be 1x1 character for arg 1\r\n",fname);
          Error(999); goto END;
        }
      else
        {
          crit = *cstk(l1);
          if (crit == 'r' || crit == 'R')
            {
              crit_nag = Nag_Select_RealPart;
            }
          else if( crit == 'm' || crit == 'M')
            {
              crit_nag = Nag_Select_Modulus;
            }
          else
            {
              sciprint("%s: Arg 1 should be either 'r' or 'm'. \r\n",fname);
              Error(999); goto END;
            }
        }

      // Read in n
      //===========
      if (m2!=1 || n2!=1)
        {
          sciprint("%s: Dimension should be 1x1 for arg 2\r\n",fname);
          Error(999); goto END;
        }
      else
        {
          n = *istk(l2);
        }

      // Read in amat
      //==============
      if (m3!=n || n3!=n)
        {
          sciprint("%s: Dimension should be nxn for arg 3\r\n",fname);
          Error(999); goto END;
        }
      else
        {
          // Convert Complex data types Scilab to NAG
          // Matrices are stored in column order in Scilab
          for(i=0; i<n; i++)
            {
              for(j=0; j<n; j++)
                {
                  AMAT(i,j) = *zstk(l3 + i + j * n);
                  A_NAG(i,j).re = AMAT(i,j).r;
                  A_NAG(i,j).im = AMAT(i,j).i;
                }
            }
        }

      // Read in wl
      //============
      if (m4!=1 || n4!=1)
        {
          sciprint("%s: Dimension should be 1x1 for arg 4\r\n",fname);
          Error(999); goto END;
        }
      else
        {
          wl = *stk(l4);
        }

      // Read in wu
      //============
      if (m5!=1 || n5!=1)
        {
          sciprint("%s: Dimension should be 1x1 for arg 5\r\n",fname);
          Error(999); goto END;
        }
      else
        {
          wu = *stk(l5);
        }

      // Read in mest
      //==============
      if (m6!=1 || n6!=1)
        {
          sciprint("%s: Dimension should be 1x1 for arg 6\r\n",fname);
          Error(999); goto END;
        }
      else
        {
          mest = *istk(l6);
        }
      //***********************************
      // End of reading in input variables
      //***********************************

      // Allocate memory for NAG output v_nag(n,mest) and w_nag
      //========================================================
      if ( !(v_nag = NAG_ALLOC(n*mest,NagComplex) ))
        {
          sciprint("Allocation failure\n");
          Error(999); goto END;
        }

      if ( !(w_nag = NAG_ALLOC(MAX(1,n),NagComplex) ))
        {
          sciprint("Allocation failure\n");
          Error(999); goto END;
        }

      // Define the other arguments for the NAG routine
      //================================================
      tda = n;
      tdv = mest;

      // Call NAG routine f02gcc
      //=========================
      nag_complex_eigensystem_sel(crit_nag,n,a_nag,tda,wl,wu,mest,&m,w_nag,v_nag,tdv,&ifail);



      // Create output variables for Scilab
      //====================================
      CreateVar(7, "z", &n, &mest, &l7); // vmat
      CreateVar(8, "z", &n, &m1, &l8); // wmat
      CreateVar(9, "i", &n1, &m1, &l9); // m
      CreateVar(10, "i", &n1, &m1, &l10); // ifail

      // Allocate memory for output vmat(n,m) and wmat
      //===============================================
      if ( !(vmat = NAG_ALLOC(n*mest,SciComplex)) ||
           !(wmat = NAG_ALLOC(MAX(1,n),SciComplex)))
        {
          sciprint("Allocation failure\n");
          Error(999); goto END;
        }

      // Assign scilab output variables
      //================================

      // Convert complex types NAG to Scilab
      for (i=0; i<n; i++)
        {
          for (j=0; j<mest; j++)
            {
              VMAT(i,j).r = V_NAG(i,j).re;
              VMAT(i,j).i = V_NAG(i,j).im;
              *zstk(l7 +i + j*n) = VMAT(i,j);
            }
        }
      for (i=0; i<n; i++)
        {
          wmat[i].r = w_nag[i].re;
          wmat[i].i = w_nag[i].im;
          *zstk(l8+i) = wmat[i];
        }
      *istk(l9) = m;
      *istk(l10) = ifail.code;

      // Print NAG error message if routine fails
      //==========================================
      if (ifail.code)
        {
          sciprint("ifail message %s",ifail.message);
        }

      // Assign order for output variables
      //===================================
      LhsVar(1) = 8;
      LhsVar(2) = 7;
      LhsVar(3) = 9;
      LhsVar(4) = 10;

      // Free allocated memory
      //=======================
     END:
      if (a_nag) NAG_FREE(a_nag);
      if (amat) NAG_FREE(amat);
      if (v_nag) NAG_FREE(v_nag);
      if (vmat) NAG_FREE(vmat);
      if (w_nag) NAG_FREE(w_nag);
      if (wmat) NAG_FREE(wmat);

      return(0);
    }

    Points to note about this code:

    • Matrices are stored in column order on the Scilab stacks which is the same as Fortran code but not C code, where matrices are stored in row order.
    • We demonstrate the use of taking in a character argument as an option. This then gets changed to a NAG enum type nag_select_eigenvalues.
    • The Scilab complex type is originally doublecomplex and is a struct with .r being the real part and .i being the imaginary part. The NAG Complex type on the other hand stores the real part and imaginary part as .re and .im. Therefore we have to convert between the two different Complex types in this interface.
    • 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 zstk()

  3. Compiling with the Builder Script

    To compile and link this function, you can run the build script nag_builder3.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_builder3.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_example3.sce.

        --> exec nag_example3.sce
      v  =

      - 0.3865491 + 0.1732346i  - 0.0356136 - 0.1782180i
      - 0.3539288 + 0.4528810i    0.1263743 + 0.2666324i
        0.6123701                 0.0129326 - 0.2965682i
      - 0.0859284 - 0.3283626i    0.8898240

      w  =

      - 5.0000335 + 2.0060272i
      3.0022643 - 3.9998187i
      7.9981945 - 0.9963651i
      - 6.0004253 - 6.9998434i

    Tip: If you get any error messages in Scilab check the troubleshooting section.