Module helixfit_fitllsq_mod

module helixfit_fitllsq_mod

        ! Uses
    use precision_mod
    use namelist_mod
    use track_mod
    use helixfit_common_mod
    use helixfit_fitpar_mod
    use helixfit_hits_mod
    use helixfit_track_mod
    use helixfit_tref_mod

        ! Types
    private type gradient_type

        ! Variables
    integer, private, PARAMETER :: kMaxGradPoints = kMaxPar + 1
    integer, private, EXTERNAL :: solve_linear
    logical, private, dimension (kMaxPar) :: haveUnitsimplex = .false.
    real (kind=R8), private, dimension (kMaxPar,kMaxPar+1,kMaxPar) :: unitsimplex

        ! Subroutines and functions
    private subroutine SetupSimplex (numpar)
    public function HelixFit2gaussNewton (hits, kinks, lastplane, fitparMC, fitpar0, restype, maxiter, parstep, ambig0, doFitTref, res0, chi0, fitparB, resB, chiB, numiter) result (ierror)
    private subroutine AngleTest0 (grad, maxcosa, iexcl)
    private subroutine AngleTest (grad, maxcosa)
    private function iterate_llsq (hits, kinks, lastplane, fitparMC, fitpar0, maxiter, doStopChebyshev, doFitTref, cutoff0, fitpar, res, grad, chi, numiter, costrms, remove2) result (ierror)
    private function llsqstep (grad, res, dpar) result (ierror)
    private function addGradient (grad, res, chi) result (ierror)
    private function setupGradient (hits, kinks, lastplane, fitpar0, res0, parstep, grad) result (ierror)
    private function fitGradient (hits, fitpar, grad) result (ierror)
    private subroutine printGradient (grad)
    private subroutine AmbigCmpMC (res, cost)
    private function sqr (a) result (b)
    private function myReal8 (a4) result (a8)

end module helixfit_fitllsq_mod

Description of Types

private type gradient_type
    type (residuals_type) :: base
    integer :: numPoints
    type (residuals_type), dimension (kMaxGradPoints) :: points
    integer :: numResUsed
    real (kind=R8), dimension (kMaxPar) :: parSize
    real (kind=R8) :: parSizeMax
    real (kind=R8), dimension (kMaxPar*kMaxResiduals) :: cgrad
    logical, dimension (kMaxResiduals) :: enabled
end type gradient_type

Description of Variables

integer, private, PARAMETER :: kMaxGradPoints = kMaxPar + 1

integer, private, EXTERNAL :: solve_linear
 external solve linear equations

logical, private, dimension (kMaxPar) :: haveUnitsimplex = .false.
 Simplex template

real (kind=R8), private, dimension (kMaxPar,kMaxPar+1,kMaxPar) :: unitsimplex

Description of Subroutines and Functions

private subroutine SetupSimplex (numpar)
    integer, INTENT(in) :: numpar
            
 Description: Create a unit simplex of given dimensionality
    ! Calls: makesimplex, printsimplex
end subroutine SetupSimplex

public function HelixFit2gaussNewton (hits, kinks, lastplane, fitparMC, fitpar0, restype, maxiter, parstep, ambig0, doFitTref, res0, chi0, fitparB, resB, chiB, numiter) result (ierror)
    type (hits_type), INTENT(in) :: hits
            
 Fit data using the Gauss-Newton method for non-linear least squares

 hits data
    type (kinks_type), INTENT(in) :: kinks
             kink data
    integer, INTENT(in) :: lastplane
             last plane to use
    type (fitpar_type), INTENT(in) :: fitparMC
             "true" parameters
    type (fitpar_type), INTENT(in) :: fitpar0
             starting point parameters
    integer, INTENT(in) :: restype
             type of residuals
    integer, INTENT(in) :: maxiter
             maximum number of iterations
    real (kind=R8), INTENT(in), dimension (:) :: parstep
    integer, INTENT(in), dimension (:) :: ambig0
             ambiguity settings
    logical, INTENT(in) :: doFitTref
             enable external fitting of tref
    type (residuals_type), INTENT(out) :: res0
    type (chisqr_type), INTENT(out) :: chi0
    type (fitpar_type), INTENT(out) :: fitparB
    type (residuals_type), INTENT(out) :: resB
    type (chisqr_type), INTENT(out) :: chiB
    integer, INTENT(out) :: numiter
    integer :: ierror
             return value
    ! Calls: FitTref, clearRes, hf1, makeChiSqr, makeChisqr, printres, setAmbig, setResiduals
end function HelixFit2gaussNewton

private subroutine AngleTest0 (grad, maxcosa, iexcl)
    type (gradient_type), INTENT(in) :: grad
            
 Description: apply the Dennis-Broyden convergence test from
 p.273 in paper by J.E.Dennis "Non-Linear Least Squares and Equations"
 published in "The State of the Art in Numerical Analysis",
 edited by D.Jacobs, Academic Press, 1977.
 ISBN 0-12-378650-9
 QA297.C646 1976
    real (kind=R8), INTENT(out) :: maxcosa
    integer, INTENT(in) :: iexcl
