fitting.f90

Go to the documentation of this file.

!> \file fitting.f90  Procedures to fit functions to data
 
 
!  Copyright (c) 2002-2025  Marc van der Sluys - Nikhef/Utrecht University - marc.vandersluys.nl
!   
!  This file is part of the libSUFR package, 
!  see: http://libsufr.sourceforge.net/
!   
!  This is free software: you can redistribute it and/or modify it under the terms of the European Union
!  Public Licence 1.2 (EUPL 1.2).  This software is distributed in the hope that it will be useful, but
!  WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
!  PURPOSE.  See the EU Public Licence for more details.  You should have received a copy of the European
!  Union Public Licence along with this code.  If not, see <https://www.eupl.eu/1.2/en/>.
!  
!  
 
 
 
 
!***********************************************************************************************************************************
!> \brief  Procedures to fit functions to data
 
module sufr_fitting
  implicit none
  save
  
  private
  
  public linear_fit_yerr, nonlin_fit_yerr, nonlin_fit_example_myfunc, basefunc_polynomial
  
contains
  
  !*********************************************************************************************************************************
  !> \brief  Fit the model coefficients using a chi^2 method
  !!
  !! Do a chi^2-minimisation fit of nDat data points xDat and yDat, the latter with errors yErr, to function func.  The fit uses
  !! the nCoef coefficients stored in fCoef() that should be varied if the corresponding entry in freeCoef=1, fixed if =0.
  !!
  !! The fitting function must be provided as the subroutine myFunc which returns a value y = myFunc(x, basefunc, nCoef), where
  !!   y is a linear combination of x, i.e.  y = SUM(i=1,nCoef) fCoef(i) * basefunc(i)(x).
  !!   For example, for a polynomial, basefunc(i) = x**(i-1).
  !!
  !! The routine returns the fit coefficients in fCoef(), the variance-covariance matrix in covMat sith size nCov, and the chi^2
  !!   value in chiSq.
  !!
  !!
  !! \param  nDat      Number of data points
  !! \param  xDat      X values of the data (nDat)
  !! \param  yDat      Y values of the data (nDat)
  !! \param  yErr      Errors (standard deviations) of the y values (nDat)
  !!
  !! \param  nCoef     Number of coefficients of the fitting function
  !! \param fCoef     Coefficients of the fitting function (nCoef) (output)
  !! \param  freeCoef  Fix coefficient coef(i) if freeCoef(i)=0, otherwise let if vary freely and fit it (nCoef)
  !!
  !! \param  nCov      Size of both dimensions of covMat()
  !! \param covMat    Covariance matrix (nCov,nCov) (output)
  !!
  !! \param chiSq     Chi^2 - chi squared (output)
  !!
  !! \param  myFunc    External subroutine that describes the model value of Y for given value X
  !!
  !! \note
  !! - Needs sort_var_covar_matrix() and solve_linear_equations_Gauss_Jordan()
  !!
  !! \see Numerical Recipes in Fortran 77, Sect.15.4
  !!
  
 subroutine linear_fit_yerr(nDat, xDat,yDat, yErr,  nCoef,fCoef,freeCoef,  nCov,covMat, chiSq, myFunc)
    use sufr_kinds, only: double
    use sufr_system, only: warn, quit_program_error
    implicit none
    integer, parameter :: mmax = 50
    integer, intent(in) :: ndat, ncoef,freecoef(ncoef), ncov
    real(double), intent(in) :: xdat(ndat),ydat(ndat), yerr(ndat)
    real(double), intent(inout) :: fcoef(ncoef)
    real(double), intent(out) :: chisq, covmat(ncov,ncov)
    
    integer :: ii,ij,ik,il,im, nfit
    real(double) :: err2i,bf_e2i, tot,ym, basefunc(mmax),beta(mmax)
    external myfunc
    
    if(minval(abs(yerr)).lt.10*tiny(yerr)) call quit_program_error('libSUFR linear_fit_yerr(): errors cannot be zero', 0)
    
    ! Determine number of parameters that need to be fit:
    nfit = 0
    do ij=1,ncoef
       if(freecoef(ij).ne.0) nfit = nfit+1
    end do
    if(nfit.eq.0) call warn('libSUFR linear_fit_yerr():  No parameters to be fit for', 0)
    
    covmat(1:nfit,1:nfit) = 0.d0
    beta(1:nfit) = 0.d0
    
    do ii=1,ndat
       call myfunc(xdat(ii),ncoef, basefunc)
       ym = ydat(ii)
       if(nfit.lt.ncoef) then
          do ij=1,ncoef
             if(freecoef(ij).eq.0) ym = ym - fcoef(ij)*basefunc(ij)
          end do
       end if
       err2i = 1.d0/yerr(ii)**2  ! 1/sigma_i^2
       
       ij = 0
       do il=1,ncoef
          if(freecoef(il).ne.0) then
             ij = ij + 1
             bf_e2i = basefunc(il)*err2i
             
             ik = 0
             do im=1,il
                if(freecoef(im).ne.0) then
                   ik = ik + 1
                   covmat(ij,ik) = covmat(ij,ik) + bf_e2i*basefunc(im)
                end if
             end do
             beta(ij) = beta(ij) + ym*bf_e2i
          end if
       end do
    end do
    
    do ij=2,nfit
       do ik=1,ij-1
          covmat(ik,ij) = covmat(ij,ik)
       end do
    end do
    
    ! Solve the linear algebraic equations using Gauss-Jordan elimination:
    call solve_linear_equations_gauss_jordan(nfit, ncov,covmat,  1, 1,beta)
    
    ij = 0
    do il=1,ncoef
       if(freecoef(il).ne.0) then
          ij = ij+1
          fcoef(il) = beta(ij)
       end if
    end do
    
    chisq = 0.d0
    do ii=1,ndat
       call myfunc(xdat(ii),ncoef, basefunc)
       
       tot = sum(fcoef(1:ncoef)*basefunc(1:ncoef))
       chisq = chisq  +  ( (ydat(ii)-tot) / yerr(ii) )**2
    end do
    
    ! Sort the resulting variance-covariance matrix:
    call sort_var_covar_matrix(ncov,covmat, ncoef,freecoef, nfit)
    
 subroutine linear_fit_yerr(nDat, xDat,yDat, yErr,  nCoef,fCoef,freeCoef,  nCov,covMat, chiSq, myFunc) …
  end subroutine linear_fit_yerr
  !*********************************************************************************************************************************
  
  
  
  
  !*********************************************************************************************************************************
  !> \brief  Levenberg-Marquardt method to reduce the value of chi-squared of a fit of a set of data points with errors in Y to a
  !!         non-linear function
  !!
  !! You will need to call this routine repeatedly, with different values for iterStat, until chiSq no longer decreases
  !!   (significantly)
  !!
  !! \param nY        Number of Y values for each X value (normally 1, but e.g. 2 for a sky position with two coordinates)
  !! \param nDat      Number of data points in X and Y
  !!
  !! \param xDat      X data points to fit
  !! \param yDat      Y data points to fit
  !! \param yErr      Errors (standard deviations) for yDat
  !!
  !! \param nCoef     Number of coefficients used to describe the non-linear function myFunc
  !! \param fCoef     Coefficients for the non-linear function myFunc, updated after each call
  !! \param freeCoef  Determines which coefficients for the non-linear function myFunc should be fitted for.  Set freeCoef(i) = 0 
  !!                    in order to keep fCoef(i) fixed
  !!
  !! \param nMat      Size of the matrices; > nCoef
  !! \param covMat    Variance-covariance matrix - returned on last call with iterStat.  Fixed parameters return zero (co)variances
  !! \param curvMat   Hessian/curvature matrix - double partial derivative of chi squared to two coefficients fCoef
  !!
  !! \param chiSq     Chi squared
  !! \param iterStat  Iteration status;  Set to <0 for initialisation in the first call, and let it vary in subsequent calls,
  !!                    until the fit converges.  After that, set to 0 to return the variance-covariance and curvature matrices
  !!                    on the final call.  IterStat decreases 10x if chiSq becomes smaller, and increases 10x otherwise.
  !! \param myFunc    External subroutine that describes the model value of Y and partial derivatives dY/dXi for given value X and
  !!                    function coefficients fCoef
  !!
  !! \see  Numerical Recipes in Fortran, 15.5: Modelling of Data / Non-linear models
  !! 
  !! \note  Uses sort_var_covar_matrix(), solve_linear_equations_Gauss_Jordan() and nonlin_fit_eval()
  
  subroutine nonlin_fit_yerr(nY,nDat, xDat,yDat, yErr,  nCoef,fCoef,freeCoef,  nMat,covMat,curvMat,  chiSq, iterStat, myFunc)
    use sufr_kinds, only: double
    use sufr_numerics, only: deq0
    
    implicit none
    integer, intent(in) :: ny,ndat, ncoef,freecoef(ncoef), nmat
    real(double), intent(in) :: xdat(ndat), ydat(ny,ndat), yerr(ndat)
    real(double), intent(inout) :: fcoef(ncoef), iterstat
    real(double), intent(out) :: curvmat(nmat,nmat), covmat(nmat,nmat), chisq
    
    integer, parameter :: mmax = 20
    integer :: ifit,mfit, icoef,cntcoef
    real(double) :: ochisq,trycoef(mmax),dchisq(mmax),dfcoef(mmax)
    
    save ochisq,trycoef,dchisq,dfcoef,mfit
    external myfunc
    
    
    ! An EXTERNAL subroutine must be CALLed - use a dummy call here, since the subroutine name is only passed on:
    if(.false.) call myfunc()  ! myFunc(nY, xDat(iDat), nCoef,fCoef, yMod,dyda)
    
    
    ! Initialisation on the first call:
    if(iterstat.lt.0.d0) then
       mfit = 0
       do icoef=1,ncoef
          if(freecoef(icoef).ne.0) mfit = mfit + 1
       end do
       
       iterstat = 0.001d0
       
       ! Evaluate the linearized curvature matrix and vector dChiSq/dfCoef, and compute chi squared:
       call nonlin_fit_eval(ny,ndat, xdat,ydat, yerr,  ncoef,fcoef,freecoef, nmat,curvmat,  dchisq,chisq, myfunc)
       
       ochisq = chisq
       trycoef(1:ncoef) = fcoef(1:ncoef)
    end if
    
    
    ! Augment the diagonal elements of the variance-covariance matrix:
    covmat(1:mfit,1:mfit) = curvmat(1:mfit,1:mfit)
    dfcoef(1:mfit)        = dchisq(1:mfit)
    
    do ifit=1,mfit
       covmat(ifit,ifit) = curvmat(ifit,ifit) * (1.d0+iterstat)
    end do
    
    
    ! Solve the linear algebraic equations using Gauss-Jordan elimination:
    call solve_linear_equations_gauss_jordan(mfit, nmat,covmat,  1, 1,dfcoef)
    
    
    ! Compute the variance-covariance and curvature matrices on the last call, and return:
    if(deq0(iterstat)) then
       call sort_var_covar_matrix(nmat, covmat,   ncoef, freecoef,  mfit)
       call sort_var_covar_matrix(nmat, curvmat,  ncoef, freecoef,  mfit)
       return
    end if
    
    
    ! Update the trial coefficients:
    cntcoef = 0
    do icoef=1,ncoef
       if(freecoef(icoef).ne.0) then
          cntcoef = cntcoef+1
          trycoef(icoef) = fcoef(icoef) + dfcoef(cntcoef)
       end if
    end do
    
    
    ! Evaluate the linearized curvature matrix and vector dChiSq/dfCoef, and compute chi squared:
    call nonlin_fit_eval(ny,ndat, xdat,ydat, yerr,  ncoef,trycoef,freecoef,  nmat,covmat,  dfcoef, chisq,  myfunc)
    
    
    ! If the new chi squared is better (smaller) than the old value, accept this solution and decrease iterStat...
    if(chisq.lt.ochisq) then
       iterstat = 0.1d0*iterstat
       ochisq = chisq
       
       curvmat(1:mfit,1:mfit) = covmat(1:mfit,1:mfit)
       dchisq(1:mfit) = dfcoef(1:mfit)
       
       fcoef(1:ncoef) = trycoef(1:ncoef)
       
    else                                        !  ...otherwise keep the old chi squared, and, if larger, increase iterStat:
       
       if(chisq.gt.ochisq) iterstat = 10*iterstat
       chisq = ochisq
       
    end if
    
  subroutine nonlin_fit_yerr(nY,nDat, xDat,yDat, yErr,  nCoef,fCoef,freeCoef,  nMat,covMat,curvMat,  chiSq, iterStat, myFunc) …
  end subroutine nonlin_fit_yerr
  !*********************************************************************************************************************************
  
  
  !*********************************************************************************************************************************
  !> \brief  Evaluate the linearized fitting matrix curvMat and vector dChiSq, and calculate the chi squared (chiSq)
  !!
  !! \param nY       Number of Y values for each X value (normally 1, but e.g. 2 for a sky position with two coordinates)
  !! \param nDat     Number of data points in X and Y
  !!
  !! \param xDat     X data points to fit
  !! \param yDat     Y data points to fit
  !! \param yErr     Errors (standard deviations) for yDat
  !!
  !! \param nCoef    Number of coefficients used to describe the non-linear function myFunc
  !! \param fCoef    Coefficients for the non-linear function myFunc, updated after each call
  !! \param freeCoef Determines which coefficients for the non-linear function myFunc should be fitted for.  Set freeCoef(i) = 0 in 
  !!                   order to keep fCoef(i) fixed
  !!
  !! \param nMat     Size of the matrices; > nCoef
  !! \param curvMat  Hessian/curvature matrix - double partial derivative of chi squared to two coefficients fCoef
  !! \param dChiSq   Vector containing the partial derivatives of chiSq to each coefficient in fCoef
  !!
  !! \param chiSq    Chi squared
  !! \param myFunc   External subroutine that describes the model value of Y and partial derivatives dY/dXi for given value X and
  !!                   function coefficients fCoef
  !!
  !!    Since the 1D chi squared is the sum of the squares of the weighted differences,
  !!    \f$\chi^2 = \sum_j d_j^2\f$,  and the nY-D distance \f$d_j\f$ is defined as \f$d_j = \sqrt{\sum_i d_{ij}^2}\f$,
  !!    where \f$d_i\f$ is the distance in each dimension, the nY-D chi squared is the sum of
  !!    the squares of these distances, i.e.
  !!    \f$\chi^2 = \sum_j d_j^2  =  \sum_j \sum_i d_{ij}^2 = \sum_i \sum_j d_{ij}^2 = \sum_i \chi^2_i\f$,
  !!    simply the sum of the 1D chi squared values.  The same holds for the derivatives of
  !!    \f$\chi^2\f$, where the multidimensional derivative is the sum of the 1D values.
  !!
  !! \see  Numerical Recipes in Fortran, 15.5: Modelling of Data / Non-linear models
  
  subroutine nonlin_fit_eval(nY,nDat, xDat,yDat, yErr,  nCoef,fCoef,freeCoef, nMat,curvMat,  dChiSq,chiSq, myFunc)
    use sufr_kinds, only: double
    implicit none
    integer, intent(in) :: ny,ndat, ncoef,freecoef(ncoef), nmat
    real(double), intent(in) :: xdat(ndat),ydat(ny,ndat),yerr(ndat), fcoef(ncoef)
    real(double), intent(out) :: chisq,curvmat(nmat,nmat),dchisq(ncoef)
    
    integer, parameter :: mmax = 20
    integer :: ifit,jfit,mfit, iy,idat, icoef,jcoef, cntcoefi,cntcoefj
    real(double) :: err2i,dyda_e2i, dy,ymod(ny),dyda(ny,mmax)
    external myfunc
    
    
    ! Determine the number of non-fixed parameters:
    mfit = 0
    do icoef=1,ncoef
       if(freecoef(icoef).ne.0) mfit = mfit + 1
    end do
    
    
    ! Initialise the chi squared, the chiSq-derivative vector and the curvature matrix:
    chisq   = 0.d0
    dchisq  = 0.d0
    curvmat = 0.d0
    
    
    ! Sum over all data points:
    do idat=1,ndat
       call myfunc(ny, xdat(idat),  ncoef,fcoef,  ymod,dyda)
       
       err2i = 1.d0/(yerr(idat)**2)  ! 1/sigma^2
       
       do iy=1,ny  ! Loop over multiple Y variables (if nY>1)
          dy = ydat(iy,idat) - ymod(iy)
          
          cntcoefi = 0
          do icoef=1,ncoef
             if(freecoef(icoef).ne.0) then
                cntcoefi = cntcoefi + 1
                dyda_e2i = dyda(iy,icoef)*err2i  ! dy/dfCoef_i / sigma^2
                
                cntcoefj = 0
                do jcoef=1,icoef
                   if(freecoef(jcoef).ne.0) then
                      cntcoefj = cntcoefj + 1
                      curvmat(cntcoefi,cntcoefj) = curvmat(cntcoefi,cntcoefj) + &
                           dyda_e2i*dyda(iy,jcoef) ! Sum dy/dfCoef_i dy/dfCoef_j / sig^2
                   end if
                end do  ! jCoef
                
                dchisq(cntcoefi) = dchisq(cntcoefi) + dy*dyda_e2i
             end if
          end do  ! iCoef
          
          chisq = chisq + dy**2 * err2i  ! Compute chi squared
          
       end do  ! iY
       
    end do  ! iDat
    
    
    ! Fill in the symmetric side of the curvature matrix:
    do ifit=2,mfit
       do jfit=1,ifit-1
          curvmat(jfit,ifit) = curvmat(ifit,jfit)
       end do  ! jFit
    end do  ! iFit
    
  end subroutine nonlin_fit_eval
  !*********************************************************************************************************************************
  
  
  !*********************************************************************************************************************************
  !> \brief  Dummy example function myFunc for nonlin_fit_yerr(): return the value and partial derivatives
  !!
  !! \param nY     Number of Y values for each X value (normally 1, but e.g. 2 for a sky position with two coordinates)
  !! \param xDat   Input X value for the current data point
  !!
  !! \param nCoef  Number of coefficients
  !! \param fCoef  Vector of coefficients that describe the function
  !!
  !! \param yDat  Y values for the current data point (output)
  !! \param dyda  Partial derivatives for yDat:  1: dy/dfCoef(1),  ...,  n: dy/dfCoef(n) (output)
  !!
  !!
  !! \note  Write a subroutine with the same interface to use with nonlin_fit_yerr()
  
  subroutine nonlin_fit_example_myfunc(nY, xDat,  nCoef,fCoef,  yDat,dyda)
    use sufr_kinds, only: double
    use sufr_system, only: quit_program_error
    
    implicit none
    integer, intent(in) :: ny, ncoef
    real(double), intent(in) :: xdat, fcoef(ncoef)
    real(double), intent(out) :: ydat(ny), dyda(ny,ncoef)
    
    if(ny.ne.1)    call quit_program_error('nonlin_fit_eval_example_myFunc():  nY must be equal to 1 ',    1)
    if(ncoef.ne.3) call quit_program_error('nonlin_fit_eval_example_myFunc():  nCoef must be equal to 3 ', 1)
    
    ydat(1)   = fcoef(1)*xdat**2 + fcoef(2)   ! Replace with desired function
    dyda(1,1) = fcoef(2) * xdat               ! Replace with partial derivative w.r.t. first variable - dyDat/dfCoef(1)
    ! ...
    dyda(1,ncoef) = fcoef(1)                  ! Replace with partial derivative w.r.t. n-th variable - dyDat/dfCoef(nCoef)
    
  subroutine nonlin_fit_example_myfunc(nY, xDat,  nCoef,fCoef,  yDat,dyda) …
  end subroutine nonlin_fit_example_myfunc
  !*********************************************************************************************************************************
  
  
  !*********************************************************************************************************************************
  !> \brief  Sort covariance matrix returned by linear_fit_yerr() and nonlin_fit_yerr() to true order of fitting coefficients
  !!
  !! \param  nCov      Size of both dimensions of covMat()
  !! \param  covMat    Variance-covariance matrix
  !!
  !! \param  nCoef     Number of coefficients of the fitting function
  !! \param  freeCoef  Fix coefficient coef(i) if freeCoef(i)=0, otherwise fit for it
  !!
  !! \param  nFit      Number of parameters that need to be fit for
  !!
  !! \see Numerical Recipes in Fortran 77, Sect.15.4
  !!
  
  pure subroutine sort_var_covar_matrix(nCov,covMat, nCoef,freeCoef, nFit)
    use sufr_kinds, only: double
    use sufr_system, only: swapdbl
    
    implicit none
    integer, intent(in) :: ncov,ncoef, freecoef(ncoef),nfit
    real(double), intent(inout) :: covmat(ncov,ncov)
    
    integer :: ii,ij,ik
    
    do ii=nfit+1,ncoef
       do ij=1,ii
          covmat(ii,ij) = 0.d0
          covmat(ij,ii) = 0.d0
       end do
    end do
    
    ik = nfit
    
    do ij=ncoef,1,-1
       if(freecoef(ij).ne.0) then
          do ii=1,ncoef
             call swapdbl(covmat(ii,ij), covmat(ii,ik))  ! Swap the two numbers
          end do
          
          do ii=1,ncoef
             call swapdbl(covmat(ij,ii), covmat(ik,ii))  ! Swap the two numbers
          end do
          
          ik = ik-1
       end if
    end do
    
  end subroutine sort_var_covar_matrix
  !*********************************************************************************************************************************
  
  
  
  !*********************************************************************************************************************************
  !> \brief  Solve a set of linear algebraic equations using Gauss-Jordan elimination.
  !!
  !! \param nMat     Dimension of the equation set
  !! \param nMatArr  Physical dimensions of the array matArr  (nMatArr x nMatArr)
  !! \param matArr   Input: array containing the matrix to be inverted - output: array containing the inverted matrix  (I/O)
  !!
  !! \param nVec     Number of vectors in vecArr
  !! \param nVecArr  One of the physical dimensions of the array vecArr  (nMatArr x nVecArr)
  !! \param vecArr   Input: array containing nVec right-hand side vectors - output: solution vectors  (I/O)
  !!
  !! \see Numerical Recipes in Fortran 77, Sect.2.1
  !!
  
  subroutine solve_linear_equations_gauss_jordan(nMat, nMatArr,matArr,  nVec, nVecArr,vecArr)
    use sufr_kinds, only: double
    use sufr_system, only: warn, swapdbl
    
    implicit none
    integer, intent(in) :: nmat,nmatarr, nvec,nvecarr
    real(double), intent(inout) :: matarr(nmatarr,nmatarr),vecarr(nmatarr,nvecarr)
    
    integer, parameter :: nmax = 50  ! Maximum allowed value for nMat
    integer :: ii,icol,irow,ij,ik,il,ll,indxc(nmax),indxr(nmax),ipiv(nmax)
    real(double) :: arrmax,tmpdbl,pivinv
    
    
    icol = 0
    irow = 0
    ipiv(1:nmat) = 0
    
    
    do ii=1,nmat
       ! Find maximum array value:
       arrmax = 0.d0
       do ij=1,nmat
          if(ipiv(ij).ne.1) then
             do ik=1,nmat
                if(ipiv(ik).eq.0) then
                   if(abs(matarr(ij,ik)).ge.arrmax) then
                      arrmax = abs(matarr(ij,ik))
                      irow = ij
                      icol = ik
                   end if
                else if(ipiv(ik).gt.1) then
                   call warn('libSUFR solve_linear_equations_Gauss_Jordan():  singular matrix',0)
                end if
             end do
          end if
       end do
       
       ipiv(icol) = ipiv(icol) + 1
       
       if(irow.ne.icol) then
          do il=1,nmat
             call swapdbl(matarr(irow,il), matarr(icol,il))
          end do
          
          do il=1,nvec
             call swapdbl(vecarr(irow,il), vecarr(icol,il))
          end do
       end if
       
       indxr(ii) = irow
       indxc(ii) = icol
       
       if(abs(matarr(icol,icol)).lt.10*tiny(matarr)) call warn('libSUFR solve_linear_equations_Gauss_Jordan():  singular matrix', 0)
       
       pivinv = 1.d0/matarr(icol,icol)
       matarr(icol,icol) = 1.d0
       
       matarr(icol,1:nmat) = matarr(icol,1:nmat)*pivinv
       vecarr(icol,1:nvec) = vecarr(icol,1:nvec)*pivinv
       
       do ll=1,nmat
          if(ll.ne.icol) then
             tmpdbl = matarr(ll,icol)
             matarr(ll,icol) = 0.d0
             matarr(ll,1:nmat) = matarr(ll,1:nmat) - matarr(icol,1:nmat)*tmpdbl
             vecarr(ll,1:nvec) = vecarr(ll,1:nvec) - vecarr(icol,1:nvec)*tmpdbl
          end if
       end do
    end do
    
    do il=nmat,1,-1
       if(indxr(il).ne.indxc(il)) then
          do ik=1,nmat
             call swapdbl(matarr(ik,indxr(il)), matarr(ik,indxc(il)))
          end do
       end if
    end do
    
  end subroutine solve_linear_equations_gauss_jordan
  !*********************************************************************************************************************************
  
  
  
  
  !*********************************************************************************************************************************
  !> \brief  User-provided base function to fit a polynomial to
  !!
  !! \param  xVal      X value to evaluate yVal = basefunc_polynimial(xVal) for
  !! \param  nCoef     Number of coefficients for the base function
  !!
  !! \param baseFunc  Base function i for coefficent i.  For a polynomial, BF(i) = x**(i-1) (output)
  
  pure subroutine basefunc_polynomial(xVal,nCoef, baseFunc)
    use sufr_kinds, only: double
    implicit none
    integer, intent(in) :: ncoef
    real(double), intent(in) :: xval
    real(double), intent(out) :: basefunc(ncoef)
    
    integer :: coefi
    
    basefunc(1) = 1.d0
    do coefi = 2,ncoef
       basefunc(coefi) = basefunc(coefi-1) * xval
    end do
    
  pure subroutine basefunc_polynomial(xVal,nCoef, baseFunc) …
  end subroutine basefunc_polynomial
  !*********************************************************************************************************************************
  
  
  
  
end module sufr_fitting
!***********************************************************************************************************************************