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