/* nag_matop_real_gen_matrix_fun_num (f01elc) Example Program.
 *
 * Copyright 2014 Numerical Algorithms Group.
 *
 * Mark 24, 2013.
 */

#include <nag.h>
#include <nag_stdlib.h>
#include <nagf01.h>
#include <nagx02.h>
#include <nagx04.h> 
#include <math.h>

#ifdef __cplusplus
extern "C" {
#endif
   static void NAG_CALL f(Integer *iflag, Integer n, const Complex z[],
                          Complex fz[], Nag_Comm *comm);
#ifdef __cplusplus
}
#endif 

int main(void)
{ 
   /* Scalars */
   Integer        exit_status = 0;
   Integer        i, iflag, j, n, pda;
   double         imnorm;

   /* Arrays */
   static double ruser[1] = {-1.0};
   double         *a = 0;

   /* Nag Types */
   Nag_Comm       comm;
   Nag_OrderType  order;
   NagError       fail; 

   INIT_FAIL(fail);

   #define A(I, J)  a[(J-1)*pda + I-1]

   order = Nag_ColMajor;

   /* Output preamble */
   printf("nag_matop_real_gen_matrix_fun_num (f01elc) ");
   printf("Example Program Results\n\n");
   fflush(stdout);

  /* For communication with user-supplied functions: */
  comm.user = ruser;

   /* Skip heading in data file */
   scanf("%*[^\n]");

   /* Read in the problem size */
   scanf("%ld%*[^\n]", &n);

   pda = n; 
   if (!(a = NAG_ALLOC(pda*n, double))) {
     printf("Allocation failure\n");
     exit_status = -1;
     goto END;
   }

   /* Read in the matrix a from data file */
   for (i = 1; i <= n; i++) 
     for (j = 1; j <= n; j++) scanf("%lf", &A(i, j));
   scanf("%*[^\n]"); 
   /* Find the matrix function using
    * nag_matop_real_gen_matrix_fun_num (f01elc)
    * Function of a real matrix
    */
   nag_matop_real_gen_matrix_fun_num(n, a, pda, f, &comm, &iflag, &imnorm,
                                     &fail);
   if (fail.code != NE_NOERROR) {
     printf("Error from nag_matop_real_gen_matrix_fun_num (f01elc)\n%s\n",
            fail.message);
     exit_status = 1;
     goto END;
   }

   if (fabs(imnorm) > pow(nag_machine_precision,0.8)) {
     printf("\nWARNING: the error estimate returned in imnorm is larger than"
            " expected;\n");
     printf("         the matrix solution printed below is suspect:\n");
     printf("         imnorm = %13.4e.\n\n", imnorm);
   }

   /* Print solution using
    * nag_gen_real_mat_print (x04dac)
    * Print real general matrix (easy-to-use)
    */
   nag_gen_real_mat_print(order, Nag_GeneralMatrix, Nag_NonUnitDiag, n, n, a,
                          pda, "f(A) = cos(2A)", NULL, &fail);
   if (fail.code != NE_NOERROR) {
     printf("Error from nag_gen_real_mat_print (x04dac)\n%s\n", fail.message);
     exit_status = 2;
     goto END;
   }

  END:
   NAG_FREE(a);
   return exit_status;
}

static void NAG_CALL f(Integer *iflag, Integer n, const Complex z[],
                       Complex fz[], Nag_Comm *comm)
{
  /* Scalars */
  Integer j;
#pragma omp master
  if (comm->user[0] == -1.0)
    {
      printf("(User-supplied callback f, first invocation.)\n");
      fflush(stdout);
      comm->user[0] = 0.0;
    }
  for (j = 0; j < n; j++) {
    /* Complex representation of cos 2z */
    fz[j].re =  cos(2.0*z[j].re)*cosh(2.0*z[j].im);
    fz[j].im = -sin(2.0*z[j].re)*sinh(2.0*z[j].im);
  }
  /* Set iflag nonzero to terminate execution for any reason. */
  *iflag = 0;
}