c
c $Id: lls.f,v 1.2 1999/08/20 19:44:14 aai Exp $.
c
c{
      subroutine lls(n,x,y,k,s,a,b)
      integer n
      real x(n)
      real y(n)
      integer k
      real s
      real a
      real b
c
c=============================================================================
c
c
c     Lls - Calculate linear least square parameters.
c
c     Usage:   call lls(n,x,y,k,s,a,b)
c
c     Where:   n = Dimension of parametric arrays x and y.
c	       x = X values.
c	       y = Y values.
c	       k = Returned response indication.
c                   0: No problems - s, a, and b are also returned.
c                  -1: Not enough data for least squares.
c                  -2: Not enough data for calculation of s.
c                  -3: Indefinite slope.
c	       s = Returned standard deviation of y(exp) about y(calc).
c	       a = Returned intercept value.
c	       b = Returned slope.
c
c     The routine returns the values a and b in the equation
c     y(i) = a + b*x(i).
c
c     If(n.le.1) no action is taken by the routine.
c     If(n.le.2) no value is returned for s.
c     If(n*ssx.eq.sx*sx) a, b, and s are not returned (indef or inf
c			 slope) where ssx is the sum of squares of x
c			 and sx is the sum of x values.
c
c     See also: rdot.f and rsum.f - both compiled into the library.
c
c
c     Phillip A. Cheeseman, Purdue University Computing Center.
c
c-----------------------------------------------------------------------------
c
      real pt
      real sx
      real sxx
      real sxy
      real sy
      real syy
      real txx
      real txy
      real tyy
      real xb
      real yb
c
      k = -1
      if(n.le.1) return
c
      sx = rsum(n,x)
      sy = rsum(n,y)
      sxx = rdot(n,x,x)
      syy = rdot(n,y,y)
      sxy = rdot(n,x,y)
c
c     Finish up.
c
      pt = n
      xb = sx/pt
      yb = sy/pt
      txy = pt*sxy-sx*sy
      k = -3
      txx = pt*sxx-sx*sx
      if(txx.le.0.0) return
      tyy = pt*syy-sy*sy
      b = txy/txx
      a = yb-b*xb
      k = -2
      if(n.le.2) return
      s = sqrt(abs(tyy-b*txy)/(pt*float(n-2)))
      k = 0
c
      return
c
      end
c}