!! Barlow and Barnett, Computing for Scientists 1998
!! page 223                 C++ and Fortran 90

PROGRAM LS_fit_poly

  IMPLICIT NONE

  !  First construct the normal matrix = k x k from A = N x k
  !  part of programs to solve least-squares polynomial fits to data
  !  N = number of data points , k-1 = polynomial order  
  !  y = data at x

  INTEGER :: i , j , k , N
  PARAMETER(k=3 , N=8)     ! the user could read these in before allocation
  REAL(KIND=4) , DIMENSION(: , :) , ALLOCATABLE :: A , Nmat , L
  REAL(KIND=4) , DIMENSION(:) ,     ALLOCATABLE :: x, y, rhs, p, q 
  CHARACTER :: fmt*20
  LOGICAL   :: details = .FALSE.
  ALLOCATE ( A(N,k) , Nmat(k,k), rhs(k) , x(N ) , y(N))
   
  x = (/40., 80., 120., 160., 200., 240., 280., 320. /)    !==thermocouple deg C
  y = (/0.38, 0.64, 0.84, 0.96, 0.98, 0.93, 0.80, 0.61/)   !==thermocouple voltage
  fmt = '(3f20.4)'
  x = x / 100.0      ! scale the temperature vector
  
  A(: , 1) = 1.0;  A(: , 2) = x;  A(: , 3) = x*x        ! build up A by columns
  
  IF (details) THEN
     PRINT fmt , A(:,1); PRINT * , '/A(:,1)'
     PRINT fmt , A ;     PRINT * , '/A '
  ENDIF

  Nmat = MATMUL(TRANSPOSE(A) , A)
  rhs  = MATMUL(TRANSPOSE(A) , y)

  IF (details) THEN
     PRINT fmt , Nmat;   PRINT * , '/Nmat '
     PRINT fmt , rhs;    PRINT * , '/rhs '
  ENDIF

  ! ...
  !  Next part of program to  solve least-squares polynomial fits to data
  !  Form Cholesky lower triangle matrix  L from normal matrix V = A'A 

  ALLOCATE(  L(k , k) );    L = 0.0      ! essential to zero the new array
  IF (details) THEN
     PRINT '(A,/9E12.4)' , 'L= 0.0' , L
  ENDIF
  DO i=1,k
     L(i,i) = SQRT(Nmat(i,i) - DOT_PRODUCT(L(i , 1:i-1) , L(i , 1:i-1)))
     L( i+1:k , i) = (Nmat(i , i+1:k) - MATMUL(L(i+1:k , 1:i-1) , L(i , 1:i-1)))
     L(i+1:k ,i) = L(i+1:k , i) / L(i,i)
       IF (details) THEN
          PRINT fmt , TRANSPOSE(L); PRINT * , " "
       ENDIF
  ENDDO
  IF (details) THEN
     PRINT * , " /Press L"; PRINT fmt , TRANSPOSE(L)
  ENDIF


  ALLOCATE( q(k ) ) ; q = 0.0      ! solve for q using  Lq = rhs
  IF (details) PRINT * , "1"
  DO j=1,k
     q(j)=( rhs(j) - DOT_PRODUCT(L(j , 1:j-1), q(1:j-1) )  )/ L(j,j)
  ENDDO
  IF (details) THEN
     PRINT fmt , q; PRINT * , '/q '
  ENDIF

  ALLOCATE( p(k ) ); p = 0.0        ! solve for p using  L'p = q
  IF (details) PRINT * , "2"
  DO j = k , 1 , -1
     p(j)=( q(j) - DOT_PRODUCT(L(J+1:k , j), p(j+1:k) )  )/ L(j,j)
  ENDDO
  IF (details) THEN
     PRINT fmt, p; PRINT * , '/p'
  ENDIF

  PRINT * , ' Fit to the x , y data vectors '
  PRINT * , '       x :      ' , ' y :'
  PRINT '(1X , F10.2 , F10.3)' , ( x(j) , y(j) ,j=1,N )
  PRINT * , ' solution to y = a + bx + cx*x is ' 
  PRINT '(3(A , F12.5/))' , ' a = ',  p(1)  , ' b = ' , p(2) ,  ' c = ' , p(3)
 
  DEALLOCATE( A , Nmat ,   L , q )
 
  !...

END PROGRAM LS_fit_poly
