UCL WIKI

UCL Logo
Child pages
  • fitting.f90: Using LAPACK routines to interpolate a 1D function by polynomials
Skip to end of metadata
Go to start of metadata

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'
  read*, filename
  print*, 'Enter max polynomial degree'
  read*, maxdeg
  print*, 'Enter expansion point (x0)'
  read*, x0

  !!!!!!!!!!!!!!!!!!!!!!!
  !   read data points  !
  !!!!!!!!!!!!!!!!!!!!!!!
  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
    write(*,*) 'data file', trim(filename), 'is not found'
    stop
  endif
  ! first we estimate the number of data points by reading till the EOF
  ipoint = 0
  do while(.true.)
    read(iounit,*,iostat=info) x, y
    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
  ! read data points
  x_data = 0.0
  y_data = 0.0
  ipoint = 0
  do while(.true.)
    read(iounit,*,iostat=info) x, y
    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