! lsfit_lect.f90  from lecture notes
!               compile  gfortran -o lsfit_lect lsfit_lect.f90 simeq.f90 -lm
!               run      lsfit_lect > lsfit_lect_f90.out

program lsfit_lect
  implicit none
  double precision, dimension(3,6) :: xd = reshape((/1.0, 2.5, 3.7, &
                                                     2.0, 2.5, 3.6, &
                                                     3.0, 2.7, 3.5, &
	 	                                     2.2, 2.1, 3.1, &
		                                     1.5, 2.0, 2.6, &
		                                     1.6, 2.0, 3.1/),(/3,6/))
  double precision, dimension(6) :: yd = (/32.5, 7.2, 6.9, 22.4, 10.4, 11.3/) ! Y actual
  double precision, dimension(3,3) :: A
  double precision, dimension(3) :: X
  double precision, dimension(3) :: Y ! |A| * |X| = |Y|  |X| will be  a, b, c
  integer :: n = 3 ! number of terms in fit equation
  integer :: i, j, k
  interface
    subroutine simeq(n, sz, A, Y, X)
      implicit none
      integer, intent(in) :: n  ! number of equations 
      integer, intent(in) :: sz ! dimension of arrays
      double precision, dimension(sz,sz), intent(in) :: A
      double precision, dimension(sz), intent(in) :: Y
      double precision, dimension(sz), intent(out) :: X
    end subroutine simeq
  end interface
  
  print *, "lsfit_lect.f90 run to make lsfit_lect_f90.out"
  do i=1,n
    Y(i) = 0.0 ! must start sum with zero
    do k=1,6
      Y(i) = Y(i) + xd(i,k)*yd(k) ! Xi * Yk
    end do
    do j=1,3
      A(i,j) = 0.0  ! must start sum with zero
      do k=1,6 ! run through observations
        A(i,j) = A(i,j) + xd(j,k)*xd(i,k)  ! Xj * Xi
      end do ! k
    end do ! j
  end do  ! i
  call simeq(n, n, A, Y, X)
  print *, "a=",X(1),", b=",X(2),", c=",X(3)
  ! check
  print *, "really bad fit"
  do k=1,6  ! k is observation
    print *, "Y_actual=",Yd(k),", Y_approximate=", &
             X(1)*xd(1,k)+X(2)*xd(2,k)+X(3)*xd(3,k)
  end do

end program lsfit_lect