```/* nag_matop_real_gen_matrix_actexp_rcomm (f01gbc) Example Program.
*
* Copyright 2014 Numerical Algorithms Group.
*
* Mark 24, 2013.
*/

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

int main(void)
{

/* Scalars */
Integer       exit_status = 0, irevcm = 0;
Integer       i, j, m, n, nnz, ldb, ldb2, ldx, ldy;
double        one = 1.0, zero = 0.0;
double        t, tr;

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

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

INIT_FAIL(fail);
INIT_FAIL(fail2);

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

/* Output preamble */
printf("nag_matop_real_gen_matrix_actexp_rcomm (f01gbc) ");
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("%ld %ld %lf %*[^\n] ", &n, &m, &t);

ldb = ldb2 = ldx = ldy = n;

if (!(b = NAG_ALLOC(n*m, double)) ||
!(b2 = NAG_ALLOC(n*m, double)) ||
!(comm = NAG_ALLOC(n*m+3*n+12, double)) ||
!(p = NAG_ALLOC(n, double)) ||
!(r = NAG_ALLOC(n, double)) ||
!(x = NAG_ALLOC(n*2, double)) ||
!(y = NAG_ALLOC(n*2, double)) ||
!(z = NAG_ALLOC(n, 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("%ld", &icolzp[i]);
scanf("%*[^\n]");

nnz = icolzp[n] - 1;

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

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

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

/* Compute the trace of A */
tr = 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 += a[i];
}

/* Find exp(tA) B using
* nag_matop_real_gen_matrix_actexp_rcomm (f01gbc):
*     Action of the exponential of a complex matrix on a complex matrix.
*
* Sparse matrix multiplies are performed using
* nag_superlu_matrix_product (f11mkc):
*     Real sparse nonsymmetric matrix matrix multiply, compressed column
*     storage.
*/
while (1) {
nag_matop_real_gen_matrix_actexp_rcomm(&irevcm, n, m, b, ldb, t, tr, b2,
ldb2, x, ldx, y, ldy, p, r, z, comm,
icomm, &fail);
switch(irevcm) {
case 0:
/* Final exit.*/
goto END_REVCM;
break;
case 1 :
/* Compute AB and store in B2 */
nag_superlu_matrix_product(order, notran, n, m, one, icolzp, irowix, a, b,
ldb, zero, b2, ldb2, &fail2);
break;
case 2:
/* Compute AX and store in Y */
nag_superlu_matrix_product(order, notran, n, 2, one, icolzp, irowix, a, x,
ldx, zero, y, ldy, &fail2);
case 3:
/* Compute A^T Y and store in X */
nag_superlu_matrix_product(order, tran, n, 2, one, icolzp, irowix, a, y,
ldy, zero, x, ldx, &fail2);
case 4:
/* Compute AZ and store in P */
nag_superlu_matrix_product(order, notran, n, 1, one, icolzp, irowix, a, z,
n, zero, p, n, &fail2);
case 5:
/* Compute A^T Z and store in R */
nag_superlu_matrix_product(order, tran, n, 1, one, icolzp, irowix, a, z,
n, zero, r, n, &fail2);
}
if (fail2.code != NE_NOERROR) {
printf("Error from nag_superlu_matrix_product (f11mkc).\n%s\n",
fail2.message);
exit_status = 1;
goto END;
}
}
END_REVCM:

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

/* Print solution using
* nag_gen_real_mat_print (x04cac)
* Print real general matrix (easy-to-use)
*/
nag_gen_real_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_real_mat_print (x04cac)\n%s\n", fail.message);
exit_status = 3;
goto END;
}

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