```/* nag_opt_sparse_convex_qp (e04nkc) Example Program.
*
* Copyright 2014 Numerical Algorithms Group.
*
* Mark 5, 1998.
* Mark 6 revised, 2000.
* Mark 7 revised, 2001.
* Mark 8 revised, 2004.
*
*/

#include <nag.h>
#include <stdio.h>
#include <string.h>
#include <nag_stdlib.h>
#include <nage04.h>

#ifdef __cplusplus
extern "C" {
#endif
static void NAG_CALL qphess(Integer ncolh, const double x[], double hx[],
Nag_Comm *comm);
#ifdef __cplusplus
}
#endif

/* Declare a data structure for passing sparse Hessian data to qphess */
typedef struct
{
double  *hess;
Integer *khess;
Integer *hhess;
} HessianData;

#define NAMES(I, J) names[(I)*9+J]
#define MAXHESSNNZ MAXNNZ

int main(void)
{
HessianData hess_data;
Integer     exit_status = 0, *ha = 0, *hhess = 0, i, icol, iobj, j, jcol;
Integer     *ka = 0, *khess = 0, m, n, nbnd, ncolh, ninf, nnz, nnz_hess;
Nag_Comm    comm;
Nag_E04_Opt options;
char        **crnames = 0, *names = 0;
double      *a = 0, *bl = 0, *bu = 0, *hess = 0, obj, sinf, *x = 0;
NagError    fail;

INIT_FAIL(fail);

printf("nag_opt_sparse_convex_qp (e04nkc) Example Program Results\n");
fflush(stdout);

/* Skip heading in data file */
scanf(" %*[^\n]");
/* Read the problem dimensions */
scanf(" %*[^\n]");
scanf("%ld%ld", &n, &m);

/* Read nnz, iobj, ncolh */
scanf(" %*[^\n]");
scanf("%ld%ld%ld", &nnz, &iobj, &ncolh);

if (n >= 1 && m >= 1 && nnz >= 1 && nnz <= n*m)
{
nbnd = n + m;
if (!(a = NAG_ALLOC(nnz, double)) ||
!(bl = NAG_ALLOC(nbnd, double)) ||
!(bu = NAG_ALLOC(nbnd, double)) ||
!(x = NAG_ALLOC(nbnd, double)) ||
!(ha = NAG_ALLOC(nnz, Integer)) ||
!(ka = NAG_ALLOC(n+1, Integer)) ||
!(khess = NAG_ALLOC(n+1, Integer)) ||
!(crnames = NAG_ALLOC(nbnd, char *)) ||
!(names = NAG_ALLOC(nbnd*9, char))
)
{
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
}
else
{
printf("Invalid n or m or nnz.\n");
exit_status = 1;
return exit_status;
}

/* Read the matrix and set up ka */
jcol = 1;
ka[jcol-1] = 0;
scanf(" %*[^\n]");
for (i = 0; i < nnz; ++i)
{
/* a[i] stores the (ha[i], icol) element of matrix */
scanf("%lf%ld%ld", &a[i], &ha[i], &icol);

/* Check whether we have started a new column */
if (icol == jcol+1)
{
ka[icol-1] = i;  /* Start of icol-th column in a */
jcol = icol;
}
else if (icol > jcol+1)
{
/* Index in a of the start of the icol-th column
* equals i, but columns jcol+1, jcol+2, ...,
* icol-1 are empty.  Set the corresponding elements
* of ka to i.
*/
for (j = jcol+1; j < icol; ++j)
ka[j-1] = i;

ka[icol-1] = i;
jcol = icol;
}
}
ka[n] = nnz;

/* If the last columns are empty, set ka accordingly */
if (n > icol)
{
for (j = icol; j <= n - 1; ++j)
ka[j] = nnz;
}

/* Read the bounds */
nbnd = n+m;
scanf(" %*[^\n]"); /* Skip heading in data file */
for (i = 0; i < nbnd; ++i)
scanf("%lf", &bl[i]);
scanf(" %*[^\n]");
for (i = 0; i < nbnd; ++i)
scanf("%lf", &bu[i]);

/* Read the column and row names */
scanf(" %*[^\n]"); /* Skip heading in data file */
scanf(" %*[^']");
for (i = 0; i < nbnd; ++i)
{
scanf(" '%8c'", &NAMES(i, 0));
NAMES(i, 8) = '\0';
crnames[i] = &NAMES(i, 0);
}

/* Read the initial estimate of x */
scanf(" %*[^\n]"); /* Skip heading in data file */
for (i = 0; i < n; ++i)
scanf("%lf", &x[i]);

/* Read nnz_hess */
scanf(" %*[^\n]");
scanf("%ld", &nnz_hess);

if (!(hess = NAG_ALLOC(nnz_hess, double)) ||
!(hhess = NAG_ALLOC(nnz_hess, Integer)))
{
printf("Allocation failure\n");
exit_status = -1;
goto END;
}

/* Read the hessian matrix and set up khess */
jcol = 1;
khess[jcol-1] = 0;
scanf(" %*[^\n]");
for (i = 0; i < nnz_hess; ++i)
{
/* hess[i] stores the (hhess[i], icol) element of matrix */
scanf("%lf%ld%ld", &hess[i], &hhess[i], &icol);

/* Check whether we have started a new column */
if (icol == jcol+1)
{
khess[icol-1] = i;  /* Start of icol-th column in hess */
jcol = icol;
}
else if (icol > jcol+1)
{
/* Index in hess of the start of the icol-th column
* equals i, but columns jcol+1, jcol+2, ...,
* icol-1 are empty.  Set the corresponding elements
* of khess to i.
*/
for (j = jcol+1; j < icol; ++j)
khess[j-1] = i;

khess[icol-1] = i;
}
}
khess[ncolh] = nnz_hess;

/* Initialize options structure */
/* nag_opt_init (e04xxc).
* Initialization function for option setting
*/
nag_opt_init(&options);
options.crnames = crnames;

/* Package up Hessian data for communication via comm */
hess_data.hess = hess;
hess_data.khess = khess;
hess_data.hhess = hhess;

comm.p = (Pointer)&hess_data;

/* Solve the problem */
/* nag_opt_sparse_convex_qp (e04nkc), see above. */
nag_opt_sparse_convex_qp(n, m, nnz, iobj, ncolh, qphess, a, ha, ka, bl, bu,
x, &ninf, &sinf, &obj, &options, &comm, &fail);
if (fail.code != NE_NOERROR)
{
printf("Error from nag_opt_sparse_convex_qp (e04nkc).\n%s\n",
fail.message);
exit_status = 1;
goto END;
}

printf("\nPerturb the problem and re-solve with warm start.\n");
fflush(stdout);
for (i = 0; i < nnz_hess; ++i)
hess[i] *= 1.001;

options.start = Nag_Warm;
options.print_level = Nag_Soln;
/* nag_opt_sparse_convex_qp (e04nkc), see above. */
nag_opt_sparse_convex_qp(n, m, nnz, iobj, ncolh, qphess, a, ha, ka, bl, bu,
x, &ninf, &sinf, &obj, &options, &comm, &fail);
if (fail.code != NE_NOERROR)
{
printf("Error from nag_opt_sparse_convex_qp (e04nkc).\n%s\n",
fail.message);
exit_status = 1;
goto END;
}

/* Free memory NAG-allocated to members of options */
/* nag_opt_free (e04xzc).
* Memory freeing function for use with option setting
*/
nag_opt_free(&options, "", &fail);
if (fail.code != NE_NOERROR)
{
printf("Error from nag_opt_free (e04xzc).\n%s\n", fail.message);
exit_status = 1;
goto END;
}

END:
NAG_FREE(a);
NAG_FREE(bl);
NAG_FREE(bu);
NAG_FREE(x);
NAG_FREE(hess);
NAG_FREE(ha);
NAG_FREE(ka);
NAG_FREE(hhess);
NAG_FREE(khess);
NAG_FREE(crnames);
NAG_FREE(names);

return exit_status;
}

static void NAG_CALL qphess(Integer ncolh, const double x[], double hx[],
Nag_Comm *comm)
{
Integer     i, j, jrow;
HessianData *hd = (HessianData *)(comm->p);
double      *hess = hd->hess;
Integer     *hhess = hd->hhess;
Integer     *khess = hd->khess;

for (i = 0; i < ncolh; ++i)
hx[i] = 0.0;

for (i = 0; i < ncolh; ++i)
{
/* For each column of Hessian */
for (j = khess[i]; j < khess[i+1]; ++j)
{
/* For each non-zero in column */

jrow = hhess[j] - 1;

/* Using symmetry of hessian, add contribution
* to hx of hess[j] as a diagonal or upper
* triangular element of hessian.
*/
hx[i] += hess[j]*x[jrow];

/* If hess[j] not a diagonal element add its
* contribution to hx as a lower triangular
* element of hessian.
*/
if (jrow != i)
hx[jrow] += hess[j]*x[i];
}
}
}  /* qphess */
```