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