program least3
c
c -- parabolic least-squares fit
c -- input polynomial order from console
c -- subroutines gausj, square, and plot required
c -- figure 7.5
c     
      integer out, maxr, maxc, lines, nrow, ncol
      real x(20), y(20), ycalc(20), coef(3), correl
      real sig(3), resid(20)
      character*1 ch
      common /inout/ out, maxr, maxc
c
      out = 6
      maxr = 20
      maxc = 3
10    call input(x, y, nrow, ncol)
      call linfit(x, y, ycalc, resid, coef, sig,
     *  nrow, ncol, correl)
      call output(x, y, ycalc, resid, coef, sig,
     *  nrow, ncol, correl)
      lines = 2*(nrow - 1) + 1
      write(*,*) 'Press RETURN for plot'
      read(*,'(a1)') ch
      call plot(x, y, ycalc, nrow, out, lines)
      goto 10
      end
      subroutine input(x, y, nrow, ncol)
c -- get values for nrow and arrays x and y
c
      integer nrow, i, out, maxr, maxc, ncol
      real x(1), y(1)
      common /inout/ out, maxr, maxc
c
5     write(out, 101)
      read(*, 102) ncol
      if (ncol .gt. maxc - 1) goto 5
      if (ncol .lt. 1) goto 99
      ncol = ncol + 1
      nrow = 9
      do 10 i = 1, nrow
       x(i) = i
10    continue
      y(1) =  2.07
      y(2) =  8.6
      y(3) = 14.42
      y(4) = 15.8
      y(5) = 18.92
      y(6) = 17.96
      y(7) = 12.98
      y(8) =  6.45
      y(9) =  0.27
      return
99    stop
101   format(' order? ' )
102   format(i2)
      end
      subroutine output(x, y, ycalc, resid, coef, sig,
     *  nrow, ncol, correl)
c
c -- print out the answers
c
      integer out, nrow, ncol, i
      real x(1), y(1), ycalc(1), coef(1), correl
      real resid(1), sig(1)
      common /inout/ out, maxr, maxc
c
      write(out, 101)
      write(out, 102) (i, x(i), y(i), ycalc(i),
     *  resid(i), i = 1, nrow)
      write(out, 103)
      write(out, 106) coef(1), sig(1)
      write(out, 104) (coef(i), sig(i), i = 2, ncol)
      write(out, 105) correl
      return
101   format('   i      x       y     y calc    resid')
102   format(i4, f8.1, 3f9.2)
103   format(/' coefficients    errors')
104   format(0pf8.3, 3x, 1pe12.3)
105   format(/' correlation coefficient is', f7.4)
106   format(0pf8.3, 3x, 1pe12.3, '  constant term')
      end
      subroutine linfit(x, y, ycalc, resid, coef, sig,
     *  nrow, ncol, cor)
c
c -- least squares fit nrow sets of x-y points
c -- using gauss-jordan elimination
c -- subroutines square and gaussj needed
c
      logical error
      integer nrow, ncol, i, j, maxr, maxc, out
      integer index(20,3), nvec
      real x(1), y(1), ycalc(1), coef(1), resid(1)
      real a(3,3), xmatr(20,3), sig(1)
      real sumy, sumy2, xi, yi, yc, res, cor, srs, see
      common /inout/ out, maxr, maxc
      data nvec/1/
c
      do 10 i = 1, nrow
       xi = x(i)
       xmatr(i,1) = 1.0
       xmatr(i,2) = xi
       xmatr(i,3) = xi * xi
10    continue
      call square(xmatr, y, a, coef, nrow, ncol, maxr, maxc)
      call gaussj(a, coef, index, ncol, maxc, nvec, error, out)
      sumy = 0.0
      sumy2 = 0.0
      srs = 0.0
      do 20 i = 1, nrow
       yi = y(i)
       yc = 0.0
       do 15 j = 1, ncol
15      yc = yc + coef(j) * xmatr(i,j)
       ycalc(i) = yc
       res = yc - yi
       resid(i) = res
       srs = srs + res*res
       sumy = sumy + yi
       sumy2 = sumy2 + yi * yi
20    continue
      cor = sqrt(1.0 - srs/(sumy2 - sumy*sumy/nrow))
      if (nrow .eq.ncol) see = sqrt(srs) 
      if (nrow .ne.ncol) see = sqrt(srs / (nrow - ncol))
      do 30 i = 1, ncol
30     sig(i) = see * sqrt(a(i,i))
      return
      end