/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:33:11 */ /*FOR_C Options SET: ftn=u io=c no=p op=aimnv pf=,p_djacg1 s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include #include #include #include "p_djacg1.h" #include /* PARAMETER translations */ #define LDFJAC M #define LIWK 21 #define LWK (3*M + 18) #define M 1 #define N 2 /* end of PARAMETER translations */ int main( ) { byte uv; long int _l0, i, iopt[4], iwk[LIWK], mode; double a, actual[N][LDFJAC], b, c, f[M], f2, fac[N], fjac[N][LDFJAC], maxerr, re[2][4], u, wk[LWK], xscale[N], y[N]; static char what[4][56]={"One sided partials, default settings ", "One sided partial, second partial known and skipped ","One sided partials, first partial accumulated ", "Central difference partials "}; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const F = &f[0] - 1; double *const Fac = &fac[0] - 1; long *const Iopt = &iopt[0] - 1; long *const Iwk = &iwk[0] - 1; double *const Wk = &wk[0] - 1; double *const Xscale = &xscale[0] - 1; double *const Y = &y[0] - 1; /* end of OFFSET VECTORS */ /*>> 2006-04-12 DRDJACG1 Hanson -- Reduced lengths of djacg work arrays. *>> 2003-07-08 DRDJACG1 Hanson -- Check for MODE < 0. *>> 2003-07-07 DRDJACG1 Krogh -- Changed 3 arg max to 2 arg for C conv. *>> 2002-06-21 DRDJACG1 R. J. Hanson Example 1 Code, with Download * Demo driver for DJACG, using numerical derivatives for a gradient. * ------------------------------------------------------------------ *-- D replaces "?": DR?JACG1, ?JACG */ /* The function used for this demo is f(y_1,y_2) = a*exp(b*y_1)+c*y_1 * *(y_2)**2. * Its gradient vector is (df/dy_1,df/dy_2)= * (a*b*exp(b*y_1)+c*(y_2)**2, 2*c*y_1*y_2). * This is used for comparison of the computed results with actual * results. */ /* This driver shows how to: */ /* A. Compute approximate derivatives using one-sided divided * differences * B. Use one-sided differences on the first component and * analytically compute the second component * C. Accumulate a known term of the first component with a * differenced term that is not known a priori * D. Use central differences on both components * (No checking for error conditions, i.e. MODE < 0.) * Define sizes and parameters. */ /* Define data and point of evaluation: */ a = 2.5e0; b = 3.4e0; c = 4.5e0; Y[1] = 2.1e0; Y[2] = 3.2e0; /* Machine precision, for measuring errors */ u = DBL_EPSILON; /* Set defaults for increments and scaling: */ Fac[1] = 0.e0; Xscale[1] = 0.e0; /* Compute true values of partials. */ actual[0][0] = a*b*exp( b*Y[1] ) + c*SQ(Y[2]); actual[1][0] = 2*c*Y[1]*Y[2]; /* A. No variable gets special treatment */ Iopt[1] = 0; /* Start each problem with MODE=0. Other starting * values are errors. Values < 0 or > N are caught. */ mode = 0; L_10: ; Wk[1] = a*exp( b*Y[1] ) + c*Y[1]*SQ(Y[2]); /* This sets the function value used in forming one-sided differences. */ if (mode == 0) { F[1] = Wk[1]; } djacg( &mode, M, N, y, f, (double*)fjac, LDFJAC, xscale, fac, iopt, wk, LWK, iwk, LIWK ); if (mode > 0) goto L_10; /* Check for an error condition. */ if (mode < 0) { printf(" Initial error in argument number %2ld\n", -mode); goto L_60; } /* Check the relative accuracy of one-sided differences. * They should be good to about half-precision. */ fjac[0][0] = (fjac[0][0] - actual[0][0])/actual[0][0]; fjac[1][0] = (fjac[1][0] - actual[1][0])/actual[1][0]; re[0][0] = fjac[0][0]/sqrt( u ); re[1][0] = fjac[1][0]/sqrt( u ); /* B. Skip variable number 2. */ Iopt[1] = 1; Iopt[2] = -4; Iopt[3] = 2; mode = 0; L_20: ; Wk[1] = a*exp( b*Y[1] ) + c*Y[1]*SQ(Y[2]); if (mode == 0) { F[1] = Wk[1]; /* The second component partial is skipped, * since it is known analytically */ fjac[1][0] = 2.e0*c*Y[1]*Y[2]; } djacg( &mode, M, N, y, f, (double*)fjac, LDFJAC, xscale, fac, iopt, wk, LWK, iwk, LIWK ); if (mode > 0) goto L_20; /* Check for an error condition. */ if (mode < 0) { printf(" Initial error in argument number %2ld\n", -mode); goto L_60; } fjac[0][0] = (fjac[0][0] - actual[0][0])/actual[0][0]; fjac[1][0] = (fjac[1][0] - actual[1][0])/actual[1][0]; re[0][1] = fjac[0][0]/sqrt( u ); re[1][1] = fjac[1][0]/sqrt( u ); /* C. Accumulate a part of the first partial. */ Iopt[1] = -3; Iopt[2] = 1; /* Shift to using one-sided differences for the * rest of the variables. */ Iopt[3] = -1; Iopt[4] = 0; mode = 0; L_30: ; /* Since part of the partial is known, evaluate what is * to be differenced. */ if (mode != 2) Wk[1] = a*exp( b*Y[1] ); if (mode == 0) { /* Start with part of the derivative that is known. */ F[1] = Wk[1]; fjac[0][0] = c*SQ(Y[2]); /* This is the function value for the partial wrt y_2. */ f2 = c*Y[1]*SQ(Y[2]); } if (mode == 2) { /* The function value for the second partial has the part removed * that depends on the first variable only. */ F[1] = f2; Wk[1] = c*Y[1]*SQ(Y[2]); } djacg( &mode, M, N, y, f, (double*)fjac, LDFJAC, xscale, fac, iopt, wk, LWK, iwk, LIWK ); if (mode > 0) goto L_30; /* Check for an error condition. */ if (mode < 0) { printf(" Initial error in argument number %2ld\n", -mode); goto L_60; } fjac[0][0] = (fjac[0][0] - actual[0][0])/actual[0][0]; fjac[1][0] = (fjac[1][0] - actual[1][0])/actual[1][0]; re[0][2] = fjac[0][0]/sqrt( u ); re[1][2] = fjac[1][0]/sqrt( u ); /* D. Use central differences and get more accuracy. * Twice the function evaluations are needed. */ Iopt[1] = -2; Iopt[2] = 0; /* Set the increment used at the default value. * This value must be assigned when using central differences. */ Wk[3*M + 3] = 0.e0; mode = 0; L_40: ; Wk[1] = a*exp( b*Y[1] ) + c*Y[1]*SQ(Y[2]); if (mode == 0) { F[1] = Wk[1]; } djacg( &mode, M, N, y, f, (double*)fjac, LDFJAC, xscale, fac, iopt, wk, LWK, iwk, LIWK ); if (mode > 0) goto L_40; /* Check for an error condition. */ if (mode < 0) { printf(" Initial error in argument number %2ld\n", -mode); goto L_60; } /* Check the relative accuracy of central differences. * They should be good to about two thirds-precision. */ fjac[0][0] = (fjac[0][0] - actual[0][0])/actual[0][0]; fjac[1][0] = (fjac[1][0] - actual[1][0])/actual[1][0]; f2 = pow(3.e0*u,2.e0/3.e0); re[0][3] = fjac[0][0]/f2; re[1][3] = fjac[1][0]/f2; /* Output the results and what is expected. */ printf("Rel Err of partials, f= a*exp(b*y_1)+c*y_1*(y_2)**2.\nCase df/dy_1 df/dy_2, u=macheps**.5, v=(3*macheps)**(2/3)\n"); maxerr = 0.e0; uv = 'u'; for (i = 1; i <= 4; i++) { maxerr = fmax( maxerr, fmax( fabs( re[0][i - 1] ), fabs( re[1][i - 1] ) ) ); if (i == 4) uv = 'v'; printf("%3ld%7.2f%c %7.2f%c %55.55s\n", i, re[0][i - 1], uv, re[1][i - 1], uv, what[i - 1]); } /* All expected relative errors (in units of truncation error) * should not exceed 8. If they do there may be an error. */ if (maxerr <= 8.e0) { printf(" Numbers above with absolute values .le. 8 are considered acceptable.\n"); } else { printf(" Numbers above with absolute values .gt. 8 are considered unacceptable.\n"); } L_60: ; exit(0); } /* end of function */