/* nag_opt_qp (e04nfc) Example Program. * * Copyright 1991 Numerical Algorithms Group. * * Mark 2, 1991. * Mark 6 revised, 2000. * Mark 7 revised, 2001. * Mark 8 revised, 2004. * */ #include #include #include #include #include #ifdef __cplusplus extern "C" { #endif static void qphess1(Integer n, Integer jthcol, double h[], Integer tdh, double x[], double hx[], Nag_Comm *comm); static void qphess2(Integer n, Integer jthcol, double h[], Integer tdh, double x[], double hx[], Nag_Comm *comm); static void qphess3(Integer n, Integer jthcol, double h[], Integer tdh, 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_qp (e04nfc) Example Program Results.\n"); exit_status_ex1 = ex1(); exit_status_ex2 = ex2(); return exit_status_ex1 == 0 && exit_status_ex2 == 0 ? 0 : 1; } #define A(I,J) a[(I)*tda + J] #define H(I,J) h[(I)*tdh + J] static int ex1(void) { Integer exit_status=0, i, j, n, nbnd, nclin, tda, tdh; NagError fail; double *a=0, bigbnd, *bl=0, *bu=0, *cvec=0, *h=0, objf, *x=0; INIT_FAIL(fail); Vprintf("\nExample 1: default options used.\n"); /* Define the problem. This example is due to Bunch and Kaufman, * `A computational method for the indefinite quadratic programming * problem ', Linear Algebra and its Applications, 34, 341-370 (1980). * * h = the QP Hessian matrix. * a = the general constraint matrix. * bl = the lower bounds on x and A*x. * bu = the upper bounds on x and A*x. * x = the initial estimate of the solution. * * Set the actual problem dimensions. * n = the number of variables. * nclin = the number of general linear constraints (may be 0). */ n = 8; nclin = 7; if (n>0 && nclin>=0) { nbnd = nclin + n; if ( !( a = NAG_ALLOC(nclin*n, double)) || !( h = NAG_ALLOC(n*n, double)) || !( x = NAG_ALLOC(n, double)) || !( cvec = NAG_ALLOC(n, double)) || !( bl = NAG_ALLOC(nbnd, double)) || !( bu = NAG_ALLOC(nbnd, double)) ) { Vprintf("Allocation failure\n"); exit_status = -1; goto END; } tda = n; tdh = n; } else { Vprintf("Invalid n or nclin.\n"); exit_status = 1; return exit_status; } /* Define the value used to denote ``infinite'' bounds. */ bigbnd = 1.0e20; for (i = 0; i < nclin; ++i) for (j = 0; j < n; ++j) A(i,j) = 0.0; for (i = 0; i < nclin; ++i) { A(i,i) = -1.0; A(i,i+1) = 1.0; bl[n + i] = -1.0 - 0.05*(double)i; bu[n + i] = bigbnd; } for (j = 0; j < n; ++j) { bl[j] = -(double)(j+1) - 0.1*(double)(j); bu[j] = (double)(j+1); cvec[j] = (double)(7 - j); } for (i = 0; i < n; ++i) { for (j = i+1; j < n; ++j) H(i,j) = (double)(ABS(i-j)); H(i,i) = 1.69; } /* Set the initial estimate of the solution. */ x[0] = -1.0; x[1] = -2.0; x[2] = -3.0; x[3] = -4.0; x[4] = -5.0; x[5] = -6.0; x[6] = -7.0; x[7] = -8.0; /* Solve the QP problem. */ /* nag_opt_qp (e04nfc). * Quadratic programming */ nag_opt_qp(n, nclin, a, tda, bl, bu, cvec, h, tdh, NULLFN, x, &objf, E04_DEFAULT, NAGCOMM_NULL, &fail); if (fail.code != NE_NOERROR) { Vprintf("Error from nag_opt_qp (e04nfc).\n%s\n", fail.message); exit_status = 1; goto END; } Vprintf("Re-solve problem with the Hessian defined by function qphess1.\n"); /* Set a new initial estimate of the solution. */ x[0] = -1.0; x[1] = 12.0; x[2] = -3.0; x[3] = 14.0; x[4] = -5.0; x[5] = 16.0; x[6] = -7.0; x[7] = 18.0; /* Solve the QP problem. */ /* nag_opt_qp (e04nfc), see above. */ nag_opt_qp(n, nclin, a, tda, bl, bu, cvec, (double *)0, tdh, qphess1, x, &objf, E04_DEFAULT, NAGCOMM_NULL, &fail); if (fail.code != NE_NOERROR) { Vprintf("Error from nag_opt_qp (e04nfc).\n%s\n", fail.message); exit_status = 1; goto END; } END: if (a) NAG_FREE(a); if (h) NAG_FREE(h); if (x) NAG_FREE(x); if (cvec) NAG_FREE(cvec); if (bl) NAG_FREE(bl); if (bu) NAG_FREE(bu); return exit_status; } /* ex1 */ static void qphess1(Integer n, Integer jthcol, double h[], Integer tdh, double x[], double hx[], Nag_Comm *comm) { /* qphess1 computes the vector Hx = (H)*x for some matrix H * that defines the Hessian of the required QP problem. * * In this version qphess the Hessian matrix is implicit. * The array h[] is not accessed. There is no special coding * for the case jthcol > 0 */ Integer i, j; double sum; for (i = 0; i < n; ++i) { sum = 1.69*x[i]; for (j = 0; j < n; ++j) sum += x[j]*(double)ABS(i - j); hx[i] = sum; } } /* qphess1 */ static int ex2(void) { Nag_Boolean print; Integer exit_status=0, i, j, n, nbnd, nclin, tda, tdh; NagError fail; Nag_E04_Opt options; double *a=0, *bl=0, *bu=0, *cvec=0, *h=0, objf, *x=0; INIT_FAIL(fail); Vprintf("\nExample 2: some optional parameters are set.\n"); Vscanf(" %*[^\n]"); /* Skip heading in data file */ /* Set the actual problem dimensions. * n = the number of variables. * nclin = the number of general linear constraints (may be 0). */ n = 7; nclin = 7; if (n>0 && nclin>=0) { nbnd = n + nclin; if ( !( x = NAG_ALLOC(n, double)) || !( cvec = NAG_ALLOC(n, double)) || !( a = NAG_ALLOC(nclin*n, double)) || !( h = NAG_ALLOC(n*n, double)) || !( bl = NAG_ALLOC(nbnd, double)) || !( bu = NAG_ALLOC(nbnd, double)) ) { Vprintf("Allocation failure\n"); exit_status = -1; goto END; } tda = n; tdh = n; } else { Vprintf("Invalid n or nclin.\n"); exit_status = 1; return exit_status; } /* cvec = the coefficients of the explicit linear term of f(x). * a = the linear constraint matrix. * bl = the lower bounds on x and A*x. * bu = the upper bounds on x and A*x. * x = the initial estimate of the solution. */ /* Read the coefficients of the explicit linear term of f(x). */ Vscanf(" %*[^\n]"); /* Skip heading in data file */ for (i = 0; i < n; ++i) Vscanf("%lf",&cvec[i]); /* Read the linear constraint matrix A. */ Vscanf(" %*[^\n]"); /* Skip heading in data file */ for (i = 0; i < nclin; ++i) for (j = 0; j < n; ++j) Vscanf("%lf",&A(i,j)); /* Read the bounds. */ nbnd = n + nclin; Vscanf(" %*[^\n]"); /* Skip heading in data file */ for (i = 0; i < nbnd; ++i) Vscanf("%lf", &bl[i]); Vscanf(" %*[^\n]"); /* Skip heading in data file */ 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]); /* nag_opt_init (e04xxc). * Initialization function for option setting */ nag_opt_init(&options); /* Initialise options structure */ /* Set one option directly * Bounds >= inf_bound will be treated as plus infinity. * Bounds <= -inf_bound will be treated as minus infinity. */ options.inf_bound = 1.0e21; /* Read remaining option values from file */ print = Nag_TRUE; /* nag_opt_read (e04xyc). * Read options from a text file */ nag_opt_read("e04nfc", "stdin", &options, print, "stdout", &fail); if (fail.code != NE_NOERROR) { Vprintf("Error from nag_opt_read (e04xyc).\n%s\n", fail.message); exit_status = 1; goto END; } /* Solve the problem from a cold start. * The Hessian is defined implicitly by function qphess2. */ /* nag_opt_qp (e04nfc), see above. */ nag_opt_qp(n, nclin, a, tda, bl, bu, cvec, (double *)0, tdh, qphess2, x, &objf, &options, NAGCOMM_NULL, &fail); if (fail.code != NE_NOERROR) { Vprintf("Error from nag_opt_qp (e04nfc).\n%s\n", fail.message); exit_status = 1; goto END; } /* The following is for illustrative purposes only. We do a warm * start with the final working set of the previous run. * This time we store the Hessian explicitly in h[][], and use * the corresponding function qphess3(). * Only the final solution from the results is printed. */ Vprintf("\nA run of the same example with a warm start:\n"); options.start = Nag_Warm; options.print_level = Nag_Soln; for (i = 0; i < n; ++i) { for (j = 0; j < n; ++j) H(i,j) = 0.0; if (i <= 4) H(i,i) = 2.0; else H(i,i) = -2.0; } H(2,3) = 2.0; H(3,2) = 2.0; H(5,6) = -2.0; H(6,5) = -2.0; /* Solve the problem again. */ /* nag_opt_qp (e04nfc), see above. */ nag_opt_qp(n, nclin, a, tda, bl, bu, cvec, h, tdh, qphess3, x, &objf, &options, NAGCOMM_NULL, &fail); if (fail.code != NE_NOERROR) { Vprintf("Error from nag_opt_qp (e04nfc).\n%s\n", fail.message); exit_status = 1; } /* Free memory allocated by nag_opt_qp (e04nfc) to pointers in options */ /* nag_opt_free (e04xzc). * Memory freeing function for use with option setting */ nag_opt_free(&options, "all", &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 (x) NAG_FREE(x); if (cvec) NAG_FREE(cvec); if (a) NAG_FREE(a); if (h) NAG_FREE(h); if (bl) NAG_FREE(bl); if (bu) NAG_FREE(bu); return exit_status; } /* ex2 */ static void qphess2(Integer n, Integer jthcol, double h[], Integer tdh, double x[], double hx[], Nag_Comm *comm) { /* In this version of qphess the Hessian matrix is implicit. * The array h[] is not accessed. There is no special coding * for the case jthcol > 0. */ hx[0] = 2.0*x[0]; hx[1] = 2.0*x[1]; hx[2] = 2.0*(x[2] + x[3]); hx[3] = hx[2]; hx[4] = 2.0*x[4]; hx[5] = -2.0*(x[5] + x[6]); hx[6] = hx[5]; } /* qphess2 */ #undef H static void qphess3(Integer n, Integer jthcol, double h[], Integer tdh, double x[], double hx[], Nag_Comm *comm) { /* In this version of QPHESS, the matrix H is stored in h[] * as a full two-dimensional array. */ #define H(I,J) h[(I)*tdh + (J)] Integer i, j; if (jthcol != 0) { /* Special case -- extract one column of H. */ j = jthcol - 1; for (i = 0; i < n; ++i) hx[i] = H(i,j); } else { /* Normal Case. */ for (i = 0; i < n; ++i) hx[i] = 0.0; for (i = 0; i < n; ++i) for (j = 0; j < n; ++j) hx[i] += H(i,j)*x[j]; } } /* qphess3 */