# UCL WIKI

##### Child pages
• fitting.f90: Using LAPACK routines to interpolate a 1D function by polynomials
Go to start of banner

# fitting.f90: Using LAPACK routines to interpolate a 1D function by polynomials

Using Code::Blocks walk through the source code and check if you understand all Fortran language structures.

Exercise 1: Compile and run the program using the fitting3.dat (see also fitting.dat, fitting2.dat),  data set as an input file, setting the maximum polynomial degree to 10 and the reference expansion point to 1.4.
Exercise 2: Improve the performance of the part where the linear equations matrix and vector are composed by pre-computing the intermediates "x_data(1:npoints)-x0)**ideg" prior the nested loops.
Exercise 3: OMP-parallelize loops computing the linear equations matrix and vector.

```! Least-squares fit of a univariate data by a polynomial.
!
! Go through the source code and check if you understand all Fortran language structures.
!
! Exercise1: compile and run the program for "fitting3.dat" data set file, setting maximum polynomial degree to 10 and the reference expansion point to 1.4
!
! Exercise2: improve performance of the program at a part where linear equations matrix and vector are composed
!            by pre-computing the intermediates "x_data(1:npoints)-x0)**ideg" prior the nested loops.
!
! Exercise3: OMP-parallelize loops computing the linear equations matrix and vector.

program fitting
implicit none
integer :: iounit, info, ipoint, npoints, maxdeg, ideg, jdeg, lwork, nsv
real(8) x, y, x0, rms, dev, svtol
real(8), allocatable :: x_data(:), y_data(:), matrix(:,:), vector(:), work(:), sv(:)
character(80) filename

print*, 'Enter name of data file'
print*, 'Enter max polynomial degree'
print*, 'Enter expansion point (x0)'

!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!
write(*,'(/1x,a,1x,a)') 'read data points from file:', trim(filename)
! open data file
iounit = 10
open(iounit, form='formatted', file=filename, action='read', status='old', iostat=info)
if (info/=0) then
stop
endif
! first we estimate the number of data points by reading till the EOF
ipoint = 0
do while(.true.)
if (info/=0) exit   ! exit loop if EOF
ipoint = ipoint + 1
enddo
npoints = ipoint
write(*,'(1x,a,1x,i3)') 'number of data points:', npoints
! allocate arrays to store data points
allocate(x_data(npoints), y_data(npoints), stat=info)
! check if allocation was successful
if (info/=0) stop 'stop: failed to allocate arrays x_data, y_data'
rewind(iounit)  ! next "read(iounit,..)" will start from the beginning of the file
x_data = 0.0
y_data = 0.0
ipoint = 0
do while(.true.)
if (info/=0) exit
ipoint = ipoint + 1
x_data(ipoint) = x
y_data(ipoint) = y
enddo
npoints = ipoint
! close file
close(iounit)
! print data on screen
write(*,'(1x,a)') 'data points'
do ipoint=1, npoints
write(*,'( 1x, i4, 3x, f12.6, 3x, f12.6 )') ipoint, x_data(ipoint), y_data(ipoint)
enddo

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!   form linear least-squares system   !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
write(*,'(/1x,a,1x,i3,1x,a,1x,i3,1x,a)') 'compose', maxdeg+1, 'by', maxdeg+1, 'linear system'
! allocate arrays to keep linear system matrix and right-hand-side vector
allocate(matrix(0:maxdeg,0:maxdeg), vector(0:maxdeg), stat=info)
if (info/=0) stop 'stop: failed to allocate matrix and vector'
matrix = 0.0
vector = 0.0

! ---- start the most computationally expensive part ---- !
!          need to be optimized and parallelized
! compute linear system matrix and right-hand-side vector
do ideg=0, maxdeg
do jdeg=0, ideg
matrix(jdeg,ideg) = sum( (x_data(1:npoints)-x0)**ideg * (x_data(1:npoints)-x0)**jdeg )
if (jdeg/=ideg) matrix(ideg,jdeg) = matrix(jdeg,ideg)
enddo
enddo
do ideg=0, maxdeg
vector(ideg) = sum( (x_data(1:npoints)-x0)**ideg * y_data(1:npoints) )
enddo
! --------------------------------------------------------!

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!   solve linear equations system using LAPACK library   !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
write(*,'(/a)') 'solve linear system by singular value decomposition'
! solve system by singular value decomposition
lwork = (maxdeg+1)**2*32
svtol = -1.0d-12
allocate(work(lwork), sv(0:maxdeg), stat=info)
if (info/=0) stop 'stop: failed to allocate work and sv arrays'
! use linear equations solver from Lapack/MKL (documentation https://software.intel.com/en-us/node/521114)
call dgelss( maxdeg+1, maxdeg+1, 1, matrix(0:maxdeg,0:maxdeg), maxdeg+1, vector(0:maxdeg), maxdeg+1, sv(0:maxdeg), svtol, nsv, work, lwork, info )
if (info/=0) stop 'error: SVD failed'
! compute covariance (degree of correlation between two linear solutions or measure of linear dependence between two solutions)
! COVAR = V^T * S^{-1} * V,
!   where "V" is right singular matrix,
!   and "S" is diagonal matrix with singular values on diagonal
! allocate arrays to store covariance matrix and some intermediates
!allocate(covar(0:maxdeg,0:maxdeg), w()..., ... stat=info)
!if (info/=0) ...
! 1. compute w=S^{-1}
!w = 0.0
!forall (i=0:maxdeg) ...
! 2. compute S^{-1} * V
!call dgemm(...)
! 3. compute V^T * S^{-1} * V
!call dgemm(...)
! print values of the polynomial coefficients and diagonal elements of covariance matrix
write(*,'(/a)') 'resulting coefficients'
do ideg=0, maxdeg
write(*,'(1x,i3,1x,f16.8,1x,f16.8)') ideg, vector(ideg) !, covar(ideg,ideg)
enddo

! compute and print deviations between model polynomial and data points
write(*,'(/a)') 'deviations'
rms = 0.0
do ipoint=1, npoints
y = sum( (/( vector(ideg) * (x_data(ipoint)-x0)**ideg, ideg=0, maxdeg )/) )
dev = y_data(ipoint) - y
rms = rms + dev**2
!write(*,'( 1x, f12.6, 1x, f12.6, 1x, f12.6, 1x, f12.6 )') x_data(ipoint), y_data(ipoint), y, dev
enddo
write(*,'(/a,1x,f12.8)') 'root-mean-square deviation:', sqrt(rms/real(npoints,kind=8))

! deallocate all allocatable arrays
deallocate(matrix, vector, x_data, y_data, work, sv)
end program fitting

```
• No labels