/* nag_matop_complex_gen_matrix_actexp_rcomm (f01hbc) Example Program.
 *
 * Copyright 2017 Numerical Algorithms Group.
 *
 * Mark 26.2, 2017.
 */

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

static void matmul(Nag_TransType tran, Integer n, Integer m, Complex a[],
                   Integer icolzp[], Integer irowix[], Complex b[],
                   Complex c[]);

int main(void)
{

  /* Scalars */
  Integer exit_status = 0, irevcm = 0;
  Integer i, j, m, n, nnz, ldb, ldb2, ldx, ldy;
  Complex t, tr;

  /* Arrays */
  Complex *a = 0, *b = 0, *b2 = 0, *ccomm = 0;
  Complex *p = 0, *r = 0, *x = 0, *y = 0, *z = 0;
  double *comm = 0;
  Integer *icolzp = 0, *irowix = 0, *icomm = 0;

  /* Nag Types */
  NagError fail;
  Nag_OrderType order = Nag_ColMajor;
  Nag_TransType tran = Nag_ConjTrans, notran = Nag_NoTrans;

  INIT_FAIL(fail);

#define B(I, J) b[(J-1)*ldb + I-1]

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

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

  /* Read in the problem size and the value of the parameter t */
  scanf("%" NAG_IFMT " %" NAG_IFMT " (%lf , %lf ) %*[^\n]", &n, &m, &t.re,
        &t.im);

  ldb = ldb2 = ldx = ldy = n;

  if (!(b = NAG_ALLOC(n * m, Complex)) ||
      !(b2 = NAG_ALLOC(n * m, Complex)) ||
      !(ccomm = NAG_ALLOC(n * (m + 2), Complex)) ||
      !(p = NAG_ALLOC(n, Complex)) ||
      !(r = NAG_ALLOC(n, Complex)) ||
      !(x = NAG_ALLOC(n * 2, Complex)) ||
      !(y = NAG_ALLOC(n * 2, Complex)) ||
      !(z = NAG_ALLOC(n, Complex)) ||
      !(comm = NAG_ALLOC(3 * n + 14, double)) ||
      !(icolzp = NAG_ALLOC(n + 1, Integer)) ||
      !(icomm = NAG_ALLOC(2 * n + 40, Integer)))
  {
    printf("Allocation failure\n");
    exit_status = -1;
    goto END;
  }

  /* Read in the sparse matrix a in compressed column format */
  for (i = 0; i < n + 1; ++i)
    scanf("%" NAG_IFMT "", &icolzp[i]);
  scanf("%*[^\n]");

  nnz = icolzp[n] - 1;

  if (!(a = NAG_ALLOC(nnz, Complex)) || !(irowix = NAG_ALLOC(nnz, Integer)))
  {
    printf("Allocation failure\n");
    exit_status = -2;
    goto END;
  }

  for (i = 0; i < nnz; ++i)
    scanf(" ( %lf , %lf ) %" NAG_IFMT "%*[^\n]", &a[i].re, &a[i].im,
          &irowix[i]);

  /* Read in the matrix b from data file */
  for (i = 1; i <= n; i++)
    for (j = 1; j <= m; j++)
      scanf(" ( %lf , %lf ) ", &B(i, j).re, &B(i, j).im);
  scanf("%*[^\n]");

  /* Compute the trace of A */
  tr.re = 0.0;
  tr.im = 0.0;
  j = 1;
  for (i = 0; i < nnz; ++i) {
    /* new column? */
    if (icolzp[j] <= i + 1)
      j++;
    /* diagonal element? */
    if (irowix[i] == j) {
      tr.re += a[i].re;
      tr.im += a[i].im;
    }
  }

  /* Find exp(tA) B using
   * nag_matop_complex_gen_matrix_actexp_rcomm (f01hbc)
   * Action of the exponential of a complex matrix on a complex matrix
   */
  while (1) {
    nag_matop_complex_gen_matrix_actexp_rcomm(&irevcm, n, m, b, ldb, t, tr,
                                              b2, ldb2, x, ldx, y, ldy, p, r,
                                              z, ccomm, comm, icomm, &fail);

    switch (irevcm) {
    case 0:
      /* Final exit. */
      goto END_REVCM;
      break;
    case 1:
      /* Compute AB and store in B2 */
      matmul(notran, n, m, a, icolzp, irowix, b, b2);
      break;
    case 2:
      /* Compute AX and store in Y */
      matmul(notran, n, 2, a, icolzp, irowix, x, y);
      break;
    case 3:
      /* Compute A^H Y and store in X */
      matmul(tran, n, 2, a, icolzp, irowix, y, x);
      break;
    case 4:
      /* Compute AZ and store in P */
      matmul(notran, n, 1, a, icolzp, irowix, z, p);
      break;
    case 5:
      /* Compute A^H Z and store in R */
      matmul(tran, n, 1, a, icolzp, irowix, z, r);
    }
  }

END_REVCM:
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_matop_complex_gen_matrix_actexp_rcomm (f01hbc)\n"
           "%s\n", fail.message);
    exit_status = 2;
    goto END;
  }

  /* Print solution using
   * nag_gen_complx_mat_print (x04dac)
   * Print complex general matrix (easy-to-use)
   */
  nag_gen_complx_mat_print(order, Nag_GeneralMatrix, Nag_NonUnitDiag, n, m,
                           b, ldb, "exp(tA) B", NULL, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_gen_complx_mat_print (x04dac)\n%s\n",
           fail.message);
    exit_status = 3;
    goto END;
  }

