c clsq3.f c c
      program leastsq
      implicit none
c
c   Program designed to perform a least squares fit of a quadratic equation
c   to some data.  In this specific example the data is the result of
c   measuring the location of a falling object at various times.  The result
c   of interest is gravitational acceleration obtained from the coefficient
c   off t**2.
c
      real t(50),z(50),coef(0:2),g
      integer np,npoly
c
c   Variables
c
c   t   -  array containing times (seconds) at which positions are measured
c   z   -  measured distance of the falling object from the nominal point of
c          release (fit may let you estimate any offset in the actual release
c          point
c   np  -  number of data points
c   coef - Coefficients in the second order polynomial approximating the data
c   g   -  approximation to the gravitational acceleration deduced from data
c
c
      npoly=2
      call getdata(t,z,np)
      call fit(t,z,np,npoly,coef(0))
      g=2*coef(2)
      write(6,2000) g
 2000 format (' Predicted value of g = ', f6.3,'m/s**2')
      write(6,2001) coef(0),coef(1)
 2001 format(' Predicted initial offset = ',f7.4,'m/s'/
     #         ' Predicted initial velocity = ',f7.4,'m')
      stop
      end
      subroutine getdata(t,s,n )
      implicit none
c
c   Get Experimental Data
c
c   Input Arguments: NONE
c   Output Arguments:
c     t     -    time (sec) of measurement
c     s     -    distance of fall measured (meters)
c     n     -    number of data points
c
      integer n
      real t(*),s(*)
      open(11,file='fall.data')
      n=0
   10 read (11,*,end=20 ) t(n+1),s(n+1)
      n=n+1
      go to 10
   20 return
      end
      subroutine fit (xd,yd,ndata,npoly,c)
      implicit none
c
c    This version of the subroutine has been modified to improve speed
c    on vector or parallel machines, at the cost of extra memory (see
c    the added array xdp)
c
c
c       Take corresponding data points from the arrays xd and yd and fit
c       them with the following equation:
c
c       y = c(1) + c(2) * x + c(3) * x**2  + ... + c(npoly+1)*x**npoly
c
      real xd(ndata),yd(ndata),aa(npoly+1,npoly+1),c(*),fs,
     & sqrt,rsid,xdsum(0:(npoly+npoly)),xdp(ndata,npoly+npoly)
      integer ipvt(npoly+1),ir,is,ij,ndata, info, npoly,i,j,npow, neq
c
c    NOTE:  The use of aa(npoly+1,npoly+1), ipvt(npoly+1), etc. above is
c           an example of the FORTRAN 90 AUTOMATIC ARRAY feature.  It is
c           a short-hand way of using the ALLOCATE statement.  Space for
c           aa, c, and ipvt is assigned when SUBROUTINE FIT is entered and
c           is DEALLOCATED automatically when a RETURN is executed.
c           Note the flexibility available in assigning the array sizes.
c           A compiler without FORTRAN 90 will not permit using npoly to
c           specify the sizes of these arrays unless they (aa, xdsum, ipvt)
c           also appear in the argument list.  Beware of use of this
c           convenient FORTRAN 90 feature in frequently used subroutines.
c           A fair amount of machine time is associated with each ALLOCATE
c           and DEALLOCATE.
c
c    Input Arguments:
c      xd   -   x values for data points
c      yd   -   y valuess for data point pairs
c      ndata -  number of data points
c      npoly -  order of polynomial fit (highest power of x in polynomial)
c
c    Output Arguements:
c      c    -   array containing the coefficients in the chosen order
c               polynomial that provide the "best" fit to the data from a
c               least squares method.  Note that it also temporarily holds
c               the values of the right hand sides of each Least Squares
c               equation.
c    Other Key Variables:
c      xdsum -  array containing the sums of all needed powers of the elements
c               in xd
c      aa   -   matrix containing coefficients of the system of equations
c               generated by the least squares method
c      ipvt -   an array used by Linpack routines to keep track of "pivoting"
c               operations during Gauss Elimination.
c      neq   -  number of linear equations to be solved
c
c      Zero Sum counters
c
      npow=npoly+npoly
      neq=npoly+1
      xdsum(0)=ndata
      do 20 j=1,npow
 20   xdsum(j)=0.
      do 21 j=1,neq
  21  c(j)=0.
c
c   Generate the needed powers of xd(i)
c
      do 30 i=1,ndata
        xdp(i,1)=xd(i)
        do 30 j=2,npow
 30        xdp(i,j)=xdp(i,j-1)*xd(i)
c
c   Sum the Powers of xd
c
      do 40 j=1,npow
c
40       xdsum(j)=sum(xdp(1:ndata,j))
c
c   Generate the right hand sides for all equations
c
      c(1)=sum(yd(1:ndata))
      do 45 j=2,neq
c
   45       c(j)=dot_product(yd(1:ndata),xdp(1:ndata,j-1))
c
c       DO loop 50 generates the coefficients a given equation
c
      do 50 i=1,neq
         do 50 j=1,neq
  50        aa(i,j)=xdsum(i+j-2)
c
c    Solve the Least Squares Equations
c
      call sgefa(aa,neq,neq,ipvt,info)
      call sgesl(aa,neq,neq,ipvt,c ,0)
c
c    Calculate a measure of mean error between the data and the curve
c
      rsid=0.
      do 65 ir=1,ndata
        fs=0.
        do 60 is=1,neq
   60      fs=fs+c(is)*xd(ir)**(is-1)
   65   rsid=rsid+(yd(ir)-fs)**2
      rsid=sqrt(rsid/float(ndata-1))
      write(6,2001) rsid
 2001 format(' Fit to 2nd order polynomial has a mean error of',
     $         1p,e12.5)
      return
      end
c
c c