/* nag_opt_qp (e04nfc) Example Program * * Copyright 1991 Numerical Algorithms Group. * * Mark 2, 1991. * Mark 6 revised, 2000. * Mark 7 revised, 2001. * */ #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 void ex1(void); static void ex2(void); #define MAXN 10 #define MAXLIN 7 #define MAXBND MAXN+MAXLIN int main(void) { /* 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("e04nfc Example Program Results.\n"); ex1(); ex2(); return EXIT_SUCCESS; } static void ex1(void) { double a[MAXLIN][MAXN], h[MAXN][MAXN]; double x[MAXN], cvec[MAXN]; double bl[MAXBND], bu[MAXBND]; double bigbnd, objf; Integer i, j, n, nclin, tda, tdh; static NagError fail; Vprintf("\nExample 1: default options used.\n"); fail.print = TRUE; /* 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; tda = MAXN; tdh = MAXN; /* 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. */ e04nfc(n, nclin, &a[0][0], tda, bl, bu, cvec, &h[0][0], tdh, NULLFN, x, &objf, E04_DEFAULT, NAGCOMM_NULL, &fail); if (fail.code == NE_NOERROR) { 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. */ e04nfc(n, nclin, &a[0][0], tda, bl, bu, cvec, (double *)0, tdh, qphess1, x, &objf, E04_DEFAULT, NAGCOMM_NULL, &fail); } if (fail.code != NE_NOERROR) exit(EXIT_FAILURE); } /* 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 void ex2(void) { double x[MAXN], cvec[MAXN]; double a[MAXLIN][MAXN], h[MAXN][MAXN]; double bl[MAXBND], bu[MAXBND]; double objf; Integer tda, tdh; Integer i, j, n, nclin, nbnd; Boolean print; Nag_E04_Opt options; static NagError fail, fail2; Vprintf("\nExample 2: some optional parameters are set.\n"); Vscanf(" %*[^\n]"); /* Skip heading in data file */ fail.print = TRUE; fail2.print = TRUE; /* Set the actual problem dimensions. * n = the number of variables. * nclin = the number of general linear constraints (may be 0). */ tda = MAXN; tdh = MAXN; n = 7; nclin = 7; /* 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]); e04xxc(&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 */ fail.print = TRUE; print = TRUE; e04xyc("e04nfc", "stdin", &options, print, "stdout", &fail); /* Solve the problem from a cold start. * The Hessian is defined implicitly by function qphess2. */ if (fail.code == NE_NOERROR) e04nfc(n, nclin, &a[0][0], tda, bl, bu, cvec, (double *)0, tdh, qphess2, x, &objf, &options, NAGCOMM_NULL, &fail); if (fail.code == NE_NOERROR) { /* 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. */ e04nfc(n, nclin, &a[0][0], tda, bl, bu, cvec, &h[0][0], tdh, qphess3, x, &objf, &options, NAGCOMM_NULL, &fail); } /* Free memory allocated by e04nfc to pointers in options */ e04xzc(&options, "all", &fail2); if (fail.code != NE_NOERROR || fail2.code != NE_NOERROR) exit(EXIT_FAILURE); } /* 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 */ 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 */