/* nag_opt_sparse_convex_qp (e04nkc) Example Program. * * Copyright 1998 Numerical Algorithms Group. * * Mark 5, 1998. * Mark 6 revised, 2000. * Mark 7 revised, 2001. * Mark 8 revised, 2004. * */ #include #include #include #include #ifdef __cplusplus extern "C" { #endif static void qphess1(Integer ncolh, double x[], double hx[], Nag_Comm *comm); static void qphess2(Integer ncolh, double x[], double hx[], Nag_Comm *comm); #ifdef __cplusplus } #endif static int ex1(void); static int ex2(void); int main(void) { Integer exit_status_ex1=0; Integer exit_status_ex2=0; /* Two examples are called: ex1() uses the * default settings to solve a problem while * ex2() solves another problem with some * of the optional parameters set by the user. */ Vprintf("nag_opt_sparse_convex_qp (e04nkc) Example Program Results.\n"); exit_status_ex1 = ex1(); exit_status_ex2 = ex2(); return exit_status_ex1 == 0 && exit_status_ex2 == 0 ? 0 : 1; } static int ex1(void) { #define MAXNNZ 50 Integer exit_status=0, *ha=0, i, icol, iobj, j, jcol, *ka=0, m, n, nbnd; Integer ncolh, ninf, nnz; NagError fail; double *a=0, *bl=0, *bu=0, obj, sinf, *x=0; INIT_FAIL(fail); Vprintf("\nExample 1: default options used.\n"); Vscanf(" %*[^\n]"); /* Skip headings in data file */ Vscanf(" %*[^\n]"); /* Read the problem dimensions */ Vscanf(" %*[^\n]"); Vscanf("%ld%ld", &n, &m); /* Read nnz, iobj, ncolh */ Vscanf(" %*[^\n]"); Vscanf("%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(MAXNNZ, Integer)) || !( ka = NAG_ALLOC(n+1, Integer)) ) { Vprintf("Allocation failure\n"); exit_status = -1; goto END; } } else { Vprintf("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; Vscanf(" %*[^\n]"); for (i = 0; i < nnz; ++i) { /* a[i] stores the (ha[i], icol) element of matrix */ Vscanf("%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; Vscanf(" %*[^\n]"); /* Skip heading in data file */ for (i = 0; i < nbnd; ++i) Vscanf("%lf", &bl[i]); Vscanf(" %*[^\n]"); for (i = 0; i < nbnd; ++i) Vscanf("%lf", &bu[i]); /* Read the initial estimate of x */ Vscanf(" %*[^\n]"); /* Skip heading in data file */ for (i = 0; i < n; ++i) Vscanf("%lf", &x[i]); /* Solve the problem */ /* nag_opt_sparse_convex_qp (e04nkc). * Solves sparse linear programming or convex quadratic * programming problems */ nag_opt_sparse_convex_qp(n, m, nnz, iobj, ncolh, qphess1, a, ha, ka, bl, bu, x, &ninf, &sinf, &obj, E04_DEFAULT, NAGCOMM_NULL, &fail); if (fail.code != NE_NOERROR) { Vprintf("Error from nag_opt_sparse_convex_qp (e04nkc).\n%s\n", fail.message); exit_status = 1; goto END; } END: if (a) NAG_FREE(a); if (bl) NAG_FREE(bl); if (bu) NAG_FREE(bu); if (x) NAG_FREE(x); if (ha) NAG_FREE(ha); if (ka) NAG_FREE(ka); return exit_status; } /* ex1 */ static void qphess1(Integer ncolh, double x[], double hx[], Nag_Comm *comm) { /* Function to compute H*x. */ hx[0] = 2.0*x[0]; /* 2.0 0.0 0.0 0.0 0.0 0.0 0.0 */ hx[1] = 2.0*x[1]; /* 0.0 2.0 0.0 0.0 0.0 0.0 0.0 */ hx[2] = 2.0*(x[2] + x[3]); /* 0.0 0.0 2.0 2.0 0.0 0.0 0.0 */ hx[3] = hx[2]; /* 0.0 0.0 2.0 2.0 0.0 0.0 0.0 */ hx[4] = 2.0*x[4]; /* 0.0 0.0 0.0 0.0 2.0 0.0 0.0 */ hx[5] = 2.0*(x[5] + x[6]); /* 0.0 0.0 0.0 0.0 0.0 2.0 2.0 */ hx[6] = hx[5]; /* 0.0 0.0 0.0 0.0 0.0 2.0 2.0 */ } /* qphess1 */ /* Example 2 */ /* Declare a data structure for passing sparse Hessian data to qphess2 */ typedef struct { double *hess; Integer *khess; Integer *hhess; } HessianData; static int ex2(void) { #define NAMES(I,J) names[(I)*9+J] #define MAXHESSNNZ MAXNNZ HessianData hess_data; Integer exit_status=0, *ha=0, *hhess=0, i, icol, iobj, j, jcol, *ka=0; Integer *khess=0, m, n, nbnd, ncolh, ninf, nnz, nnz_hess; NagError fail; 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; INIT_FAIL(fail); Vprintf("\nExample 2: some optional parameters are set.\n"); Vscanf(" %*[^\n]"); /* Read the problem dimensions */ Vscanf(" %*[^\n]"); Vscanf("%ld%ld", &n, &m); /* Read nnz, iobj, ncolh */ Vscanf(" %*[^\n]"); Vscanf("%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)) || !( hess = NAG_ALLOC(nnz, double)) || !( ha = NAG_ALLOC(nnz, Integer)) || !( ka = NAG_ALLOC(n+1, Integer)) || !( hhess = NAG_ALLOC(nnz, Integer)) || !( khess = NAG_ALLOC(n+1, Integer)) || !(crnames = NAG_ALLOC(nbnd, char *)) || !(names = NAG_ALLOC(nbnd*9, char )) ) { Vprintf("Allocation failure\n"); exit_status = -1; goto END; } } else { Vprintf("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; Vscanf(" %*[^\n]"); for (i = 0; i < nnz; ++i) { /* a[i] stores the (ha[i], icol) element of matrix */ Vscanf("%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; Vscanf(" %*[^\n]"); /* Skip heading in data file */ for (i = 0; i < nbnd; ++i) Vscanf("%lf", &bl[i]); Vscanf(" %*[^\n]"); for (i = 0; i < nbnd; ++i) Vscanf("%lf", &bu[i]); /* Read the column and row names */ Vscanf(" %*[^\n]"); /* Skip heading in data file */ Vscanf(" %*[^']"); for (i = 0; i < nbnd; ++i) { Vscanf(" '%8c'", &NAMES(i,0)); NAMES(i,8) = '\0'; crnames[i] = &NAMES(i,0); } /* Read the initial estimate of x */ Vscanf(" %*[^\n]"); /* Skip heading in data file */ for (i = 0; i < n; ++i) Vscanf("%lf", &x[i]); /* Read nnz_hess */ Vscanf(" %*[^\n]"); Vscanf("%ld", &nnz_hess); /* Read the hessian matrix and set up khess */ jcol = 1; khess[jcol-1] = 0; Vscanf(" %*[^\n]"); for (i = 0; i < nnz_hess; ++i) { /* hess[i] stores the (hhess[i], icol) element of matrix */ Vscanf("%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, qphess2, a, ha, ka, bl, bu, x, &ninf, &sinf, &obj, &options, &comm, &fail); if (fail.code != NE_NOERROR) { Vprintf("Error from nag_opt_sparse_convex_qp (e04nkc).\n%s\n", fail.message); exit_status = 1; goto END; } Vprintf("\nPerturb the problem and re-solve with warm start.\n"); 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, qphess2, a, ha, ka, bl, bu, x, &ninf, &sinf, &obj, &options, &comm, &fail); if (fail.code != NE_NOERROR) { Vprintf("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) { Vprintf("Error from nag_opt_free (e04xzc).\n%s\n", fail.message); exit_status = 1; goto END; } END: if (a) NAG_FREE(a); if (bl) NAG_FREE(bl); if (bu) NAG_FREE(bu); if (x) NAG_FREE(x); if (hess) NAG_FREE(hess); if (ha) NAG_FREE(ha); if (ka) NAG_FREE(ka); if (hhess) NAG_FREE(hhess); if (khess) NAG_FREE(khess); if (crnames) NAG_FREE(crnames); if (names) NAG_FREE(names); return exit_status; } /* ex2 */ static void qphess2(Integer ncolh, 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]; } } } /* qphess2 */