```/* nag_matop_complex_gen_matrix_actexp_rcomm (f01hbc) Example Program.
*
* Copyright 2020 Numerical Algorithms Group.
*
* Mark 27.1, 2020.
*/

#include <math.h>
#include <nag.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_file_print_matrix_complex_gen (x04dac)
* Print complex general matrix (easy-to-use)
*/
nag_file_print_matrix_complex_gen(order, Nag_GeneralMatrix, Nag_NonUnitDiag,
n, m, b, ldb, "exp(tA) B", NULL, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_file_print_matrix_complex_gen (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;
}
}
}
}
```