***POLYNOMIAL REGRESSION***********************************************
*                 PROGRAM FOR POLYNOMIAL REGRESSION                   *
***********************************************************************
*
      program POLREGRESSION
      implicit doubleprecision (a-h,o-z)
      parameter(nd=40)
      doubleprecision X(nd),Y(nd),S(nd),A(nd,nd)
      integer l,r
*
      do i = 1, 20
         S(i) = 0
      enddo
*
      sy = 0
      st = 0
      sr = 0
*
      print*, 'Enter the number of data pairs, n:'
      read(*,*) n
      print*, 'Enter the degree of regression polynomial, m:'
      read(*,*) m
10    if (m >= n) then
         write(6,*) 'WARNING: *The degree of polynomial must be less tha
     $n number of data points.*'
         write(6,*)
         print*, 'Enter the degree of regression polynomial, m:'
         read(*,*) m
         go to 10
      endif
*
      print*, 'Enter X data:'
      read(*,*) (X(i), i = 1,n)
      print*, 'Enter Y data:'
      read(*,*) (Y(i), i = 1,n)
      write(6,*)'------------------------------------------------------'
*
      do i = 1, n
         sy = sy + Y(i)
      enddo
      ymean = sy/n
*
      mm = m + 1
      r = mm + 1
*
      do i = 1, mm
      do j = 1, r
         A(i,j) = 0
      enddo
      enddo
*
      do i = 1, mm
         do j = 1, mm
            k = i + j - 2
            do l = 1, n
               A(i,j) = A(i,j) + X(l)**k
            enddo 
         enddo
         do l = 1, n
            A(i,r) = A(i,r) + (Y(l) * X(l)**(i - 1))
         enddo 
      enddo 
*
      CALL GAUSS(A,mm,S)
*
      write(6,*) 'The coefficients of the least square polynomial equati
     $on are:'
      write(6,*)
      write(6,101) (S(i), i = 1,mm)
101   format('',2x,f9.4)
      write(6,*)
      write(6,*) 'THE POLYNOMIAL FIT IS:'
      write(6,102) 'Y  =',S(1),'+',S(2),'X','+',S(3),'X **2','+',S(4
     $),'X **3','+',S(5),'X **4','+',S(6),'X **5','+',S(7),'X **6','+
     $',S(8),'X **7','+',S(9),'X **8'
102   format('',1x,a4,f8.4,a3,f8.4,a2,7(a3,f8.4,a6)) 
*
      do i = 1, n
         st = st + (Y(i) - ymean) **2
         sr = sr + (Y(i) - S(1) - S(2)*X(i) - S(3)*X(i)**2 - S(4)
     $        *X(i)**3 - S(5)*X(i)**4 - S(6)*X(i)**5 - S(7)*X(i)**6 - 
     $         S(8)*X(i)**7 - S(8)*X(i)**8)**2
      enddo
      write(6,*)
      if (n.ne.(m+1)) then
         serror = sqrt(sr/(n-m-1))
         write(6,*) 'Standard error of estimate:',serror
      endif
      cdet   = (st - sr)/st
      correl = dsqrt(dabs(cdet))
      write(6,*) 'Coefficient of correlation:',correl
      write(6,*) 'Coefficient of determination:',cdet
      write(6,*)'------------------------------------------------------'
*
      stop
      end
*
***GAUSS ELIMINATION**************************************************
*
      subroutine GAUSS(A,n,Y)
      implicit doubleprecision(a-h,o-z)
      parameter (nd = 40)
      doubleprecision A(nd,nd),X(nd),Y(nd),ORDC(nd),ORD(nd)
      integer p,r
      parameter(epsil = 1d-9)
*
      m = n + 1
      nn = n - 1
c
c     establish initial order in the column order vector
      do i = 1, n
         ORDC(i) = i
      enddo
c
c     segment for partial pivoting
      do p = 1, nn
         CALL PIVOT(A,n,ORD,ORDC,p)
c
c     triangularization by eliminating the variables
         kk = p + 1
         do i = kk, n
            if (dabs(A(p,p)) <= epsil) then
               write(*,*)'-----------------------------------------------'
               write(*,*)'Coefficient matrix singular. No Solution exist.'
               write(*,*)'-----------------------------------------------'
               stop      ! terminating the program
            else
               qt = A(i,p)/A(p,p)
            endif 
            do j = p, m
               A(i,j) = A(i,j) - qt *A(p,j)
            enddo
         enddo
      enddo
c
c     checking for the singularity of coefficient matrix
      do i = 1, n
         if (dabs(A(i,i)) <= epsil) then
            write(*,*)'-----------------------------------------------'
            write(*,*)'Coefficient matrix singular. No Solution exist.'
            write(*,*)'-----------------------------------------------'
            stop      ! terminating the program
         endif
      enddo
c
c     back substitution
      X(n) = A(n,m)/A(n,n)
      do i = nn, 1, -1
         sum = 0
         k = i + 1
         do j = k, n
            sum = sum + A(i,j) *X(j)
         enddo
         X(i) = (A(i,m) - sum)/A(i,i)
      enddo
c
c     rearrange the solution vector
      do i = 1, n
         j = ORDC(i)
         Y(j) = X(i)
      enddo
*
      return
      end
*
***COMPLETE PIVOTING (Loop in the calling program)********************
*
      subroutine PIVOT(A,n,ORD,ORDC,i)
      implicit doubleprecision(a-h,o-z)
      parameter (nd = 40)
      doubleprecision A(nd,nd),ORD(nd),ORDC(nd)
      integer col,row,p,r
*
*     complete pivoting
      row = i
      col = i
      do p = i, n
      do r = i, n
         if (dabs(A(row,col)) < dabs(A(p,r))) then
            row = p
            col = r
         endif
      enddo
      enddo
*
      if (col.ne.i) then
         do ii = 1, n
            temp = A(ii,i)
            A(ii,i) = A(ii,col)
            A(ii,col) = temp
         enddo
         t = ORDC(i)
         ORDC(i) = ORDC(col)
         ORDC(col) = t
      endif
*
      m = n + 1 
      if (row.ne.i) then
         do jj = 1, m
            temp = A(i,jj)
            A(i,jj) = A(row,jj)
            A(row,jj) = temp
         enddo
         t = ORD(i)
         ORD(i) = ORD(row)
         ORD(row) = t
      endif
*     
      return
      end