FUNCTION FLOAT XITERATE(X0,F0,X1,F1,X2,F2,ICount,ICvg) /*********************************************************************/ /* */ /* SUBROUTINE: XITERATE */ /* LANGUAGE: FORTRAN 77 */ /* PURPOSE: Iterately solves for the value of X which satisfies */ /* F(X)=0. Given Xi,F(Xi) pairs, the subroutine tests */ /* for convergence and provides a new guess for */ /* the value of the independent variable X. */ /* CODE IS TAKEN FROM ASHRAE HVAC 2 TOOLKIT */ /* */ /*********************************************************************/ /* INPUT VARIABLES */ /* F0 Current value of the function F(X) */ /* X0 Current value of X */ /* F1,F2 Two previous values of F(Xi) */ /* X1,X2 Two previous values of X */ /* ICount Number of iterations */ /* */ /* NOTE: F1,X1,F2,X2 MUST BE STORED AND SAVED IN CALLING */ /* ROUTINE. THEY NEED NO INITIALIZATION */ /* */ /* OUTPUT VARIABLES */ /* XIterate New estimate of X for F(X)=0 */ /* ICvg Convergence flag ICvg = 0: Not converged */ /* ICvg = 1: Converged */ /*********************************************************************/ /* WRITTEN BY: Michael J. Brandemuehl */ /* */ /* ADDAPTED TO NMF BY: Jari Hyttinen */ /* */ /* */ /* DATE: 9 th January 1997 */ /* */ /* CALLS: None */ /* */ /* REFERENCE: ASHRAE HVAC 2 TOOLKIT */ /* Pages 9-3 - 9-7 */ /*********************************************************************/ /* INTERNAL VARIABLES */ /* small Small number used in place of zero */ /* mode Number of points used in fit */ /* mode = 1: Use XPerburb to get new X */ /* mode = 2: Linear equation to get new X */ /* mode > 2: Quadratic equation to get new X */ /* coef(i) Coefficients for quadratic fit */ /* F(X) = coef(1) + coef(2)*X + coef(3)*X*X */ /* check Term under radical in quadratic solution */ /* FiQ,XiQ Double precision values of Fi,Xi */ /* slope Slope for linear fit */ /* tolRel Relative error tolerance */ /* xPerturb Perturbation applied to X to initialize iteration */ /*********************************************************************/ LANGUAGE F77 INPUT FLOAT X0,F0,X1,F1,X2,F2,ICount; OUTPUT FLOAT ICvg; CODE DOUBLE PRECISION FUNCTION XITERATE(X0,F0,X1,F1,X2,F2,ICount,ICvg) DOUBLE PRECISION X0,F0,X1,F1,X2,F2 DOUBLE PRECISION coef(3),check,F0Q,F1Q,F2Q,X0Q,X1Q,X2Q DOUBLE PRECISION small,slope,tolRel,xPerturb,xOther INTEGER ICount,ICvg,mode DATA tolRel/1.E-5/,xPerturb/0.1/,small/1.E-9/ C1*** Check for convergence by comparing change in X IF ((ABS(X0-X1) .LT. tolRel*MAX(ABS(X0),small) .AND. & ICount .NE. 1) .OR. F0 .EQ. 0.) THEN XIterate = X0 ICvg=1 RETURN ENDIF C1*** Not converged. C2*** If after the second iteration there are enough previous points to C2 fit a quadratic for the new X. If the quadratic fit is not C2 applicable, mode will be set to 1 or 2 and a new X will be C2 determined by incrementing X from xPerturb or from a linear fit. ICvg=0 mode=ICount 10 IF (mode .EQ. 1) THEN C1*** New guess is specified by xPerturb IF (ABS(X0) .GT. small) THEN XIterate = X0*(1.+xPerturb) ELSE XIterate = xPerturb ENDIF ELSEIF (mode .EQ. 2) THEN C1*** New guess calculated from LINEAR FIT of most recent two points SLOPE=(F1-F0)/(X1-X0) IF(slope.EQ.0) THEN mode=1 GO TO 10 ENDIF XIterate=X0-F0/SLOPE ELSE C1*** New guess calculated from QUADRATIC FIT C1*** If two Xi are equal, set mode for linear fit and return to top IF (X0 .EQ. X1) THEN X1=X2 F1=F2 mode=2 GO TO 10 ELSEIF (X0 .EQ. X2) THEN mode=2 GO TO 10 ENDIF C1*** Determine quadratic coefficients from the three data points C1*** using double precision. F2Q=F2 F1Q=F1 F0Q=F0 X2Q=X2 X1Q=X1 X0Q=X0 coef(3)=((F2Q-F0Q)/(X2Q-X0Q)-(F1Q-F0Q)/(X1Q-X0Q))/(X2Q-X1Q) coef(2)=(F1Q-F0Q)/(X1Q-X0Q)-(X1Q+X0Q)*coef(3) coef(1)=F0-(coef(2)+coef(3)*X0Q)*X0Q C1*** If points are colinear, set mode for linear fit and return to top IF (ABS(coef(3)) .LT. 1.D-10) THEN mode=2 GO TO 10 ENDIF C1*** Check for precision. If the coefficients do not accurately C1*** predict the given data points due to round-off errors, set C1*** mode for a linear fit and return to top. IF (ABS((coef(1)+(coef(2)+coef(3)*X1Q)*X1Q-F1Q)/F1Q) .GT. & 1.D-4) THEN mode=2 GO TO 10 ENDIF C1*** Check for imaginary roots. If no real roots, set mode to C1*** estimate new X by simply incrementing by xPerturb check=coef(2)**2-4*coef(1)*coef(3) IF (check .LT. 0) THEN C1*** Imaginary roots -- go back to linear fit mode=2 GO TO 10 ELSEIF (check .GT. 0) THEN C1*** Real unequal roots -- determine root nearest to most recent guess XIterate=(-coef(2)+SQRT(check))/coef(3)/2 xOther=-XIterate-coef(2)/coef(3) IF (ABS(XIterate-X0) .GT. ABS(xOther-X0)) XIterate=xOther ELSE C1*** Real Equal Roots -- one solution XIterate=-coef(2)/coef(3)/2 ENDIF ENDIF C1*** Set previous variable values for the next iteration IF (mode .LT. 3) THEN C1*** No valid previous points to eliminate. X2=X1 F2=F1 X1=X0 F1=F0 ELSE C1*** Eliminate one previous point based on sign and magnitude of F(X) C2*** Keep the current point and eliminate one of the previous ones. IF (F1*F0 .GT. 0 .AND. F2*F0 .GT. 0) THEN C2*** All previous points of same sign. Eliminate one with biggest F(X) IF (ABS(F2) .GT. ABS(F1)) THEN X2=X1 F2=F1 ENDIF ELSE C1*** Points of different sign. C1*** Eliminate the previous one with the same sign as current F(X). IF (F2*F0 .GT. 0) THEN X2=X1 F2=F1 ENDIF ENDIF X1=X0 F1=F0 ENDIF RETURN END END_CODE END_FUNCTION