/* The easiest way to pull in all necessary definitions is to include oct.h
   Before we can do it, we need to overwrite Complex type duplicated 
   in Octave ang NAG header files */
#undef Complex
#define Complex OctComplex
#include <octave/oct.h>
/* For NAG definitions we need to include the appropriate NAG header files.
   Again, Complex type has to be overwritten */
#undef Complex
#define Complex NagComplex
// MatrixType is another duplicated definition
#define MatrixType NagMatrixType
#include <nag.h>
#include <nag_stdlib.h>
#include <nagf07.h>

/* The third parameter (nargout) is not used, so it is
   omitted from the list of arguments to DEFUN_DLD in order to avoid
   the warning from gcc about an unused function parameter. */
DEFUN_DLD (nag_cmplx_lineqn, args, ,
   "Calls nag_zgetrf (f07arc) followed by \n\
nag_zgetrs (f07asc).\n\
nag_zgetrs solves a complex system of linear equations with multiple \n\
right-hand sides, AX=B, A^TX=B or A^HX=B, where A has been factorized \n\
by nag_zgetrf.\n")
{
  // Variable to store function output values
  octave_value_list retval;
  int nargin = args.length ();
  if (nargin != 7)
    {
      error ("Insufficient input arguments.");
    }
  else
    {
      // Retrieve input arguments from args
      std::string transpose = args(0).string_value ();
      Integer n = args(1).int_value();
      Integer nrhs = args(2).int_value();
      ComplexMatrix a(args(3).complex_matrix_value());
      Integer pda = args(4).int_value();
      ComplexMatrix b(args(5).complex_matrix_value());
      Integer pdb = args(6).int_value();
      if (! error_state)
        {
          // Declare local variables
#define A(I,J) a_nag[I*pda + J]
#define B(I,J) b_nag[I*pdb + J]
          int i, j;
          NagComplex *a_nag=0, *b_nag=0;
          Integer *ipiv=0;
          Nag_TransType trans;
          Nag_OrderType order=Nag_RowMajor;
          NagError fail;
          INIT_FAIL(fail);
          // Assign the appropriate value of trans
          if (transpose.compare("Nag_NoTrans") == 0) trans=Nag_NoTrans;
          else if (transpose.compare("Nag_Trans") == 0) trans=Nag_Trans;
          else if (transpose.compare("Nag_ConjTrans") == 0) trans=Nag_ConjTrans;
          else error ("transpose not Nag_NoTrans, Nag_Trans or Nag_ConjTrans");
          
          if ( !( ipiv = NAG_ALLOC(n, Integer))  ||
               !( a_nag = NAG_ALLOC(n*pda, NagComplex)) ||
               !( b_nag = NAG_ALLOC(n*pdb, NagComplex)) )
            {
              error ("Allocation failure\n");
              goto END;
            }
          
          /* NAG Complex type differs from Octave Complex type, so we need to 
             copy the input OctComplex matrices into NagComplex arrays */
          for (i = 0; i < n; i++)
            for (j = 0; j < pda; j++)
              {
                A(i,j).re = a.checkelem(i,j).real();
                A(i,j).im = a.checkelem(i,j).imag();
              } 
          for (i = 0; i < n; i++)
            for (j = 0; j < pdb; j++)
              {
                B(i,j).re = b.checkelem(i,j).real();
                B(i,j).im = b.checkelem(i,j).imag();
              }
      
          // Call NAG routines
          /* Factorize A */
          /* nag_zgetrf (f07arc).
           * LU factorization of complex m by n matrix
           */
          nag_zgetrf(order,n,n,a_nag,pda,ipiv,&fail);
          if (fail.code!=NE_NOERROR)
            {
              error ("nag_zgetrf (f07arc).\n%s\n", fail.message);
              goto END;
            }
          /* Compute solution */
          /* nag_zgetrs (f07asc).
           * Solution of complex system of linear equations, multiple
           * right-hand sides, matrix already factorized by nag_zgetrf
           * (f07arc)
           */
          nag_zgetrs(order,trans,n,nrhs,a_nag,pda,ipiv,b_nag,pdb,&fail);
          if (fail.code!=NE_NOERROR)
            {
              error ("nag_zgetrs (f07asc).\n%s\n",fail.message);
              goto END;
            }
          /* The result is store in NAG Complex type, so we need to 
             copy the output NagComplex array back into OctComplex matrix */
          for (i = 0; i < n; i++)
            for (j = 0; j < pdb; j++)
              {
                b.checkelem(i,j).real() = B(i,j).re;
                b.checkelem(i,j).imag() = B(i,j).im; 
              }
          if (a_nag) NAG_FREE(a_nag);
          if (b_nag) NAG_FREE(b_nag);
          if (ipiv) NAG_FREE(ipiv);
          // Assign output arguments to retval
          retval(0) = b;
          retval(1) = fail.code;
        }
    }
  
 END:
  return retval;
}