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