lsq.f

网络带宽测试工具

      Subroutine lsq ( k, x, y, slope, incpt, perc )
! ----------------------------------------------------------------------
! --- Forms Least Squares fit of straight line to data received so
!     far.
! --- Input data in pairs (x,y). Best straight line fitted is described
!     by the output parameters slope and incpt.
!
! --- k      - If k = 0 initialize, else to fit line.
!     x      - Abscissa.
!     y      - Ordinate.
!     slope  - Inverse slope of best line of all data up to now.
!     incpt  - Negative intercept of best line on X-axis.
!     perc   - Percentage of error.
! ----------------------------------------------------------------------
       Implicit              None

       Integer            :: k
       Real(8)            :: x, y, slope ,incpt, perc

       Integer            :: ndata
       Real(8), Parameter :: zero = 0.0d0, one = 1.0d0
       Real(8)            :: permax, permin
       Real(8)            :: dx, dy, rmse, denom, arg
       Real(8)            :: xndata, sumx, sumx2, sumy, sumy2, sumxy

       Save                  xndata, sumx, sumx2, sumy, sumy2, sumxy
! ----------------------------------------------------------------------
! --- Initialise.

      If ( k == 0 ) Then
         xndata = zero
         sumx   = zero
         sumx2  = zero
         sumy   = zero
         sumy2  = zero
         sumxy  = zero
         permin = Huge( one )
         permax = zero
         perc   = zero
         Return
      End If
! ----------------------------------------------------------------------
! --- Update sums with new data x and y.

      xndata = xndata + one
      ndata = Int( xndata )
      sumx  = sumx + x
      sumx2 = sumx2 + x*x
      sumy  = sumy + y
      sumy2 = sumy2 + y*y
      sumxy = sumxy + x*y
      If ( ndata < 3 ) Then
         slope = zero
         incpt = zero
         rmse  = zero
      Else                           ! --- Calculate new slope and incpt.
         denom = xndata*sumx2 - sumx*sumx
         If ( Abs( denom ) < Epsilon( one ) ) Then
            slope = zero
            incpt = zero
            rmse  = zero
         Else
            slope = ( xndata*sumxy - sumx*sumy )/denom
            incpt = ( sumy*sumx2 - sumxy*sumx )/denom
            arg   = ( sumy2 + xndata*incpt**2 + sumx2*slope**2 -
     &                2.0d0*( incpt*sumy - incpt*slope*sumx +
     &                slope*sumxy ) )/( xndata - 2.0d0 )
            If ( arg >= 0 ) Then
               rmse =  Sqrt( arg*sumx2/denom )
            Else
               rmse = -Sqrt(-arg*sumx2/denom )
            End If
          End If
      End If 
! ----------------------------------------------------------------------
! --- Update errors of percentages.

       If ( ndata > 2  ) Then
          perc = 100.0d0*rmse/y
          If ( ndata == 3 ) permin = perc
          permax = Max( permax, perc )
       End If
! ---------------------------------------------------------------------
        End Subroutine lsq