/*****************************************************************************/ /* obstacle.c */ /* */ /* Implementation of the obstacle problem as described in "MCPLIB: A */ /* Collection of Nonlinear Mixed Complementarity Problems" by Dirkse and */ /* Ferris in Optimization Methods and Software, 5, 1995, 319-345. */ /*****************************************************************************/ #include #include #include #include "Standalone.h" /*****************************************************************************/ /* We will use a constant force. */ /*****************************************************************************/ #define FORCE_CONSTANT 1.0 static int N; /* The number of elements on the X and Y axes*/ static int fill_structure; /* Do we need to fill in the structure of */ /* the Jacobian? */ int funcEval(int n, double *z, double *f) { /************************************************************************/ /* Evaluate the function at z, placing the result in f. */ /************************************************************************/ int i, j, l; double dt; /************************************************************************/ /* The number of variables is N*N, the interior vertices of the */ /* triangularization. */ /************************************************************************/ if (n != N*N) { fprintf (stderr, "Bad dimension argument 'n' to F: exiting . . .\n"); exit (-1); } dt = -(FORCE_CONSTANT) / ((N + 1.0)*(N + 1.0)); /************************************************************************/ /* Calculate the function value. */ /************************************************************************/ f[0] = dt + 4*z[0] - z[1] - z[N]; for (i = 1; i < N-1; i++) { f[i] = dt + 4*z[i] - z[i-1] - z[i+1] - z[i+N]; } f[N-1] = dt + 4*z[N-1] - z[N-2] - z[N-1+N]; for (j = 1; j < N-1; j++) { l = j*N; f[l] = dt + 4*z[l] - z[l+1] - z[l-N] - z[l+N]; for (i = 1; i < N-1; i++) { f[l+i] = dt + 4*z[l+i] - z[l+i-1] - z[l+i+1] - z[l+i-N] - z[l+i+N]; } f[l+N-1] = dt + 4*z[l+N-1] - z[l+N-2] - z[l-1] - z[l+N-1+N]; } l = (N-1)*N; f[l] = dt + 4*z[l] - z[l+1] - z[l-N]; for (i = 1; i < N-1; i++) { f[l+i] = dt + 4*z[l+i] - z[l+i-1] - z[l+i+1] - z[l+i-N]; } f[l+N-1] = dt + 4*z[l+N-1] - z[l+N-2] - z[l-1]; return 0; } int jacEval(int n, int nnz, double *z, int *col_start, int *col_len, int *row, double *data) { /************************************************************************/ /* Evaluate the jacobian at z, placing the result in col_start, col_len,*/ /* row, and data. */ /* */ /* Note: We only need to fill in the structure of the jacobian */ /* (col_start, col_len, row) once. We do this by looking at the */ /* fill_structure variable. */ /************************************************************************/ int i, j, k, l; /************************************************************************/ /* The number of variables is N*N, the interior vertices of the */ /* triangularization. */ /************************************************************************/ if (n != N*N) { fprintf (stderr, "Bad dimension argument 'n' to " "sparseJ: exiting . . .\n"); exit (-1); } /************************************************************************/ /* The number of nonzeros in the difference approximation is 5*N*N - 4*N*/ /* Exit if there is not enough room allocated. */ /************************************************************************/ if (nnz < 5*N*N - 4*N) { fprintf (stderr, "Bad dimension argument 'nnz' to " "sparseJ: exiting . . .\n"); exit (-1); } if (fill_structure) { /*******************************************************************/ /* This is the first call to jacEval, fill in the structure. */ /*******************************************************************/ k = 0; l = 0; col_start[k] = 1; row[l++] = 1; row[l++] = 2; row[l++] = N + 1; col_len[k] = l; k++; for (i = 1; i < N-1; i++) { col_start[k] = l + 1; row[l++] = k; row[l++] = k + 1; row[l++] = k + 2; row[l++] = k + N + 1; col_len[k] = l - col_start[k] + 1; k++; } col_start[k] = l + 1; row[l++] = k; row[l++] = k + 1; row[l++] = k + N + 1; col_len[k] = l - col_start[k] + 1; k++; for (j = 1; j < N-1; j++) { col_start[k] = l + 1; row[l++] = k - N + 1; row[l++] = k + 1; row[l++] = k + 2; row[l++] = k + N + 1; col_len[k] = l - col_start[k] + 1; k++; for (i = 1; i < N-1; i++) { col_start[k] = l + 1; row[l++] = k - N + 1; row[l++] = k; row[l++] = k + 1; row[l++] = k + 2; row[l++] = k + N + 1; col_len[k] = l - col_start[k] + 1; k++; } col_start[k] = l + 1; row[l++] = k - N + 1; row[l++] = k; row[l++] = k + 1; row[l++] = k + N + 1; col_len[k] = l - col_start[k] + 1; k++; } col_start[k] = l + 1; row[l++] = k - N + 1; row[l++] = k + 1; row[l++] = k + 2; col_len[k] = l - col_start[k] + 1; k++; for (i = 1; i < N-1; i++) { col_start[k] = l + 1; row[l++] = k - N + 1; row[l++] = k; row[l++] = k + 1; row[l++] = k + 2; col_len[k] = l - col_start[k] + 1; k++; } col_start[k] = l + 1; row[l++] = k - N + 1; row[l++] = k; row[l++] = k + 1; col_len[k] = l - col_start[k] + 1; k++; fill_structure = 0; } /************************************************************************/ /* Fill in the data values for the jacobian. */ /************************************************************************/ l = 0; data[l++] = 4; data[l++] = -1; data[l++] = -1; for (i = 1; i < N-1; i++) { data[l++] = -1; data[l++] = 4; data[l++] = -1; data[l++] = -1; } data[l++] = -1; data[l++] = 4; data[l++] = -1; for (j = 1; j < N-1; j++) { data[l++] = -1; data[l++] = 4; data[l++] = -1; data[l++] = -1; for (i = 1; i < N-1; i++) { data[l++] = -1; data[l++] = -1; data[l++] = 4; data[l++] = -1; data[l++] = -1; } data[l++] = -1; data[l++] = -1; data[l++] = 4; data[l++] = -1; } data[l++] = -1; data[l++] = 4; data[l++] = -1; for (i = 1; i < N-1; i++) { data[l++] = -1; data[l++] = -1; data[l++] = 4; data[l++] = -1; } data[l++] = -1; data[l++] = -1; data[l++] = 4; return 0; } int main() { double *lb; /* Lower bounds on the variables */ double *ub; /* Upper bounds on the variables */ double *z; /* Solution vector */ double *f; /* Function evaluation */ int n; /* Number of variable (N*N) */ int nnz; /* Number of nonzeros (5*N*N - 4*N) */ int status; /* Termination status from PATH */ int i, j; /************************************************************************/ /* Read in the number of elements to use on the X and Y axes. */ /************************************************************************/ printf("Input n: "); status = scanf("%d", &N); if ((status != 1) || (N < 2)) { N = 2; } /************************************************************************/ /* Determine n (number of variables) and nnz (number of nonzeros) */ /************************************************************************/ n = N*N; nnz = 5*N*N - 4*N; /************************************************************************/ /* Allocate space for the bounds, starting point, and function value. */ /************************************************************************/ lb = (double *)malloc(sizeof(double)*n); ub = (double *)malloc(sizeof(double)*n); z = (double *)malloc(sizeof(double)*n); f = (double *)malloc(sizeof(double)*n); /************************************************************************/ /* Fill in the lower and upper bounds and a starting point. */ /************************************************************************/ for (j = 0; j < N; j++) { for (i = 0; i < N; i++) { lb[i + j*N] = pow(sin(9.2*i/(N+1.0))*sin(9.3*j/(N+1.0)), 3.0); ub[i + j*N] = pow(sin(9.2*i/(N+1.0))*sin(9.3*j/(N+1.0)), 2.0) + 0.02; z[i + j*N] = 0.1; } } /************************************************************************/ /* Set fill_structure to true (we ALWAYS need to fill in the structure */ /* of the jacobian at least once per call to PATH). */ /************************************************************************/ fill_structure = 1; /************************************************************************/ /* Call PATH. */ /************************************************************************/ pathMain(n, nnz, &status, z, f, lb, ub); /************************************************************************/ /* Deallocate memory. */ /************************************************************************/ free(lb); free(ub); free(z); free(f); return 0; }