PRACE NPT

      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