# Calling NAG Library Routines from Scilab - Example 3

## 3.3. Example 3

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

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
//******************************************************************

//*============*
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;
}
}

//===========
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);
}

//==============
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;
}
}
}

//============
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);
}

//============
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);
}

//==============
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);
}



• 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

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.