END:
  NAG_FREE(a);
  NAG_FREE(b);
  NAG_FREE(b2);
  NAG_FREE(comm);
  NAG_FREE(ccomm);
  NAG_FREE(p);
  NAG_FREE(r);
  NAG_FREE(x);
  NAG_FREE(y);
  NAG_FREE(z);
  NAG_FREE(icomm);
  NAG_FREE(icolzp);
  NAG_FREE(irowix);
  return exit_status;
}

static void matmul(Nag_TransType tran, Integer n, Integer m, Complex a[],
                   Integer icolzp[], Integer irowix[], Complex b[],
                   Complex c[])
{
  Integer i, ir, j, k, l;
  Complex a1, a2, t1;

  for (j = 0; j < m * n; j++) {
    c[j].re = 0.0;
    c[j].im = 0.0;
  }
  if (tran == Nag_NoTrans) {
    l = -1;
    for (j = 0; j < m; j++) {
      for (i = 0; i < n; i++) {
        l++;
        t1.re = b[l].re;
        t1.im = b[l].im;
        for (k = icolzp[i] - 1; k < icolzp[i + 1] - 1; k++) {
          ir = j * n + irowix[k] - 1;
          /* Evaluate t1*a[k] using nag_complex_multiply (a02ccc). */
          a1 = nag_complex_multiply(t1, a[k]);
          c[ir].re += a1.re;
          c[ir].im += a1.im;
        }
      }
    }
  }
  else {
    l = -1;
    for (j = 0; j < m; j++) {
      for (i = 0; i < n; i++) {
        l++;
        t1.re = 0.0;
        t1.im = 0.0;
        for (k = icolzp[i] - 1; k < icolzp[i + 1] - 1; k++) {
          ir = j * n + irowix[k] - 1;
          /* Evaluate conjg(a[k]) using nag_complex_conjg (a02cfc). */
          a2 = nag_complex_conjg(a[k]);
          /* Evaluate conjg(a[k])*b[ir] using nag_complex_multiply (a02ccc). */
          a1 = nag_complex_multiply(a2, b[ir]);
          t1.re += a1.re;
          t1.im += a1.im;
        }
        c[l].re += t1.re;
        c[l].im += t1.im;
      }
    }
  }
}