/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:23 */ /*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include "szero.h" #include #include #include /* PARAMETER translations */ #define C0 0.e0 #define C1 1.e0 #define C1P031 1.03125e0 #define C1P25 1.25e0 #define C2 2.e0 #define C4 4.e0 #define C8 8.e0 #define CP001 0.001e0 #define CP01 0.01e0 #define CP125 0.125e0 #define CP25 0.25e0 #define CP5 0.5e0 #define CP75 0.75e0 #define CP99 0.99e0 #define LINIT (-40) #define LTXTAB 8 #define LTXTAC 63 #define LTXTAD 112 #define LTXTAE 169 #define MEEMES 52 #define MERET 51 #define METEXT 53 /* end of PARAMETER translations */ void /*FUNCTION*/ szero( float *x1, float *f1, float *x2, float *f2, long *mode, float tol) { long int _l0, i, idat[2], j; static long int knkp[2], lnlp[2]; float dfdxxo, dfdxxx, dxdfff, dxdffo, fdat[4], ff, ffdfo, ffmfb, ffmfl, qfm, qxm, tolx, tp, tp1, xrnd, xxmxl, xxmxo; static float fl, flmfb, fo, small, xl, xlmxb, xo, xx, xxmxol; static char mtxtaa[1][194]={"SZERO$BBest bound for zero is [$F, $F], but tolerance is $F.$EApparent discontinuity in function near X = $F.$ECan not$ find a sign change: X1=$F, X2=$F, F1=$F, F2=$F$ECalled with MODE$ = $I.$E"}; static char mtxtab[1][47]={"In SZERO -- X1=$F F1=$F KTYP=$I DIV=$F KS=$I$E"}; static char mtxtac[1][26]={" X2=$F F2=$F$E"}; static long mact[5]={MEEMES,0,0,0,MERET}; static long mact1[2]={METEXT,MERET}; static long mloc[4]={LTXTAB,LTXTAC,LTXTAD,LTXTAE}; static float rnd = C0; static long ks = 0; static long ktyp = 0; static long lmode = 2; static float div = C0; static long lchg[30]={1,2,0,0,0,1,1,4,5,0,1,2,3,3,0,1,4,4,3,0, 1,4,5,5,0,1,5,5,5,0}; static long np = 0; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const Fdat = &fdat[0] - 1; long *const Idat = &idat[0] - 1; long *const Knkp = &knkp[0] - 1; long *const Lchg = &lchg[0] - 1; long *const Lnlp = &lnlp[0] - 1; long *const Mact = &mact[0] - 1; long *const Mact1 = &mact1[0] - 1; long *const Mloc = &mloc[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. *>> 2010-04-14 SZERO Krogh No discontinuity message if |F1-F2| small. *>> 2010-04-12 SZERO Krogh Fixed KNKP to get discontinuity diagnostic. *>> 2010-02-20 SZERO Krogh $G => $F for print out of iterations *>> 2008-03-01 SZERO Krogh Minor change in diagnostic print. *>> 2000-12-01 SZERO Krogh Removed unused variable C1P01. *>> 1998-11-01 SZERO Krogh Set so errors stop less easily. *>> 1998-11-01 SZERO Krogh For "mangle", INDIC replaced with MACT(3). *>> 1996-03-30 SZERO Krogh Added external statement. *>> 1995-11-09 SZERO Krogh Fixed so char. data at col. 72 is not ' '. *>> 1994-11-11 SZERO Krogh Declared all vars. *>> 1994-10-20 SZERO Krogh Changes to use M77CON *>> 1994-09-08 SZERO Krogh Added CHGTYP code. *>> 1993-04-27 SZERO Krogh Additions for Conversion to C. *>> 1993-04-13 SZERO Krogh Minor change for new MESS. *>> 1992-04-08 SZERO Krogh Unused label 400 removed. *>> 1992-01-09 SZERO Krogh Moved calc. of XXMXO up (for error msg.) *>> 1991-11-26 SZERO Krogh Converted to new error processor. *>> 1988-08-14 SZERO Krogh Labels runumbered. *>> 1988-03-07 SZERO Krogh Initial code. * *--S replaces "?": ?ZERO, ?MESS * * SUBROUTINE TO FIND A BOUNDED ZERO * * analysis and coding by Fred T.Krogh at the Jet Propulsion * Laboratory, Pasadena, Calif. April 25, 1972. * Modified for portability, April 1984 by Krogh. * Algorithmic changes, vars. added to save stmt., Sept. 1987 by Krogh * * Parameters in the calling sequence are defined as follows: * * X1 = independent variable * F1 = dependent variable -- initially F1=F(X1). * When MODE=1 (or 5) the user is to compute F(X1) given X1 * X2 = second value of independent variable * F2 = F(X2) on the initial entry. When MODE = 2-4, F2=F(X2) and * F1*F2 .le. 0. * MODE is a parameter used for communication between this * subroutine and the user. (The user should set MODE * only to initialize it to 0 before the first call) * =1 compute F(X1) and call $ZERO * =2 F(X1) is approximately 0, the iteration is finished * and the error criterion is satisfied. * =3 same as MODE=2, except the error criterion can * not be satisfied. * =4 apparently the function has a discontinuity * between X1 and X2 -- No zero can be found * =5 F1*F2 was greater than zero on the first call, and an attempt * to bound the zero on both sides have failed. * =6 fatal error -- $ZERO was called after mode was set .ge.2. * If $ZERO is called again, the program will be stopped. * (Unless MODE is set to 0) * <0 If MODE is set <0 and $ZERO is called, no action is taken * except that print is turned on for -MODE calls to $ZERO. * This print gives all values of X and F used in the iteration. * TOL is the error tolerance * TOL.GT.0 Iterate until values of X1 and X2 are known * for which abs(X1-X2) .le. tol and F1*F2 .le. 0. * TOL.LT.0 Iterate until a value of X1 is found for which * abs(F1) .le. abs(TOL). * TOL = 0 Iterate until the zero is determined as * precisely as possible. MODE = 3 is impossible * in this case. * * Parameters in the calling sequence have the following types * */ /* Usage is as follows (of course, variations are possible.) * Somehow one has available X1, F1, X2, and F2 such * that F1 = F(X1), F2 = F(X2) and F1*F2 .le. 0. * In addition, one should assign a value to TOL. * MODE = 0 **** In the statement below, $ is replaced by an 'S' for single **** precision and a 'D' for double. * XXX call $ZERO(X1,F1,X2,F2,MODE,TOL) * go to (N1,N2,N3,N4,N5,N6), MODE * N1 COMPUTE F1=F(X1) * go to XXX * * N4 continue * N5 continue * N6 stop * N3 If you want to -- print results to note that error * is too big. * N2 zero is found, do whatever you want to with it. * * End of comments explaining usage. * * ************************* Usage of internal variables **************** * * C0 Parameter = 0. * C1 Parameter = 1. * C1P25 Parameter = 1.25 * C2 Parameter = 2. * C4 Parameter = 4. * CP01 Parameter = 0.01 * CP125 Parameter = 1.25 * CP25 Parameter = 0.25 * CP5 Parameter = 0.5 * CP75 Parameter = 0.75 * CP99 Parameter = 0.99 * R1MACH Gets constants associated with floating point arithmetic. * DFDXXX = (XXMXL/FFMFL) * (est. deriv. of f w.r.t. x at x = XX). All * derivatives are base on a second degree polynonial that interpolates * the last three points generated. * DFDXXX = (XXMXO/FFMFL) * (est. deriv. of f w.r.t. x at x = X0). * DIV If artificial subdivision of the interval is used, determines * the amount of the sudivision. (-XXMXOO * DIV / (1. + DIV)) * SMESS Prints error messages. * DXDFFF = (FFMFL/XXMXL) * (est. deriv. of x w.r.t. f at f = FF). * DXDFFO = (FFMFO/XXMXL) * (est. deriv. of x w.r.t. f at f = FO). * F1 (formal arg.) The last value of F computed, on return the value * of F(X1). * F2 (formal arg.) The other initial value provided for F. Set to * the value of F(X2) on returns. * FDAT Temporary storage for floating point values for messages. * FF Value of F1 after F is computed. * FFDFO FF / FO * FFMFB FFMFL + FLMFB = FF - value of FF 2 iterations back. * FFMFL FF - FL * FL Value of FF from the previous iteration. * FLMFB Value of FFMFL from the previous iteration * FO F(XO) * I Comments for LCHNG define how I is set. * IDAT Temporary storage for integer values for messages. * J This is 1 if FF .le. 0., and is 2 if FF > 0. * KNKP KNKP(J) (see J above) is decreased by 3 when there are signs of * decent convergence. It is counted up when convergence is slow. * KS =-1 initially, * = 0 whenever F changes sign, otherwise * = number of times the sign of F has remained the same * KTYP = 1 if interpolation was used to get the last iterate, = 0 if * an artificial subdivision was used. * LCHG the J-th continuation in the data statement for LCHG below gives * new states for the case when the state number is J-1. State 0 is the * initial state. The I-th entry on a row gives the state for case on I * as follows: (TP is the ratio (new f) / (last f of same sign) * I = 1 TP < 0.01 * I = 2 .01 <= TP < 1 * I = 3 TP = 1 * I = 4 1 < TP <= 4 * I = 5 TP > 4. * States are as follows: * 0 initial state, or big increase, or small increase in state 0 * 1 after big decrease, perhaps followed by small decreases * 2 after one or more small decreases from state 0 * 3 one or more small increases from state 2 * 4 one or more small decreases from state 3 * 5 decision made that noise is a problem on this side of zero. * LINIT - the maximum number of iterations that can be taken with out * getting a sign change in F. * LMODE The value of MODE the last time in this routine. * LNLP * LTXTxx Names of this form are used in setting up data statements for * error messages. These names are generated automatically by PMESS, * the program that makes up these messages. * MACT This array difines the actions to be taken by the error message * program. See comments in MESS for details. MODE is set to MACT(3) * on exit. * MACT1 As for MACT except used for the diagnostic print. * MExxxx Parameters defining constants used for interaction with the * error message program MESS. See comments there for definitions. * MLOC Contains locations in MTXTAA for error messages. * MODE (formal) See comments above. * MTXTAA Text for error messages. * MTXTAB Text for diagnostic message. */ /* MTXTAC Text for diagnostic message. * NP If > 0, gives number of iterations till diagnostic print stops. * QFM * QXM * RND Largest relative difference between succesive floating point * numbers. * SMALL .5 / (RND * largest floating point number) * TOL (Formal) See description above. * TOLX Actually tolerance required for accuracy in X. Usually = * max(TOL, XRND). It can be cut by a factor of 2 for use in setting * bounds on an acceptable interval. * TP Ordinarily the ratio (FF / prev. FF of the same sign. * TP1 Used for temporary storage. * X1 (Formal) Value of x where F is to be computed, and value * returned for the zero after convergence. * X2 (Formal) Initially other value of x where F is given. After * convergence gives the other closest x which gives an F of opposite * sign from that given by x1. * XL Value of XX from the previous iteration. * XLMXB Value of XXMXL from the previous iteration. * XO Value of x on opposite side of the zero from the current x. * XRND Best accuracy that one could hope for based on the finite * precision of floating point numbers. * XX Current x, the last value of X1 where F was computed. * XXMXL XX - XL * XXMXO XX - XO = length of interval in which 0 lies. * XXMXOL Value of XXMXO from a previous iteration. * */ /* Declarations for error message processing. * */ /* ********* Error message text *************** *[Last 2 letters of Param. name] [Text generating message.] *AA SZERO$B *AB Best bound for zero is [$F, $F], but tolerance is $F.$E *AC Apparent discontinuity in function near X = $F.$E *AD Can not find a sign change: X1=$F, X2=$F, F1=$F, F2=$F$E *AE Called with MODE = $I.$E * $ *AF In SZERO -- X1=$F F1=$F KTYP=$I DIV=$G KS=$I$E * $ *AG X2=$F F2=$F$E */ /* **** End of automatically generated text * 1 2 3 4 5 */ /* INITIALIZE * */ if (*mode < 0) { np = -1 - *mode; return; } if (np > 0) { np -= 1; Fdat[1] = *x1; Fdat[2] = *f1; Fdat[3] = div; Idat[1] = ktyp; Idat[2] = ks; smess( mact1, (char*)mtxtab,47, idat, fdat ); if (*mode != 0) switch (IARITHIF(lmode - 1)) { case -1: goto L_70; case 0: goto L_80; case 1: goto L_450; } Fdat[1] = *x2; Fdat[2] = *f2; smess( mact1, (char*)mtxtac,26, idat, fdat ); } else if (*mode != 0) { switch (IARITHIF(lmode - 1)) { case -1: goto L_70; case 0: goto L_80; case 1: goto L_450; } } if (rnd == C0) { rnd = FLT_EPSILON; small = CP5/(rnd*FLT_MAX); } xl = *x2; fl = *f2; L_30: tp = C1; *mode = 1; Mact[3] = 2; xxmxol = C0; Knkp[1] = 0; Knkp[2] = 0; Lnlp[1] = 0; Lnlp[2] = 0; ks = -1; xx = *x1; ff = *f1; switch (SARITHIF(fl)) { case -1: goto L_40; case 0: goto L_75; case 1: goto L_50; } L_40: switch (SARITHIF(ff)) { case -1: goto L_60; case 0: goto L_230; case 1: goto L_100; } L_50: switch (SARITHIF(ff)) { case -1: goto L_100; case 0: goto L_230; case 1: goto L_60; } L_60: lmode = 0; /* Take care of points on same side of zero. */ L_70: ff = *f1; xx = *x1; tp = ff/fl; if (tp < C0) goto L_30; lmode -= 1; if (lmode < LINIT) { Mact[3] = 5; Fdat[1] = xx; Fdat[2] = xl; Fdat[3] = ff; Fdat[4] = fl; goto L_250; } if (tp > C1) { ff = fl; xx = xl; fl = *f1; xl = *x1; } if (fabsf( ff ) >= C8*fabsf( fl - ff )) { tp = C8; } else { tp = fmaxf( -CP25*(float)( lmode ), ff/(fl - ff) ); } fl = ff; xo = xl; xl = xx; if (xx == xo) xo = C1P031*xx + signf( CP001, xx ); xx += tp*(xx - xo); *x1 = xx; *mode = 1; return; L_75: *x1 = xl; *f1 = fl; goto L_250; /* END OF INITIALIZATION * * * ENTRY AFTER COMPUTING F FOR THE LAST ITERATE */ L_80: ff = *f1; tp = ff/fl; switch (SARITHIF(tp)) { case -1: goto L_90; case 0: goto L_230; case 1: goto L_110; } L_90: tp = ff/fo; ks = 0; L_100: fo = fl; xo = xl; goto L_120; L_110: ks += 1; L_120: j = 1; if (ff > C0) j = 2; switch (SARITHIF(tp - C1)) { case -1: goto L_150; case 0: goto L_140; case 1: goto L_130; } L_130: i = 4; if (tp > C4) i = 5; goto L_160; L_140: i = 3; goto L_160; L_150: i = 2; if (tp < CP01) i = 1; if (tp < CP99) goto L_170; L_160: Knkp[j] += 1; goto L_180; L_170: Knkp[j] = 0; L_180: xxmxo = xx - xo; Lnlp[j] = Lchg[5*Lnlp[j] + i]; if (Lnlp[j] >= 4) { if (Lnlp[3 - j] >= 4) goto L_210; } /* XXMXO GIVES THE LENGTH OF THE INTERVAL INSIDE WHICH * THE ZERO IS KNOWN TO LIE. */ if (C2*fabsf( xxmxo ) < fabsf( xxmxol )) { Knkp[j] = max( 0, Knkp[1] - 3 ); } xxmxol = xxmxo; xrnd = rnd*(fabsf( xx ) + fabsf( xo ) + small); /* TEST FOR CONVERGENCE */ switch (SARITHIF(tol)) { case -1: goto L_190; case 0: goto L_200; case 1: goto L_200; } L_190: ; if (fabsf( ff ) <= fabsf( tol )) goto L_220; L_200: ; tolx = fmaxf( tol, xrnd ); if (fabsf( xxmxo ) > tolx) goto L_310; /* CONVERGENCE -- PREPARE FOR FINAL EXIT */ L_210: if ((fabsf( xxmxo ) > tol) && (tol != C0)) { Mact[3] = 3; Fdat[3] = tol; if (xxmxo > 0) { Fdat[2] = xx; Fdat[1] = xo; } else { Fdat[1] = xx; Fdat[2] = xo; } } /* SET FINAL VALUES FOR X1,F1,X2,AND F2 */ L_220: ; if (fabsf( ff ) <= fabsf( fo )) goto L_240; *f1 = fo; *x1 = xo; L_230: fo = ff; xo = xx; L_240: *x2 = xo; *f2 = fo; /* TEST FOR DISCONTINUITY */ if ((Knkp[1] > 5) || (Knkp[2] > 5)) { if (fabsf( *f1 - *f2 ) > rnd*fmaxf( *x1, 1.e0 )) { Mact[3] = 4; Fdat[1] = xx; } } L_250: *mode = Mact[3]; switch (IARITHIF(Mact[3] - 2)) { case -1: goto L_420; case 0: goto L_420; case 1: goto L_430; } /* END OF CODE FOR FINAL EXIT * * F NOT DECREASING (OR THE FIRST ITERATE) * PREPARE TO DIVIDE THE INTERVAL */ L_260: tp = C1; switch (IARITHIF(ks)) { case -1: goto L_370; case 0: goto L_280; case 1: goto L_270; } L_270: if (ktyp == 0) goto L_290; L_280: div = C2; L_290: ; div = fmaxf( div, ffdfo ); /* KTYP=0 IF AND ONLY IF THE INTERVAL WAS DIVIDED (USING DIV) * ON THE LAST ITERATION */ if (ktyp == 0) div *= C1P25/(C1P25 - tp); /* DIVIDE THE INTERVAL AS SPECIFIED BY DIV */ L_300: tp1 = -xxmxo*(div/(div + C1)); ktyp = 0; goto L_410; L_310: ; xxmxl = xx - xl; ffmfl = ff - fl; ffdfo = fabsf( ff/fo ); tolx *= CP5; if (tp >= C1) goto L_260; /* DIVIDE THE INTERVAL IF F HAS HAD THE SAME SIGN FOR * FOUR OR MORE TIMES IN SUCCESSION */ switch (IARITHIF(ks - 4)) { case -1: goto L_320; case 0: goto L_340; case 1: goto L_290; } L_320: ; if (flmfb == C0) goto L_340; /* BEGINNING OF CODE TO DETERMINE IF INVERSE QUADRATIC * INTERPOLATION IS TO BE USED. */ ffmfb = ffmfl + flmfb; if (ffmfb == C0) goto L_330; qfm = C1 - (ffmfl/flmfb)*(xlmxb/xxmxl); qxm = C1 - (xxmxl/xlmxb)*(flmfb/ffmfl); dxdfff = C1 + (ffmfl/ffmfb)*qfm; dxdffo = dxdfff + C2*((fo - ff)/ffmfb)*qfm; tp1 = xxmxl + xlmxb; dfdxxx = C1 + (xxmxl/tp1)*qxm; dfdxxo = dfdxxx + C2*((xo - xx)/tp1)*qxm; tp1 = dxdfff*dfdxxx; if ((tp1 <= CP25) || (tp1 >= C4)) goto L_330; tp1 = dxdffo*dfdxxo; if ((tp1 > CP25) && (tp1 < C4)) goto L_380; /* DERIVATIVES DO NOT MATCH WELL ENOUGH */ L_330: ; if (ks == 0) switch (SARITHIF(ffdfo - C1)) { case -1: goto L_350; case 0: goto L_370; case 1: goto L_360; } L_340: ; if ((ktyp == 0) && (tp >= CP75)) goto L_290; ; tp = C1 - tp; if (tp <= ffdfo) goto L_280; ffdfo /= tp; div = CP125; goto L_290; L_350: ; div = CP5*fmaxf( fmaxf( CP25, ffdfo ), tp/(C1P25 - fminf( tp, C1 )) ); goto L_300; L_360: ; div = fminf( C4, CP5*ffdfo ); goto L_300; /* INTERPOLATE WITH SECANT METHOD */ L_370: tp1 = -xxmxl; goto L_390; /* DERIVATIVES MATCH UP PRETTY WELL. */ L_380: ; /* INTERPOLATE USING THE INVERSE QUADRATIC */ tp1 = xxmxl*(qfm*(fl/ffmfb) - C1); L_390: tp1 *= ff/ffmfl; ktyp = 1; /* EXIT TO GET F(X) */ L_410: ; fl = ff; flmfb = ffmfl; xlmxb = xxmxl; xl = xx; /* COMPUTE X1, INSURING THAT IT IS NOT TOO CLOSE TO THE * ENDS OF THE INTERVAL */ xx = fminf( fmaxf( xl + tp1, fminf( xl, xo ) + tolx ), fmaxf( xl, xo ) - tolx ); *x1 = xx; L_420: lmode = *mode; return; L_430: Mact[2] = 11*Mact[3] - 19; L_440: Mact[4] = Mloc[Mact[3] - 2]; smess( mact, (char*)mtxtaa,194, idat, fdat ); goto L_420; /* A CALL TO THE SUBROUTINE HAS BEEN MADE WITH MODE.NE.1 */ L_450: Idat[1] = *mode; Mact[3] = 6; *mode = 6; if (lmode != 6) goto L_430; Mact[2] = 99; goto L_440; } /* end of function */