/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:41 */ /*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include "dintdl.h" #include /* COMMON translations */ struct t_dintnc { double ainit, binit, fncval, s, tp, fer, fer1, relobt, tps, xj, xjp; long int fea, fea1, kdim, inc, inc2, istop[2][2], jprint, iprint, kk, kmaxf, ndim, nfindx, nfmax, nfmaxm, reltol, reverm, revers, wherem; LOGICAL32 needh; } dintnc; struct t_dintc { double acum, pacum, result[2], aacum, local[4], abscis, ta, delta, delmin, diff, discx[2], end[2], errina, errinb, fat[2], fsave, funct[24], f2, paacum, pf1, pf2, phisum, phtsum, px, space[6], step[2], start[2], sum, t, tasave, tb, tend, worry[2], x1, x2, x, f1, count, xt[17], ft[17], phi[34], absdif, edue2a, edue2b, ep, epnoiz, epsmax, epso, epsr, epss, errat[2], errc, errf, errt[2], esold, extra, pepsmn, releps, rep, rndc, tlen, xjump, erri, err, epsmin, eps, re, reprod; long int discf, dischk, endpts, inew, iold, ip, ixkdim, j, j1, j1old, j2, j2old, kmax, kmin, l, lendt, nfjump, nsubsv, nxkdim, taloc, where2, i, k, kaimt, nsub, part, search, where, nfeval; LOGICAL32 did1, fail, fats[2], fsaved, havdif, iend, init, roundf, xcdobt[2], pad[7]; } dintc; struct t_dintec { double emeps, eepsm8, edelm2, edelm3, esqeps, ersqep, ersqe6, eminf, esmall, enzer, edelm1, eninf; } dintec; /* end of COMMON translations */ void /*FUNCTION*/ dintdl( double work[]) { long int instop; double pmax; static long dtlent[3]={5,7,15}; static long ib1[3]={1,10,30}; static double xjumpn = 0.25e0; static double xjumps = 0.5e0; static double z5 = 0.01e0; static double betas[133]={.3436491673103708451e001,.2909944487358056275e000, .1000000000000000001e001,.1549193338482966760e001,-.1878361089654305178e000, .2909944487358056232e000,.1878361089654305149e000,.1878361089654305149e000, .1000000000000000015e001,.1830891918707665993e001,.5461818853326139665e000, .1275863152521727801e001,.1877974359956861711e001,-.2908356455650225802e000, .1000000000000000003e001,.1121212539609889453e001,.1411118538882285126e001, .2061029159147656548e000,.7837831181373271228e000,.6990495472071992051e000, .6990495472071992051e000,.7764649827988754109e000,-.2654374897523893370e000, .5461818853326139587e000,.3710664836162815639e000,.2948330583120690783e000, .2654374897523893283e000,.2654374897523893275e000,.1000000000000000036e001, .2160483964289764787e001,.4628592558560086054e000,.1580721165376834188e001, .2788648704047761366e001,-.1659797647456139336e000,.1348061361080862260e001, .1938801234758797259e001,.3001371201600030120e001,.5530131183278168619e-001, .1217373546890437278e001,.1549754626219785878e001,.2074091833232865716e001, .2928362520543975785e001,-.1888472190338949389e-001,.1128427531866524394e001, .1318613567768740554e001,.1600194155392385446e001,.2022328205978644955e001, .2663864875210043425e001,.7089219156396005268e-002,.1059422356252840359e001, .1156722975732715836e001,.1303394479759645425e001,.1518037435357655518e001, .1829746290708357198e001,.2281090415029519236e001,-.3107820325615749339e-002, .1000000000000000007e001,.1028853894914960778e001,.1089361496056109992e001, .1187656139924609753e001,.1334063529367354927e001,.1544285866555422247e001, .1839592669101145570e001,.1689406778911703374e-002,.9439106076041133657e000, .9174389213758399244e000,.9174389213758399244e000,.9435097871391599084e000, .9977420273821314863e000,.1084782404725432559e001,.1211987194763725285e001, .1388740619829111076e001,-.1216502747013758565e-002,.8861889414785074481e000, .8116449626771402972e000,.7665628757411781400e000,.7453813701489683382e000, .7453813701489683365e000,.7659548081345759537e000,.8082259176409903757e000, .8747663621098320137e000,.9688039688731088254e000,.1255674817712366966e-002, .8214405533571194637e000,.7029626866105575496e000,.6238580132970762736e000, .5722253065058340695e000,.5411220158316336791e000,.5265875548985245415e000, .5265875548985245415e000,.5404122752915469377e000,.5682467824057592445e000, .6105250012592685202e000,-.2056713181478912089e-002,.7418059955358485150e000, .5827083724168839593e000,.4801712113071248704e000,.4122773863015805726e000, .3670318230087162626e000,.3375820566477249038e000,.3199261416708740619e000, .3117418541236820569e000,.3117418541236820569e000,.3192575114045456230e000, .3337942651525015294e000,.6161619285278180942e-002,.6326226420594590703e000, .4398667200206020996e000,.3286669728647044621e000,.2600621241862664580e000, .2157588962131737290e000,.1863878189966076094e000,.1668253817995868195e000, .1541355533671636686e000,.1465855112207078905e000,.1431347342606662998e000, .1431347342606662993e000,.1460769330400709832e000,-.4218064520555042749e-001, .4628592558560086436e000,.2623677271576596603e000,.1694821597218475388e000, .1200403163514239665e000,.9113109297367898293e-001,.7309959229942654766e-001, .6136503430084129995e-001,.5355473492811662567e-001,.4835640867710043200e-001, .4500777786787335424e-001,.4304768732135345497e-001,.4218064520555042532e-001, .4218064520555042532e-001,.1000000000000000051e001}; double *const alocal = (double*)dintc.local; double *const blocal = (double*)((double*)dintc.local + 1); double *const fata = (double*)dintc.fat; LOGICAL32 *const fatas = (LOGICAL32*)dintc.fats; double *const fatb = (double*)((double*)dintc.fat + 1); LOGICAL32 *const fatbs = (LOGICAL32*)((LOGICAL32*)dintc.fats + 1); long int *const ic1 = (long*)dintnc.istop; long int *const ic2 = (long*)((long*)dintnc.istop + 2); double *const phit = (double*)((double*)dintc.phi + 17); /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const Betas = &betas[0] - 1; long *const Dtlent = &dtlent[0] - 1; double *const End = &dintc.end[0] - 1; double *const Ft = &dintc.ft[0] - 1; double *const Funct = &dintc.funct[0] - 1; long *const Ib1 = &ib1[0] - 1; double *const Local = &dintc.local[0] - 1; double *const Phi = &dintc.phi[0] - 1; double *const Phit = &phit[0] - 1; double *const Start = &dintc.start[0] - 1; double *const Step = &dintc.step[0] - 1; double *const Work = &work[0] - 1; double *const Worry = &dintc.worry[0] - 1; double *const Xt = &dintc.xt[0] - 1; /* end of OFFSET VECTORS */ /* Copyright (c) 1996 California Institute of Technology, Pasadena, CA. * ALL RIGHTS RESERVED. * Based on Government Sponsored Research NAS7-03001. *>> 1996-04-27 DINTDL Krogh Changes to use .C. and C%%. *>> 1996-03-31 DINTDL Krogh Removed unused variable in common. *>> 1996-03-30 DINTDL Krogh Change specific intrinsics to generics. *>> 1995-11-28 DINTDL Krogh Converted from SFTRAN to Fortran 77. *>> 1994-11-14 DINTDL Krogh Declared all vars. *>> 1994-10-19 DINTDL Krogh Changes to use M77CON *>> 1994-08-22 DINTDL Krogh -- Modified data for C conversion. *>> 1994-07-07 DINTDL Snyder set up for CHGTYP. *>> 1994-05-24 DINTDL Krogh -- Fixes data and max's for C conversion. *>> 1994-05-02 DINTDL Snyder Structured some spaghetti code *>> 1993-05-18 DINTDL Krogh -- Changed "END" to "END PROGRAM" *>> 1991-09-20 DINTDL Krogh converted '(1)' dimensioning to '(*)'. *>> 1987-11-19 DINTDL Snyder Initial code. * * THIS SUBPROGRAM FORMS AND ANALYZES DIFFERENCE LINES FOR DINT * *--D replaces "?": ?INT, ?INTA, ?intc, ?intec, ?INTDL, ?INTNC, ?INTO * * ***** FORMAL ARGUMENT ************************************** * * WORK IS THE USER WORK ARRAY. IT IS NOT USED HERE, BUT IS * PASSED INTO DINTO. */ /* ***** EXTERNAL REFERENCES *********************************** * * DINTO USED FOR PRINTING OUTPUT. * * ***** INTERNAL AND COMMON VARIABLES ************************* * * ALOCAL IS THE LEFT END OF THE CURRENT PANEL. */ /* BETAS ARE COEFFICIENTS USED TO COMPUTE DIFFERENCE LINES */ /* BLOCAL IS THE RIGHT END OF THE CURRENT PANEL. */ /* DTLENT IS A TABLE OF INITIAL DIFFERENCE TABLE LENGTHS. */ /* FATA IS THE FUNCTION VALUE AT ALOCAL. */ /* FATAS INDICATES WHEN FATA HAS USEFUL DATA. */ /* FATB IS THE FUNCTION VALUE AT BLOCAL. */ /* FATBS INDICATES WHEN FATB HAS USEFUL DATA. */ /* IB1 ARE STARTING VALUES FOR INDEXING THE BETAS TABLE */ /* PHIT IS THE BACKWARD DIFFERENCE LINE. PHIT IS EQUIVALENCED * WITH PHI (BELOW). */ /* PMAX IS THE MAXIMUM OF THREE ADJACENT VALUES OF PHI OR PHIT. THIS * IS USED WHEN EXAMINING JUMPS TO AVOID OVERFLOWS. */ /* XJUMPN IS THE VALUE TO USE FOR XJUMP WHEN NOT DOING A SEARCH. */ /* XJUMPS IS THE VALUE TO USE FOR XJUMP WHEN DOING A SEARCH. */ /* Z5 IS AN ADJUSTABLE PARAMETER USED DURING DEVELOPMENT. */ /* SEE DINTA FOR A DESCRIPTION OF REMAINING COMMON VARIABLES * * ***** COMMON STORAGE ****************************************** * * COMMON /DINTNC/ CONTAINS VARIABLES NOT SEPARATELY SAVED FOR * EACH DIMENSION OF A MULTIPLE QUADRATURE. COMMON /DINTC/ * CONTAINS VARIABLES THAT MUST BE SAVED FOR EACH DIMENSION OF THE * QUADRATURE. THE VARIABLES IN EACH COMMON BLOCK ARE STORED IN THE * ORDER - ALWAYS DOUBLE, DOUBLE IF DOUBLE PRECISION PROGRAM, DOUBLE * IF DOUBLE PRECISION PROGRAM AND EXPONENT RANGE OF DOUBLE AND * SINGLE VERY DIFFERENT, SINGLE, INTEGER, LOGICAL. A PAD OF LOGICAL * VARIABLES IS INCLUDED AT THE END OF /DINTC/. THE DIMENSION OF * THE PAD MAY NEED TO BE VARIED SO THAT NO VARIABLES BEYOND THE END * OF THE COMMON BLOCK ARE ALTERED. * * DECLARATIONS OF COMMON /DINTNC/ VARIABLES. * */ /* DECLARATIONS OF COMMON /DINTC/ VARIABLES. * *--D Next line special: S => D, X => Q, D => D, P => D */ /* 139 $.TYPE.$ VARIABLES */ /* Note XT, FT, and PHI above are last, because they must be in adjacent * locations in DINTC. * 30 $DSTYP$ VARIABLES */ /* 29 INTEGER VARIABLES */ /* 11 TO 18 LOGICALS (7 ARE PADDING). */ /* THE COMMON BLOCKS. * */ /* 1 2 3 4 5 6 7 8 * 9 10 11 12 13 1 2 3 * 4 (2,2) 8 9 10 11 12 13 14 * 15 16 17 18 19 20 */ /* 1 2 (4) 6 7 8 9 10 11 (2) * 13 (2) 15 16 17 (2) 19 20 (24) 44 * 45 46 47 48 49 50 51 (6) * 57 (2) 59 (2) 61 62 63 64 65 * 66 (2) 68 69 70 71 72 * 73 (17) 90 (17) 107 (34) */ /* 141 142 143 144 145 146 * 147 148 149 150 (2) 152 153 * 154 (2) 156 157 158 159 160 * 161 162 163 * 164 165 166 167 168 169 */ /* 170 171 172 * 1 2 3 4 5 6 7 8 */ /* THE VARIABLES HERE DEFINE THE MACHINE ENVIRONMENT. ALL ARE SET * IN DINTOP. THE MEANING ATTACHED TO THESE VARIABLES CAN BE * FOUND BY LOOKING AT THE DEFINITIONS IN DINTOP. */ /* ***** EQUIVALENCE STATEMENTS ******************************* * */ /* ***** DATA STATEMENTS *************************************** * */ /* BETAS( 1: 9) are for computing 5-point difference line. * BETAS(10: 29) are for computing 7-point difference line. * BETAS(30:133) are for computing 15-point difference line. * *++ Save data by elements if ~.C. */ /* ***** PROCEDURES ****************************************** * */ if (dintc.where == 0) { /* PREPARE TO COMPUTE THE DIFFERENCES. * */ dintc.havdif = TRUE; dintc.xjump = xjumpn; if (fabs( Worry[dintc.part] - Start[dintc.part] ) < fabs( *blocal - Start[dintc.part] ) || dintc.dischk > 0) dintc.xjump = xjumps; /* PACK THE ABSCISSA AND FUNCTION TABLES. */ dintnc.kk = min( dintc.k, 4 ); /* KK MUST BE RECOMPUTED BECAUSE IT IS NOT SAVED IN EACH DIMENSION */ dintc.lendt = Dtlent[dintnc.kk - 1]; dintc.j1 = dintc.l; if (dintc.k == 2) dintc.j1 = 0; for (dintc.j = 1; dintc.j <= dintc.lendt; dintc.j++) { Xt[dintc.j] = Xt[dintc.j1 + 1]; Ft[dintc.j] = Ft[dintc.j1 + 1]; dintc.j1 += dintc.l; } /* COMPUTE MODIFIED DIVIDED DIFFERENCES. * */ Phit[dintc.lendt] = Ft[dintc.lendt]; Phi[1] = Ft[1]; Phit[1] = Ft[2] - Ft[1]; Phi[2] = -Phit[1]; dintc.l = Ib1[dintnc.kk - 1]; for (dintc.i = 3; dintc.i <= dintc.lendt; dintc.i++) { Phit[dintc.i - 1] = Ft[dintc.i] - Ft[dintc.i - 1]; for (dintc.j = 3; dintc.j <= dintc.i; dintc.j++) { Phit[dintc.i - dintc.j + 1] = Phit[dintc.i - dintc.j + 2] - Betas[dintc.l]*Phit[dintc.i - dintc.j + 1]; dintc.l += 1; } Phi[dintc.i] = Phit[1]*Betas[dintc.l]; dintc.l += 1; } dintc.nfjump = dintc.nfeval; /* COMPUTE THE SUMS OF THE ABSOLUTE VALUES OF THE MIDDLE DIFFERENCES. * */ dintc.l = (dintc.lendt + 1)/2; dintc.phisum = fabs( Phi[dintc.l] ); dintc.phtsum = fabs( Phit[dintc.l] ); if (dintc.k >= 3) { dintc.phisum += fabs( Phi[dintc.l - 1] ) + fabs( Phi[dintc.l + 1] ); dintc.phtsum += fabs( Phit[dintc.l - 1] ) + fabs( Phit[dintc.l + 1] ); } /* IF THE CURRENT PANEL CONTAINS THE END OF A PART, * WE ARE NOT USING A T**N TRANSFORMATION AND PHISUM IS LESS THAN * PHTSUM, TURN THE PROBLEM AROUND AND INTEGRATE IN THE OPPOSITE * DIRECTION. * */ if (dintc.iend) { if (dintc.nsub == 0) { if (dintc.phisum < dintc.phtsum) { /* REVERSE DIRECTION OF INTEGRATION STEP ACCUMULATION. */ if (Worry[dintc.part] == *blocal) Worry[dintc.part] = *alocal; Start[dintc.part] = *blocal; End[dintc.part] = *alocal; *alocal = *blocal; *blocal = End[dintc.part]; Local[3] = *alocal; Local[4] = *blocal; dintc.tend = *blocal; dintc.ta = Start[dintc.part]; dintc.tb = *blocal - dintc.ta; dintc.x1 = dintnc.xj; dintnc.xj = dintnc.xjp; dintnc.xjp = dintc.x1; dintc.j1 = dintc.l - 1; for (dintc.i = 1; dintc.i <= dintc.j1; dintc.i++) { Funct[dintc.i + 17] = Funct[dintc.i + 1] - Funct[dintc.i + 17]; dintnc.s = Xt[dintc.l - dintc.i]; Xt[dintc.l - dintc.i] = Xt[dintc.l + dintc.i]; Xt[dintc.l + dintc.i] = dintnc.s; dintnc.s = Ft[dintc.l - dintc.i]; Ft[dintc.l - dintc.i] = Ft[dintc.l + dintc.i]; Ft[dintc.l + dintc.i] = dintnc.s; dintnc.s = Phi[dintc.l - dintc.i]; Phi[dintc.l - dintc.i] = Phit[dintc.l + dintc.i]; Phit[dintc.l + dintc.i] = dintnc.s; dintnc.s = Phit[dintc.l - dintc.i]; Phit[dintc.l - dintc.i] = Phi[dintc.l + dintc.i]; Phi[dintc.l + dintc.i] = dintnc.s; } dintnc.s = Phi[dintc.l]; Phi[dintc.l] = Phit[dintc.l]; Phit[dintc.l] = dintnc.s; dintc.diff = -dintc.diff; Step[dintc.part] = -Step[dintc.part]; dintnc.s = *fata; *fata = *fatb; *fatb = dintnc.s; dintc.fsaved = *fatas; *fatas = *fatbs; *fatbs = dintc.fsaved; dintc.fsaved = FALSE; /* IF WE ARE TURNING AROUND THE END OF PART 2, HAVE NOT YET DONE * ANYTHING IN PART 2 AND HAVE NO DISCONTINUITY DIAGNOSTIC * PENDING, UNDO THE SUBDIVISION AND RETURN TO MARCHING. */ if (dintc.part == 2) { if (dintc.discf == 0) { if (Start[1] == *blocal) { Start[1] = Start[2]; dintc.part = 1; dintc.iend = FALSE; } } } if (dintnc.iprint > 3) dinto( 4, work ); } } } } /* EXAMINE THE DIFFERENCE LINES. * */ dintc.i = 2; dintc.j = 2; dintnc.kk = 1; dintnc.inc = 1; dintnc.inc2 = 2; instop = dintc.lendt; dintc.j2 = 0; if (!dintc.did1) dintc.j2 = dintc.nsub + 4; dintnc.s = 100.0e0*dintc.epsmin/dintc.delta; if (dintc.search == 4) dintnc.s = 0.0e0; /* do forever */ L_80: ; dintc.l = dintc.i; dintnc.istop[dintnc.kk - 1][0] = dintc.l; dintnc.tp = fabs( Phi[dintc.j - dintnc.inc] ) + z5*dintnc.s; dintnc.xj = dintc.xjump; dintc.j1 = 0; /* While */ L_90: if (dintc.i != instop) { dintc.i += dintnc.inc; dintnc.tp += fabs( Phi[dintc.j] ); dintc.j += dintnc.inc; if (fabs( Phi[dintc.j - dintnc.inc] ) + 2.0e0*fabs( Phi[dintc.j] ) < dintnc.s + fabs( Phi[dintc.j - dintnc.inc2] )) { dintc.l = dintc.i + (dintc.l + dintnc.inc2 - dintc.i)/ 2; if (dintc.l == dintc.i) goto L_80; if (dintc.i <= dintc.j2) goto L_80; } else { if (sign( 1.0e0, Phi[dintc.j] )*Phi[dintc.j - dintnc.inc] >= 0.0e0) goto L_120; } if (dintc.search != 4) dintnc.xj += dintnc.xj; L_120: ; if (fabs( Phi[dintc.j]*dintnc.xj ) <= dintnc.tp) goto L_90; dintnc.xj = fmax( 4.0e-4, dintnc.tp/fabs( Phi[dintc.j] ) ); dintc.j1 = dintc.i; goto L_90; } /* end while */ if (fabs( dintnc.xj ) >= dintc.xjump) dintc.j1 = 0; dintnc.istop[dintnc.kk - 1][1] = dintc.j1; if (dintnc.inc < 0) goto L_140; dintc.j2 = 0; dintnc.xjp = dintnc.xj; instop = 1; dintnc.inc2 = -2; dintnc.inc = -1; dintc.i += dintnc.inc; dintc.j = 16 + dintc.lendt; dintnc.kk = 2; goto L_80; /* end forever */ L_140: ; /* WEAKEN JUMPS IF A LOT OF CONVERGENCE OVERLAP WITH JUMP, THEN TEST * IF JUMP CAN BE IGNORED. */ dintc.j1 = dintnc.istop[0][1]; dintc.j2 = dintnc.istop[1][1]; /* do block */ if (dintc.j1 != 0) { dintc.i = dintc.j1 - *ic2 - 1; if (dintc.i > 0) { dintnc.xjp *= powi(6.0e0,dintc.i); if (dintnc.xjp >= 1.0e0) goto L_160; } dintnc.s = 1.0e-4; if (dintc.j1 != dintc.lendt) dintnc.s = powi(0.1e0,max( dintc.j1 - dintc.j2 + 1, 1 )); pmax = fmax( fabs( Phi[dintc.j1 - 2] ), fmax( fabs( Phi[dintc.j1 - 1] ), fabs( Phi[dintc.j1] ) ) ); if (fabs( dintnc.s*(Phi[dintc.j1]/pmax)*Phi[dintc.j1 - 2] ) >= (Phi[dintc.j1 - 1]/pmax)*Phi[dintc.j1 - 1]) dintnc.xjp = -dintnc.xjp; if (dintc.j1 < dintc.lendt) goto L_170; if (fabs( dintnc.xjp ) < 1.0e-3) goto L_170; if (*ic1 != 2) goto L_170; } L_160: ; dintnc.istop[0][1] = 18; dintnc.xjp = fabs( dintnc.xjp ); /* end block */ L_170: ; /* do block */ if (dintc.j2 != 0) { dintc.i = *ic1 - dintc.j2 - 1; if (dintc.i > 0) { dintnc.xj *= powi(6.0e0,dintc.i); if (dintnc.xj >= 1.0e0) goto L_190; } dintnc.s = 1.0e-4; if (dintc.j2 != 1) dintnc.s = powi(0.1e0,max( dintc.j1 - dintc.j2 + 1, 1 )); pmax = fmax( fabs( Phit[dintc.j2] ), fmax( fabs( Phit[dintc.j2 + 1] ), fabs( Phit[dintc.j2 + 2] ) ) ); if (fabs( dintnc.s*(Phit[dintc.j2]/pmax)*Phit[dintc.j2 + 2] ) >= (Phit[dintc.j2 + 1]/pmax)*Phit[dintc.j2 + 1]) dintnc.xj = -dintnc.xj; if (dintc.j2 >= 2) goto L_200; if (fabs( dintnc.xj ) < 1.0e-3) goto L_200; if (*ic2 != dintc.lendt - 1) goto L_200; L_190: dintnc.istop[1][1] = 0; } dintnc.xj = fabs( dintnc.xj ); /* end block */ L_200: ; if (dintc.j2 == 1) { if (fabs( Phit[2] ) > 100.0e0*fabs( Phit[3] )) dintnc.xj = fabs( dintnc.xj ); } if (dintc.j1 == dintc.lendt) { if (fabs( Phi[dintc.lendt - 1] ) > 100.0e0*fabs( Phi[dintc.lendt - 2] )) dintnc.xjp = fabs( dintnc.xjp ); } if (dintnc.iprint > 3) dinto( 5, work ); return; } /* end of function */