Programming Example 13

Program Least_squares
implicit none

! Global variables
integer, parameter:: m=4       ! Number of coefficient
real(8):: S(m)                 ! Array of coefficient S
real(8), parameter:: TOL=1D-6  ! Numerical tolerance
real(8), dimension(:), allocatable:: X, Y  ! Arrays of data X and Y

Call Initialization(X, Y, S)  ! read data from file
Call Processing(X, Y, S)      ! perform numerical analysis and display
Call Finalization(X, Y, S)    ! write data to file

contains

!================================================================================
Subroutine Initialization(X, Y, S)
! dummy variables
real(8), dimension(:), allocatable, INTENT(OUT):: X, Y
real(8), INTENT(OUT):: S(m)
! local variables
integer:: i, n, status=0
character(len=50):: filename
real(8):: r

write(*,'(A)') "Enter the title of input data with the file extension"
read(*,*) filename

open(1, file=filename)
n=0
do
  read(1, *, iostat=status) r
  if (status/=0) exit
  n=n+1
end do
allocate(X(n), Y(n))
rewind(1)
do i=1,n
   read(1, *) X(i), Y(i)
end do
close(1)

write(*,'(A)') "The input data"
write(*,'(2A10)') "Data X", "Data Y"
do i=1,n
   write(*,'(2F10.3)') X(i), Y(i)
end do

write(*,*) "The fitting function is Y = exp(-S1*X)*sin(S2*X)*S3+S4"
write(*,*) "In order to determine non-trivial fitting coefficients, intialial values are necessary."
do i = 1, m
 write(*, '(A, I1, A, $)') "Input the guessed value for S", i, " = "
 Read(*,*) S(i)
end do

end subroutine Initialization

!================================================================================
Subroutine Processing(X, Y, S)
! dummy variables
real(8), dimension(:), INTENT(IN):: X, Y
real(8), INTENT(INOUT):: S(m)
! local variables
integer:: i, k, k_max=50 ! iteration number and maximum iteration number
! Norm of errors, vector residuals and Hessian matrix
real(8):: err, R(m), H(m,m)

write(*, '(A)') repeat('-',10*m+30)
write(*,'(A10, $)') "Iteration"
do i=1, m
 write(*,'(A9, I1, $)') "S", i
end do
write(*,'(A20)') "Norm of residuals"
write(*, '(A)') repeat('-',10*m+30)

R=Residual_vector(X, Y, S)
err=sqrt(dot_product(R,R))
k=0
write(*,'(I10, $)') k
do i=1, m
  write(*,'(F10.3, $)') S(i)
  end do
write(*,'(E20.3)') ERR
do while *1
  k=k + 1
  H=Hessian_matrix(X, Y, S) !This caluclates the Hessian wrt X,Y,S
  S=S-matmul(Minverse(H),R) !This remakes S. S=S- [Hessian]**(-1) ×　RessidualMatrix
  R=Residual_vector(X, Y, S) !Remakes R with new S.
  err=sqrt(dot_product(R,R)) !Remakes ERR with new ResidualVector
  write(*,'(I10, $)') k
  do i=1, m
    write(*,'(F10.3, $)') S(i)
    end do
  write(*,'(E20.3)') ERR
end do
if (k==k_max) then
   write(*,*) "Divergence!"; STOP
end if
write(*, '(A)') repeat('-',10*m+30)

end Subroutine Processing

!================================================================================
Subroutine Finalization(X, Y, S)
! dummy variables
real(8), dimension(:), INTENT(IN):: X, Y
real(8), INTENT(IN):: S(m)
! local variables
Integer:: i, n
character(len=50):: filename

n=size(X) ! rows of vector X
write(*,'(A)') "Enter the title of output data with the file extension"
read(*,*) filename

open(2, file=filename)
write(*,'(A)') "The output data"
write(*,'(2A10)') "Data X", "Fitting Y"
do i=1,n
   write(*,'(2F10.3)') X(i), Fit(X(i),S)
   write(2,'(2F10.3)') X(i), Fit(X(i),S)
end do

write(*, '(A, F10.3)') "Sum square of errors of the converged result is ", SSE(X, Y, S)
write(*,'(A, A, $)') "The output data have been successfully written to the file ", filename
close(2)

end Subroutine Finalization

!--------------------------------------------------------------------------------
Function SSE(X, Y, C) ! Sum Square Errors
real(8), dimension(:):: X, Y
real(8):: SSE, C(m) ! sum square of residual and coefficient array
integer:: i, n

n = size(X)    ! The size of array X
SSE = 0
do i = 1, n
   SSE = SSE + ( fit(X(i),C) - Y(i) )**2
end do

end Function SSE

!--------------------------------------------------------------------------------
Function Fit(z, C)
real(8):: Fit, z, C(m)

Fit = exp(-C(1)*z)*sin(C(2)*z)*C(3)+C(4)

end Function Fit

!--------------------------------------------------------------------------------
Function Residual_vector(X, Y, C)
real(8), dimension(:):: X, Y
real(8):: Residual_vector(m), C(m), C_(m)
integer:: i

do i = 1, m
   C_ = C; C_(i) = C(i) + TOL
   Residual_vector(i) = ( SSE(X, Y, C_) - SSE(X, Y, C) ) / TOL
end do
end Function Residual_vector

!--------------------------------------------------------------------------------
Function Hessian_matrix(X, Y, C)
real(8), dimension(:):: X, Y
real(8):: Hessian_matrix(m,m), C(m), C_(m), v(m), v_(m)
integer:: i,j

do i = 1, m
   do j = 1, m
      C_ = C; C_(j) = C(j) + TOL
      v_ = Residual_vector(X, Y, C_)
      v  = Residual_vector(X, Y, C)
      Hessian_matrix(i,j) = ( v_(i) - v(i) ) / TOL
   end do
end do

end Function Hessian_matrix

!--------------------------------------------------------------------------------
Function Minverse(A)
real(8), dimension(:,:):: A ! Allow both static/dynamic memory allocation
real(8), dimension(:,:), allocatable:: Minverse, B, V
real(8), dimension(:), allocatable:: temp
real(8), parameter:: TOL=1D-6
real(8):: pivot
integer:: i, j, jmax, k, n

n=size(A)**0.5
allocate(Minverse(n,n), B(n,n), V(n,n), temp(n))
B=A
V(:,:)=0.
do i=1, n; V(i,i)=1.; end do

do i=1, n
   jmax=i

!===== Maximum pivoting technique =====

   do j=i+1, n
      if (abs(B(j,i))>abs(B(jmax,i))) then
         jmax=j
         temp=B(i,:); B(i,:)=B(jmax,:); B(jmax,:)=temp
         temp=V(i,:); V(i,:)=V(jmax,:); V(jmax,:)=temp
       end if
   end do

!===== Gauss-Jordan elimination =====

   pivot=B(i,i)
   if (abs(pivot)<TOL) then
      print *, "The matrix is singular!"
      stop
   end if
   B(i,:)=B(i,:)/pivot; V(i,:)=V(i,:)/pivot

   do k=1, n
      pivot=B(k,i)
      if (k/=i) then
         B(k,:)=B(k,:)-pivot*B(i,:)
         V(k,:)=V(k,:)-pivot*V(i,:)
      end if
   end do

end do

Minverse=V

end function Minverse

End Program

*1:err>TOL).and.(k<k_max