end subroutine AngleTest0

private subroutine AngleTest (grad, maxcosa)
    type (gradient_type), INTENT(in) :: grad
    real (kind=R8), INTENT(out) :: maxcosa
    ! Calls: angletest0
end subroutine AngleTest

private function iterate_llsq (hits, kinks, lastplane, fitparMC, fitpar0, maxiter, doStopChebyshev, doFitTref, cutoff0, fitpar, res, grad, chi, numiter, costrms, remove2) result (ierror)
    type (hits_type), INTENT(in) :: hits
    type (kinks_type), INTENT(in) :: kinks
    integer, INTENT(in) :: lastplane
    type (fitpar_type), INTENT(in) :: fitparMC
    type (fitpar_type), INTENT(in) :: fitpar0
    integer, INTENT(in) :: maxiter
             maximum iterations
    logical, INTENT(in) :: doStopChebyshev
             stop when Chebyshev norm is zero
    logical, INTENT(in) :: doFitTref
             enable external fitting of tref
    real (kind=R8), INTENT(in) :: cutoff0
             drift L-R ambiguity cutoff
    type (fitpar_type), INTENT(inout) :: fitpar
    type (residuals_type), INTENT(inout) :: res
    type (gradient_type), INTENT(inout) :: grad
    type (chisqr_type), INTENT(inout) :: chi
    integer, INTENT(out) :: numiter
    real (kind=R8), INTENT(out) :: costrms
    logical, INTENT(in) :: remove2
    integer :: ierror
             return value
    ! Calls: AmbigCmpMC, AngleTest0, FitTref, copyRes, getFitparPPCP, hf1, makeChiSqr, printGradient, printres, setAmbig, setResiduals, updateFitpar
end function iterate_llsq

private function llsqstep (grad, res, dpar) result (ierror)
    type (gradient_type), INTENT(in) :: grad
            
 Calculate next step for the Gauss-Newton method for non-linear least squares
    type (residuals_type), INTENT(in) :: res
    real (kind=R8), INTENT(out), dimension (kMaxPar) :: dpar
    integer :: ierror
             return value
end function llsqstep

private function addGradient (grad, res, chi) result (ierror)
    type (gradient_type), INTENT(inout) :: grad
            
 Add new data point for gradient calculation
    type (residuals_type), INTENT(in) :: res
    type (chisqr_type), INTENT(in) :: chi
    integer :: ierror
             return value
    ! Calls: copyRes, die
end function addGradient

private function setupGradient (hits, kinks, lastplane, fitpar0, res0, parstep, grad) result (ierror)
    type (hits_type), INTENT(in) :: hits
            
 Setup initial data points for gradient calculations
    type (kinks_type), INTENT(in) :: kinks
    integer, INTENT(in) :: lastplane
    type (fitpar_type), INTENT(in) :: fitpar0
    type (residuals_type), INTENT(in) :: res0
    real (kind=R8), INTENT(in), dimension (kMaxPar) :: parstep
    type (gradient_type), INTENT(out) :: grad
    integer :: ierror
             return value
    ! Calls: SetupSimplex, makeChiSqr, setResiduals, updateFitpar
end function setupGradient

private function fitGradient (hits, fitpar, grad) result (ierror)
    type (hits_type), INTENT(in) :: hits
            
 For the gradient for the Gauss-Newton method for non-linear least squares

 The math is as follows:

 dT - grad*dP -> 0

 (dP^T * dP) * grad = dP^T * dT

 grad(ires) = (dP^T * dP)^-1 * (dP^T * dT)(ires)

 where: dT is the change in the track (mtrk)
        dP is the change in the parameters (par)

 compute the dP and dT matrices
    type (fitpar_type), INTENT(in) :: fitpar
    type (gradient_type), INTENT(inout) :: grad
    integer :: ierror
             return value
    ! Calls: setResidualsDrift
end function fitGradient

private subroutine printGradient (grad)
    type (gradient_type), INTENT(in) :: grad
end subroutine printGradient

private subroutine AmbigCmpMC (res, cost)
    type (residuals_type), INTENT(in) :: res
            
 Description: compare set ambiguities with geant ambiguities
    real (kind=R8), INTENT(in) :: cost
    ! Calls: hf1, hf2
end subroutine AmbigCmpMC

private function sqr (a) result (b)
    real (kind=R8), INTENT(in) :: a
             Description: return square of a real(r8) number
    real (kind=R8) :: b
end function sqr

private function myReal8 (a4) result (a8)
    real, INTENT(in) :: a4
    real (kind=R8) :: a8
end function myReal8