/*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_djacg2 s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include #include #include "p_djacg2.h" #include "drdjacg2.h" /* PARAMETER translations */ #define LIW (30 + NEQ) #define LRW (45 + (5 + 2*MAXCON + 4)*NEQ + SQ(NEQ)) #define MAXCON 4 #define NDIG 11 #define NEQ 5 #define TOL (powi(10.e0,-NDIG)) /* end of PARAMETER translations */ /* COMMON translations */ struct t_counts { long int kr, kf, kc; } counts; /* end of COMMON translations */ int main( ) { LOGICAL32 constraint; long int _n, i, idid, info[16], iwork[LIW]; double atol[5], drl, drv, length, rtol[5], rwork[LRW], t, tout, y[5], yprime[5]; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const Atol = &atol[0] - 1; long *const Info = &info[0] - 1; long *const Iwork = &iwork[0] - 1; double *const Rtol = &rtol[0] - 1; double *const Rwork = &rwork[0] - 1; double *const Y = &y[0] - 1; double *const Yprime = &yprime[0] - 1; /* end of OFFSET VECTORS */ /*>> 2009-10-28 DRDJACG2 Krogh -- Remove implicit none. *>> 2009-10-26 DRDJACG2 Krogh -- *'s used for dims for NAG compiler. *>> 2008-11-01 DRDJACG2 Hanson -- Added row dim parameters to evaluators *>> 2006-12-19 DRDJACG2 Hanson -- Added tot. energy constraint in 3rd ex *>> 2006-04-11 DRDJACG1 Hanson -- Reduced lengths of djacg work arrays. *>> 2006-04-09 DRDJACG2 Krogh Fixed dimension of YSCALE. *>> 2003-07-08 DRDJACG2 R. J. Hanson Fix accum. option with djacg(). *>> 2002-02-01 DRDJACG2 R. J. Hanson Example 2 Code, with Download */ /* Solve a planar pendulum problem in rectangular coordinates. * The equation is transformed from "Index 3" to "Index 1" * by differentiating the length constraint twice. The system * is integrated (using DDASSFC) using these two constraints. * In the second integration the problem is transformed from * "Index 3" to "Index 0" by differentiating the length constraint * three times. The routine DDASSFD uses these three constraints. * A total energy constraint is added in routine DDASSFE. * This example shows that the constraints remain satisfied, * and the integration can succeed with either approach. * It is more efficient to use the "Index 0" problem than * the "Index 1" problem. */ /* (THIS VERSION USES NUMERICAL DIFFERENTIATION FOR THE PARTIALS * NEEDED IN ?DASLX, EXAMPLE 5, AVAILABLE WITH DOWNLOAD.) *--D replaces "?": DR?JACG2,?DASLX,?DASSFC,?DASSFD,?DASSFE,?COPY,?JACG */ /* Set number of equations: */ /* Set number of constraints. */ /* Work space sizes: */ /*++S Default NDIG = 4 *++ Default NDIG = 11 *++ Substitute for NDIG below */ /* Tolerances: */ for (i = 1; i <= NEQ; i++) { Atol[i] = TOL; Rtol[i] = TOL; } /* Setup options: */ for (i = 1; i <= 16; i++) { Info[i] = 0; } /* Use partial derivatives provided in evaluation routine: */ Info[5] = 2; /* Constrain solution, with 2 constraints: */ Info[10] = 2; /* Allow for more 10 * steps (than default): */ Info[12] = 50000; /* This is the pendulum length. */ length = 1.1e0; for (i = 1; i <= 1000; i += 10) { /* Integrate from T=I-1 to TOUT=T+1. Final TOUT=10. * When the solution first drifts away from the constraints, stop. */ t = (double)( i - 1 ); tout = t + 10.e0; ddaslx( ddassfc, NEQ, &t, y, yprime, tout, info, rtol, atol, &idid, rwork, LRW, iwork, LIW ); /* Compute residuals on length and its variation. They * should be smaller than the tolerances used. */ drl = (SQ(Y[1]) + SQ(Y[2]) - SQ(length))/(Rtol[1]*SQ(length) + Atol[1]); drv = (Y[1]*Y[3] + Y[2]*Y[4])/(Rtol[1]*(fabs( Y[1]*Y[3] ) + fabs( Y[1]*Y[4] )) + Atol[1]); constraint = fabs( drl ) <= 1.e0 && fabs( drv ) <= 1.e0; if (!constraint) goto L_40; } L_40: ; printf(" Example Results for a Constrained Pendulum Problem, Index 1\n"); printf(" (Numerical Partials)\n"); printf("\n T y_1 y_2 y_3 y_4 y_5\n%14.4e", tout); for(_n=0L; _n < sizeof(y)/sizeof(double); _n++) printf("%14.4e", y[_n]); printf("\n"); printf(" No. of Residual Evaluations Factorizations No. of User Solves\nConstraint Partials- %8ld %8ld %8ld\n\n", counts.kr, counts.kf, counts.kc); printf(" At the time value%10.2f the integration was stopped.\n", tout); if (!constraint) { printf("The pendulum length or variation has large weighted errors.\nThese should remain less than 1 is magnitude: %8.2f%8.2f\n\n", drl, drv); } else { printf(" The pendulum length and its variation have small errors.\n\n\n"); } /* Start the integration over for the index 0 problem. */ Info[1] = 0; /* Set the number of constraints to 3. */ Info[10] = 3; for (i = 1; i <= 1000; i += 10) { /* Integrate from T=I-1 to TOUT=T+1. */ t = (double)( i - 1 ); tout = t + 10.e0; ddaslx( ddassfd, NEQ, &t, y, yprime, tout, info, rtol, atol, &idid, rwork, LRW, iwork, LIW ); /* Compute residuals on length and its variation. They * should be smaller than the tolerances used. */ drl = (SQ(Y[1]) + SQ(Y[2]) - SQ(length))/(Rtol[1]*SQ(length) + Atol[1]); drv = (Y[1]*Y[3] + Y[2]*Y[4])/(Rtol[1]*(fabs( Y[1]*Y[3] ) + fabs( Y[1]*Y[4] )) + Atol[1]); constraint = fabs( drl ) <= 1.e0 && fabs( drv ) <= 1.e0; if (!constraint) goto L_60; } L_60: ; printf(" Example Results for a Constrained Pendulum Problem, Index 0\n"); printf(" (Numerical Partials)\n"); printf("\n T y_1 y_2 y_3 y_4 y_5\n%14.4e", tout); for(_n=0L; _n < sizeof(y)/sizeof(double); _n++) printf("%14.4e", y[_n]); printf("\n"); printf(" No. of Residual Evaluations Factorizations No. of User Solves\nConstraint Partials- %8ld %8ld %8ld\n\n", counts.kr, counts.kf, counts.kc); printf(" At the time value%10.2f the integration was stopped.\n", tout); if (!constraint) { printf("The pendulum length or variation has large weighted errors.\nThese should remain less than 1 is magnitude: %8.2f%8.2f\n\n", drl, drv); } else { printf(" The pendulum length and its variation have small errors.\n\n\n"); } /* Start the integration over for the index 0 problem. */ Info[1] = 0; /* Set the number of constraints to 4. * (This includes constant total energy). */ Info[10] = 4; for (i = 1; i <= 1000; i += 10) { /* Integrate from T=I-1 to TOUT=T+1. */ t = (double)( i - 1 ); tout = t + 10.e0; ddaslx( ddassfe, NEQ, &t, y, yprime, tout, info, rtol, atol, &idid, rwork, LRW, iwork, LIW ); /* Compute residuals on length and its variation. They * should be smaller than the tolerances used. */ drl = (SQ(Y[1]) + SQ(Y[2]) - SQ(length))/(Rtol[1]*SQ(length) + Atol[1]); drv = (Y[1]*Y[3] + Y[2]*Y[4])/(Rtol[1]*(fabs( Y[1]*Y[3] ) + fabs( Y[1]*Y[4] )) + Atol[1]); constraint = fabs( drl ) <= 1.e0 && fabs( drv ) <= 1.e0; if (!constraint) goto L_80; } L_80: ; printf(" Example Results for a Constrained Pendulum Problem, Index 0\n"); printf(" (Numerical Partials and Total Energy Constrained)\n"); printf("\n T y_1 y_2 y_3 y_4 y_5\n%14.4e", tout); for(_n=0L; _n < sizeof(y)/sizeof(double); _n++) printf("%14.4e", y[_n]); printf("\n"); printf(" No. of Residual Evaluations Factorizations No. of User Solves\nConstraint Partials- %8ld %8ld %8ld\n\n", counts.kr, counts.kf, counts.kc); printf(" At the time value%10.2f the integration was stopped.\n", tout); if (!constraint) { printf("The pendulum length or variation has large weighted errors.\nThese should remain less than 1 is magnitude: %8.2f%8.2f\n\n", drl, drv); } else { printf(" The pendulum length and its variation have small errors.\n\n\n"); } exit(0); } /* end of function */ /* PARAMETER translations */ #define GRAVITY 9.806650e0 #define LENGTH 1.1e0 #define MASS 1.e0 #define ZERO 0.e0 /* end of PARAMETER translations */ void /*FUNCTION*/ ddassfc( double *t, double y[], double yprime[], double delta[], double *d, long *ldd, double *cj, long *ires, double rwork[], long iwork[]) { #define D(I_,J_) (*(d+(I_)*(*ldd)+(J_))) long int iwk[21], mode; static long int iopt[5]; double f[5], wk[3*5 + 18]; static double facc[5], facj[5], lsq, mg, yscale[5]; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const Delta = &delta[0] - 1; double *const F = &f[0] - 1; double *const Facc = &facc[0] - 1; double *const Facj = &facj[0] - 1; long *const Iopt = &iopt[0] - 1; long *const Iwk = &iwk[0] - 1; long *const Iwork = &iwork[0] - 1; double *const Rwork = &rwork[0] - 1; double *const Wk = &wk[0] - 1; double *const Y = &y[0] - 1; double *const Yprime = &yprime[0] - 1; double *const Yscale = &yscale[0] - 1; /* end of OFFSET VECTORS */ /* Routine for the swinging simple pendulum problem, * with constraints on the index 2 and 3 equations. */ /* Use numerical derivatives but based on calls to DJACG. */ /* This is the setup call. YPRIME need not be solved for as we * are providing the correct initial values. */ if (*ires == 0) { mg = MASS*GRAVITY; lsq = SQ(LENGTH); Y[1] = LENGTH; Y[2] = ZERO; Y[3] = ZERO; Y[4] = ZERO; Y[5] = ZERO; Yprime[1] = ZERO; Yprime[2] = ZERO; Yprime[3] = ZERO; Yprime[4] = -GRAVITY; Yprime[5] = ZERO; counts.kr = 0; counts.kf = 0; counts.kc = 0; /* Initialize things needed for numerical differentiation. */ Facc[1] = ZERO; Facj[1] = ZERO; } /* The sytem residual value. */ if (*ires == 1) { Delta[1] = Y[3] - Yprime[1]; Delta[2] = Y[4] - Yprime[2]; Delta[3] = -Y[1]*Y[5] - MASS*Yprime[3]; Delta[4] = -Y[2]*Y[5] - MASS*Yprime[4] - mg; Delta[5] = MASS*(SQ(Y[3]) + SQ(Y[4])) - mg*Y[2] - lsq*Y[5]; /* Count residual evaluations. */ counts.kr += 1; } /* The partial derivative matrix. */ if (*ires == 2) { mode = 0; /* Accumulate all but one partial derivative columns. */ /* Set to accumulate variables 1-4: */ Iopt[1] = -3; Iopt[2] = 4; /* Shift to using one-sided differences. */ Iopt[3] = -1; /* Use on all remaining variables: number 5. */ Iopt[4] = 0; Yscale[1] = ZERO; /* Pre-compute partials wrt y'. These values are accumulated. */ D(0,0) = -*cj; D(1,1) = -*cj; D(2,2) = -MASS**cj; D(3,3) = D(2,2); L_10: ; /* The value of f(y) is normally returned in WK(1:5). * Note that all but one partial with respect to y' have been * analytically pre-computed and so the results are accumulated. */ Wk[1] = Y[3]; Wk[2] = Y[4]; Wk[3] = -Y[1]*Y[5]; Wk[4] = -Y[2]*Y[5]; Wk[5] = MASS*(SQ(Y[3]) + SQ(Y[4])) - mg*Y[2] - lsq*Y[5]; /* Compute f(y) and make a special copy step for the * evaluation point. */ if (mode == 0) { F[1] = Wk[1]; F[2] = Wk[2]; F[3] = Wk[3]; F[4] = Wk[4]; F[5] = Wk[5]; Facj[1] = ZERO; } /* This is a Math a la Carte code for numerical derivatives. */ djacg( &mode, 5, 5, y, f, d, *ldd, yscale, facj, iopt, wk, 33, iwk, 21 ); if (mode > 0) goto L_10; if (mode < 0) { printf(" DDASSFC: Initial error in argument number %2ld\n", -mode); puts( "1" ); exit(0); } /* All system partials are computed. Count an evaluation. */ counts.kf += 1; } /* The constraining equations after the corrector has converged. * Both partials and residuals are required. */ if (*ires == 5) { mode = 0; /* In this use of numerical differentiation, no special treatment * is required for the variables. The rest of IOPT(:) stays intact. */ Iopt[1] = 0; Yscale[1] = ZERO; L_20: ; /* The value of f(y) is normally returned in WK(1:2). */ Wk[1] = SQ(Y[1]) + SQ(Y[2]) - lsq; Wk[2] = Y[1]*Y[3] + Y[2]*Y[4]; /* Compute f(y) and make a special copy step for the point in question. */ if (mode == 0) { Delta[6] = Wk[1]; Delta[7] = Wk[2]; Facc[1] = ZERO; } /* This is a Math a la Carte code for numerical derivatives. */ djacg( &mode, 2, 5, y, &Delta[6], &D(0,5), *ldd, yscale, facc, iopt, wk, 33, iwk, 21 ); if (mode > 0) goto L_20; if (mode < 0) { printf(" DDASSFC: Initial error in argument number %2ld\n", -mode); puts( "2" ); exit(0); } /* All constraint partials are computed. Count an evaluation. */ counts.kc += 1; } return; #undef D } /* end of function */ /* PARAMETER translations */ #define THREE 3.e0 /* end of PARAMETER translations */ void /*FUNCTION*/ ddassfd( double *t, double y[], double yprime[], double delta[], double *p, long *ldp, double *cj, long *ires, double rwork[], long iwork[]) { #define P(I_,J_) (*(p+(I_)*(*ldp)+(J_))) long int iwk[21], mode; static long int iopt[2]; double f[5], wk[3*5 + 18]; static double facc[5], facj[5], lsq, mg, yscale[5]; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const Delta = &delta[0] - 1; double *const F = &f[0] - 1; double *const Facc = &facc[0] - 1; double *const Facj = &facj[0] - 1; long *const Iopt = &iopt[0] - 1; long *const Iwk = &iwk[0] - 1; long *const Iwork = &iwork[0] - 1; double *const Rwork = &rwork[0] - 1; double *const Wk = &wk[0] - 1; double *const Y = &y[0] - 1; double *const Yprime = &yprime[0] - 1; double *const Yscale = &yscale[0] - 1; /* end of OFFSET VECTORS */ /* Routine for the swinging simple pendulum problem, * with constraints on the index 1,2 amd 3 equations. * Use P(:,:) to distinguish from D(:,:) in routine DDASSFC. */ /* Use numerical derivatives but based on calls to DJACG. */ /* This is the setup call. */ if (*ires == 0) { mg = MASS*GRAVITY; lsq = SQ(LENGTH); Y[1] = LENGTH; Y[2] = ZERO; Y[3] = ZERO; Y[4] = ZERO; Y[5] = ZERO; Yprime[1] = ZERO; dcopy( 5, yprime, 0, yprime, 1 ); Yprime[4] = -GRAVITY; counts.kr = 0; counts.kf = 0; counts.kc = 0; Yscale[1] = ZERO; } /* The sytem residual value. */ if (*ires == 1) { Delta[1] = Y[3] - Yprime[1]; Delta[2] = Y[4] - Yprime[2]; Delta[3] = -Y[1]*Y[5] - MASS*Yprime[3]; Delta[4] = -Y[2]*Y[5] - MASS*Yprime[4] - mg; Delta[5] = -THREE*mg*Y[4] - lsq*Yprime[5]; counts.kr += 1; } /* The partial derivative matrix. */ if (*ires == 2) { mode = 0; Iopt[1] = 0; /* This is the last column number to be accumulated when * doing numerical differentiation. */ Iopt[2] = 5; L_10: ; /* The value of f(y) is normally returned in WK(1:5). * Note that the partials with respect to y' have been * analytically pre-computed and so the results are accumulated. */ Wk[1] = Y[3]; Wk[2] = Y[4]; Wk[3] = -Y[1]*Y[5]; Wk[4] = -Y[2]*Y[5]; Wk[5] = -THREE*mg*Y[4]; /* Compute f(y) and make a special copy step for the point in question. */ if (mode == 0) { F[1] = Wk[1]; F[2] = Wk[2]; F[3] = Wk[3]; F[4] = Wk[4]; F[5] = Wk[5]; Facj[1] = ZERO; } /* This is a Math a la Carte code for numerical derivatives. */ djacg( &mode, 5, 5, y, f, p, *ldp, yscale, facj, iopt, wk, 33, iwk, 21 ); if (mode > 0) goto L_10; if (mode < 0) { printf(" DDASSFD: Initial error in argument number %2ld\n", -mode); puts( "1" ); exit(0); } /* All system partials are computed. Count an evaluation. */ counts.kf += 1; /* This part of the derivative matrix is a diagonal, hence * easy to compute. */ P(0,0) += -*cj; P(1,1) += -*cj; P(2,2) += -MASS**cj; P(3,3) = P(2,2); P(4,4) += -lsq**cj; } /* The constraining equations after the corrector has converged. * Both partials and residuals are required. */ if (*ires == 5) { mode = 0; Iopt[1] = 0; Yscale[1] = ZERO; L_20: ; /* The value of f(y) is normally returned in WK(1:3). */ Wk[1] = SQ(Y[1]) + SQ(Y[2]) - lsq; Wk[2] = Y[1]*Y[3] + Y[2]*Y[4]; Wk[3] = MASS*(SQ(Y[3]) + SQ(Y[4])) - mg*Y[2] - lsq*Y[5]; /* Compute f(y) and make a special copy step for the * evaluation point. */ if (mode == 0) { Delta[6] = Wk[1]; Delta[7] = Wk[2]; Delta[8] = Wk[3]; Facc[1] = ZERO; } /* This is a Math a la Carte code for numerical derivatives. */ djacg( &mode, 3, 5, y, &Delta[6], &P(0,5), *ldp, yscale, facc, iopt, wk, 33, iwk, 21 ); if (mode > 0) goto L_20; if (mode < 0) { printf(" DDASSFD: Initial error in argument number %2ld\n", -mode); puts( "2" ); exit(0); } /* All constraint partials are computed. Count an evaluation. */ counts.kc += 1; } return; #undef P } /* end of function */ /* PARAMETER translations */ #define HALF 0.5e0 #define LSQ (SQ(LENGTH)) #define MG (MASS*GRAVITY) /* end of PARAMETER translations */ void /*FUNCTION*/ ddassfe( double *t, double y[], double yprime[], double delta[], double *p, long *ldp, double *cj, long *ires, double rwork[], long iwork[]) { #define P(I_,J_) (*(p+(I_)*(*ldp)+(J_))) long int iopt[2], iwk[21], mode; double f[5], facc[5], facj[5], wk[43], yscale[1]; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const Delta = &delta[0] - 1; double *const F = &f[0] - 1; double *const Facc = &facc[0] - 1; double *const Facj = &facj[0] - 1; long *const Iopt = &iopt[0] - 1; long *const Iwk = &iwk[0] - 1; long *const Iwork = &iwork[0] - 1; double *const Rwork = &rwork[0] - 1; double *const Wk = &wk[0] - 1; double *const Y = &y[0] - 1; double *const Yprime = &yprime[0] - 1; double *const Yscale = &yscale[0] - 1; /* end of OFFSET VECTORS */ /* Routine for the swinging simple pendulum problem, * with constraints on the index 1,2 amd 3 equations. * Also the constant energy constraint is added. * Use P(:,:) to distinguish from D(:,:) in routine DDASSFC. */ /* Use numerical derivatives but based on calls to DJACG. */ /* SAVE FACC, FACJ, YSCALE, IOPT */ /* This is the setup call. */ if (*ires == 0) { Y[1] = LENGTH; Y[2] = ZERO; Y[3] = ZERO; Y[4] = ZERO; Y[5] = ZERO; Yprime[1] = ZERO; Yprime[2] = ZERO; Yprime[3] = ZERO; Yprime[5] = ZERO; Yprime[4] = -GRAVITY; counts.kr = 0; counts.kf = 0; counts.kc = 0; Yscale[1] = ZERO; } /* The sytem residual value. */ if (*ires == 1) { Delta[1] = Y[3] - Yprime[1]; Delta[2] = Y[4] - Yprime[2]; Delta[3] = -Y[1]*Y[5] - MASS*Yprime[3]; Delta[4] = -Y[2]*Y[5] - MASS*Yprime[4] - MG; Delta[5] = -THREE*MG*Y[4] - LSQ*Yprime[5]; counts.kr += 1; } /* The partial derivative matrix. */ if (*ires == 2) { mode = 0; /* Accumulate all partial derivative columns. * Pre-compute partials wrt y' and accumulate results. */ Iopt[1] = -3; /* This is the last column number to be accumulated when * doing numerical differentiation. */ Iopt[2] = 5; /* This part of the derivative matrix is a diagonal, hence * easy to compute. */ P(0,0) = -*cj; P(1,1) = -*cj; P(2,2) = -MASS**cj; P(3,3) = P(2,2); P(4,4) = -LSQ**cj; L_10: ; /* The value of f(y) is normally returned in WK(1:5). * Note that the partials with respect to y' have been * analytically pre-computed and the results are accumulated. */ Wk[1] = Y[3]; Wk[2] = Y[4]; Wk[3] = -Y[1]*Y[5]; Wk[4] = -Y[2]*Y[5]; Wk[5] = -THREE*MG*Y[4]; /* Compute f(y) and make a special copy step for the point in question. */ if (mode == 0) { F[1] = Wk[1]; F[2] = Wk[2]; F[3] = Wk[3]; F[4] = Wk[4]; F[5] = Wk[5]; Facj[1] = ZERO; Yscale[1] = ZERO; } /* This is a Math a la Carte code for numerical derivatives. */ djacg( &mode, 5, 5, y, f, p, *ldp, yscale, facj, iopt, wk, 43, iwk, 21 ); if (mode > 0) goto L_10; if (mode < 0) { printf(" DDASSFD: Initial error in argument number %2ld\n", -mode); puts( "1" ); exit(0); } /* All system partials are computed. Count an evaluation. */ counts.kf += 1; } /* The constraining equations after the corrector has converged. * Both partials and residuals are required. */ if (*ires == 5) { mode = 0; Iopt[1] = 0; Facc[1] = ZERO; Yscale[1] = ZERO; L_20: ; /* The value of f(y) is normally returned in WK(1:3). */ Wk[1] = SQ(Y[1]) + SQ(Y[2]) - LSQ; Wk[2] = Y[1]*Y[3] + Y[2]*Y[4]; Wk[3] = MASS*(SQ(Y[3]) + SQ(Y[4])) - MG*Y[2] - LSQ*Y[5]; /* This constraint is constant total energy. */ Wk[4] = HALF*MASS*(SQ(Y[3]) + SQ(Y[4])) + MG*Y[2]; /* Compute f(y) and make a special copy step for the * evaluation point. */ if (mode == 0) { Delta[6] = Wk[1]; Delta[7] = Wk[2]; Delta[8] = Wk[3]; Delta[9] = Wk[4]; F[1] = Wk[1]; F[2] = Wk[2]; F[3] = Wk[3]; F[4] = Wk[4]; } /* This is a Math a la Carte code for numerical derivatives. */ djacg( &mode, 4, 5, y, f, &P(0,5), *ldp, yscale, facc, iopt, wk, 43, iwk, 21 ); if (mode > 0) goto L_20; if (mode < 0) { printf(" DDASSFE: Initial error in argument number %2ld\n", -mode); puts( "2" ); exit(0); } /* All constraint partials are computed. Count an evaluation. */ counts.kc += 1; } return; #undef P } /* end of function */