SUBROUTINE simlp(n, x, y, sad, alpha, beta, d, iter, inext, ifault) ! Code converted using TO_F90 by Alan Miller ! Date: 2003-06-19 Time: 09:37:54 ! ALGORITHM AS 132 APPL. STATIST. (1978) VOL.27, NO.3 ! SIMLP: Fit Y = ALPHA + BETA.X + error IMPLICIT NONE INTEGER, INTENT(IN) :: n REAL, INTENT(IN OUT) :: x(n) REAL, INTENT(IN) :: y(n) REAL, INTENT(OUT) :: sad REAL, INTENT(OUT) :: alpha REAL, INTENT(OUT) :: beta REAL, INTENT(OUT) :: d(n) INTEGER, INTENT(OUT) :: iter INTEGER, INTENT(OUT) :: inext(n) INTEGER, INTENT(OUT) :: ifault REAL, PARAMETER :: acu = 1.0E-06, big = 1.0E19, half = 0.5, zero = 0.0, & one = 1.0, two = 2.0 INTEGER :: i, ibas1, ibas2, iflag, iin, iout, isave, j REAL :: a1, a2, aaa, aaaa, ahalf, aone, bbb, bbbb, ddd, det, & ratio, rho, subt, sum, t, test, tot1, tot2, y1, y2, zzz ! Initial settings ifault = 0 iter = 0 ahalf = half + acu aone = ahalf + ahalf ! Determine initial basis d(1) = zero y1 = y(1) ibas1 = 1 a1 = x(1) DO i = 2, n IF (ABS(a1-x(i)) >= acu) THEN a2 = x(i) ibas2 = i y2 = y(i) GO TO 20 END IF END DO ifault = 1 RETURN ! Calculate initial beta value 20 det = one / (a2-a1) aaaa = (a2*y1-a1*y2) * det bbbb = (y2-y1) * det ! Calculate initial D-vector DO i = 2, n ddd = y(i) - (aaaa+bbbb*x(i)) d(i) = SIGN(one,ddd) END DO tot1 = one tot2 = x(ibas2) d(ibas2) = -one DO i = 2, n tot1 = tot1 + d(i) tot2 = tot2 + d(i) * x(i) END DO t = (a2*tot1-tot2) * det IF (ABS(t) >= aone) THEN det = -det GO TO 70 END IF ! Main iterative loop begins 50 t = (tot2-a1*tot1) * det IF (ABS(t) < aone) GO TO 130 iflag = 2 iout = ibas2 x(iout) = a1 aaa = a1 bbb = a2 GO TO 80 60 t = (tot2-a2*tot1) * det IF (ABS(t) < aone) GO TO 130 70 iflag = 1 bbb = a1 aaa = a2 iout = ibas1 80 rho = SIGN(one,t) t = half * ABS(t) det = det * rho ! Perform partial sort of ratios inext(ibas1) = ibas2 ratio = big sum = ahalf DO i = 1, n ddd = (x(i)-aaa) * det IF (ddd*d(i) > acu) THEN test = (y(i)-aaaa-bbbb*x(i)) / ddd IF (test < ratio) THEN j = ibas1 sum = sum + ABS(ddd) 90 isave = ABS(inext(j)) IF (test < d(isave)) THEN IF (sum >= t) THEN subt = ABS((x(isave)-aaa)*det) IF (sum-subt >= t) THEN sum = sum - subt d(isave) = SIGN(1,inext(j)) inext(j) = inext(isave) GO TO 90 END IF END IF 100 j = isave isave = ABS(inext(j)) IF (test < d(isave)) GO TO 100 END IF inext(i) = inext(j) inext(j) = SIGN(i,INT(d(i))) d(i) = test IF (sum >= t) THEN iin = ABS(inext(ibas1)) ratio = d(iin) END IF END IF END IF END DO ! Update basic indicators iin = ABS(inext(ibas1)) j = iin 120 isave = ABS(inext(j)) IF (isave /= ibas2) THEN zzz = SIGN(1,inext(j)) tot1 = tot1 - zzz - zzz tot2 = tot2 - two * zzz * x(isave) d(isave) = -zzz j = isave GO TO 120 END IF zzz = SIGN(1,inext(ibas1)) tot1 = tot1 - rho - zzz tot2 = tot2 - rho * bbb - zzz * x(iin) d(iout) = -rho iter = iter + 1 IF (iflag /= 1) THEN x(ibas2) = a2 ibas2 = iin d(ibas2) = -one a2 = x(iin) y2 = y(iin) det = one / (a1-a2) aaaa = (a1*y2-a2*y1) * det bbbb = (y1-y2) * det GO TO 60 END IF ibas1 = iin a1 = x(iin) d(ibas1) = zero y1 = y(iin) det = one / (a2-a1) aaaa = (a2*y1-a1*y2) * det bbbb = (y2-y1) * det GO TO 50 ! Calculate optimal sum of absolute deviations 130 sad = zero DO i = 1, n d(i) = y(i) - (aaaa+bbbb*x(i)) sad = sad + ABS(d(i)) END DO alpha = aaaa beta = bbbb RETURN END SUBROUTINE simlp