text:
    !---------------------------------------------------------------------------
    !
    !> \file    GEMLineSearch.f90
    !> \brief   Perform a line search for the GEMSolver.f90
    !> \author  M.H.A. Piro
    !> \date    Apr. 25, 2012
    !> \sa      GEMSolver.f90
    !> \sa      GEMNewton.f90
    !> \sa      CompChemicalPotential.f90
    !> \sa      CompFunctionNorm.f90
    !> \todo    Figure out a long term solution for the condition to allow the solution
    !!           to only perform 1 iteration if the functional norm is less than 1E-6.
    !! \todo    Figure out a more permanent solution to how dMaxGamma shoudl be handled.
    !!           I think that I need to get away from that and focus on ensuring that the
    !!           Wolfe conditions have been satisfied.
    !
    !
    ! Revisions:
    ! ==========
    !
    !   Date            Programmer      Description of change
    !   ----            ----------      ---------------------
    !   04/25/2012      M.H.A. Piro     Original code
    !   05/02/2012      M.H.A. Piro     Improved calculation of initial step
    !                                    length (constrain changes to the Element
    !                                    potentials).
    !   05/08/2012      M.H.A. Piro     Improved calculation of initial step length
    !                                    (constrain the maximum change to the total
    !                                    number of moles of a solutoin phase).
    !   04/11/2013      M.H.A. Piro     When computing a steplength that constrains
    !                                    the maximum change to the element potential
    !                                    to less than or equal to unity, exclude elements
    !                                    with zero moles.  In other words, exclude
    !                                    electrons corresponding to ionic phases
    !                                    that are not currently stable.  This is also
    !                                    done when updating the element potentials.
    !   05/08/2013      M.H.A. Piro     Changed tolerance for calculating the steplength
    !                                    when the number of moles of a species tends to zero
    !                                    to 1D-50 from 1D-100.
    !   06/08/2013      M.H.A. Piro     In determining an initial steplength, do not constrain
    !                                    the maximum decrease in the number of moles of a
    !                                    solution phase by a certain increment (e.g., 50%) if
    !                                    the number of moles of that phase is below a certain
    !                                    value (e.g., 10**(-9)).  The motivation for doing this
    !                                    is that a solution phase may be driving out of the system
    !                                    but it may be inhibited if this condition is not made.
    !   06/08/2013      M.H.A. Piro     Exit the Wolfe loop if the functional norm is below a
    !                                    certain tolerance (e.g., 10**(-6)).
    !   04/02/2014      M.H.A. Piro     I changed one of the conditions to satisfy the line search.
    !                                    Specifically, if the relative change of the functional
    !                                    norm is 0.95 < F_norm < 1.0 to 0.97 < F_norm < 1.0.
    !                                    Consider the following scenario: the functional norm
    !                                    is 1 at global iteration 5 and currently we are at iter 6.
    !                                    The f-norm of the first line search loop is 6.1 and then
    !                                    its 5.86 (dTemp = 0.961).  Before, the system would then
    !                                    exit, but the f-norm on the global scale is not converging.
    !   08/20/2015      M.H.A. Piro     Previously, the maximum change in the element potentials was
    !                                    constrained to 1.  Now, it is a linear function that varies
    !                                    with respect to the functional norm.  The motivation for this
    !                                    is that an initially poor guess may be far from equilibrium,
    !                                    and laxing this constraint accelerates convergence.  Another
    !                                    change is the way that dStepLength is initialized when
    !                                    dMolesSpecies tends to zero.  Previously, this would compute
    !                                    the minimum dStepLength that corresponds to reducing dMolesSpecies
    !                                    by 100.  Now, this is NOT applied when dMolesSpecies is incredibly
    !                                    small (e.g., less than numerical tolerance).
    !   12/18/2018      M.H.A. Piro     I effectively removed the constraint applied to the maximum
    !                                    change to the element potentials. I think that this is inefficient
    !                                    because it's probably best to leave it to ensuring that the
    !                                    Wolfe conditions have been satisfied. Furthermore, it seems to
    !                                    be an issue for SUBG phases, which are effectively on a knife edge.
    !
    !
    ! Purpose:
    ! ========
    !
    !> \details The purpose of this subroutine is to perform a line search using
    !! the direction vector computed by the Newton/Broyden solver.  The system
    !! is updated using an appropriate step length that satisfies the Wolfe
    !! conditions.  Specifically, values of dChemicalPotential and dMolesPhase
    !! are updated.  It is possible for the system of equations to be ill-
    !! behaved and yield inappropriate results.  An initial step-length is
    !! computed by normalizing the largest change of the system variables by a
    !! pre-defined value.  The maximum change to the element potentials is 1 and
    !! the maximum change to the number of moles of a solution phase is twice of
    !! the previous value.  For more information, refer to Chapter 6 of the
    !! above reference.
    !
    !
    ! Pertinent variables:
    ! ====================
    !
    ! nElements             The number of elements in the system.
    ! nConPhases            The number of pure condensed phases in the system.
    ! nSolnPhases           The number of solution phases in the system.
    ! dUpdateVar            A double real vector that contains updates to the
    !                        system variables (element potentials, moles of
    !                        solution phases and moles of pure condensed phases).
    ! dStepLength           Step length applied to the direction vector
    !                        (dUpdateVar)
    ! dMolesPhase           A double real vector representing the number of
    !                        moles of phases predicted to be stable.
    ! dLevel                The adjustment applied to the chemical potentials of
    !                        the elements
    ! iPhaseDampen          Integer vector that counts the number of times that
    !                        the number of moles
    !                        of a solution phase had to be dampened.
    ! dChemicalPotential    A double real vector representing the chemical
    !                        potential of each species and pure condensed phase.
    ! dGEMFunctionNorm      A double real scalar representing the norm of the
    !                        functional vector in the PGESolver.
    ! dGEMFunctionNormLast  A double real scalar representing the norm of the
    !                        functional vector from the previous iteration.
    !
    !---------------------------------------------------------------------------

subroutine GEMLineSearch

    USE ModuleThermo
    USE ModuleGEMSolver

    implicit none

    integer                       :: iterWolfe
    real(8)                       :: dStepLength, dTemp, dWolfeFunctionNormLast
    real(8), dimension(nElements) :: dElementPotentialLast
    real(8), dimension(nSpecies)  :: dMolesSpeciesLast
    logical                       :: lCompEverything

    ! Initialize variable:
    dWolfeFunctionNormLast  = 1D-12
    dMolesSpeciesLast       = dMolesSpecies
    dElementPotentialLast   = dElementPotential
    dGEMFunctionNormLast    = dGEMFunctionNorm
    dPartialExcessGibbsLast = dPartialExcessGibbs
    lCompEverything         = .FALSE.

    ! Initialize the line search method:
    call InitGEMLineSearch(dStepLength,dMolesSpeciesLast,dElementPotentialLast)

    ! Commence line search:
    LOOP_WOLFE: do iterWolfe = 1, 5

        ! Compute the fractional change in the functional norm:
        dTemp = dGEMFunctionNorm / dWolfeFunctionNormLast

        ! If the functional norm is already small, call it a day:
        if (dGEMFunctionNorm < 1D-6) exit LOOP_WOLFE

        ! If the number of moles of a solution phase is tending to zero, call it a day:
        if (MINVAL(dMolesPhase(nElements - nSolnPhases + 1: nElements)) < 1D-9) exit LOOP_WOLFE

        ! Check if the system is diverging or has sufficiently progressed, otherwise dampen:
        if (dGEMFunctionNorm < 0.999D0 * dGEMFunctionNormLast) then

            ! The system has sufficiently progressed, exit:
            exit LOOP_WOLFE

        elseif ((iterWolfe > 1).AND.(dTemp >= 0.999D0).AND.(dTemp <= 1D0)) then

            ! Return values to previous line search iteration:

            ! Update the system variables:
            dStepLength = 2D0

            call UpdateSystemVariables(dStepLength,dMolesSpeciesLast,dElementPotentialLast)

            ! Compute the chemical potentials of solution species:
            call CompChemicalPotential(lCompEverything)

            ! Compute the functional norm
            call CompFunctionNorm

            ! Values are returned to their previous values and exit:
            exit LOOP_WOLFE

        !elseif ((iterWolfe > 1).AND.(dTemp >= 0.95D0).AND.(dTemp <= 1D0)) then
        elseif ((iterWolfe > 1).AND.(dTemp >= 0.97D0).AND.(dTemp <= 1D0)) then

            ! The system has sufficiently progressed, exit:
            exit LOOP_WOLFE

        elseif ((iterWolfe > 1).AND.((dGEMFunctionNorm >= dWolfeFunctionNormLast))) then

            ! Return values to previous line search iteration:

            ! Update the system variables:
            dStepLength = 2D0

            call UpdateSystemVariables(dStepLength,dMolesSpeciesLast,dElementPotentialLast)

            ! Compute the chemical potentials of solution species:
            call CompChemicalPotential(lCompEverything)

            ! Compute the functional norm
            call CompFunctionNorm

            ! Values are returned to their previous values and exit:
            exit LOOP_WOLFE

        else
            ! Dampen the system variables:

            ! If the functional norm has only increased by a nominal amount (i.e., 1%), then exit:
            dTemp = dGEMFunctionNorm / dWolfeFunctionNormLast
            if ((dTemp > 1D0).AND.(dTemp < 1.05D0)) exit LOOP_WOLFE

            dStepLength = 0.5D0

            call UpdateSystemVariables(dStepLength,dMolesSpeciesLast,dElementPotentialLast)

            ! Compute the chemical potentials of solution species:
            call CompChemicalPotential(lCompEverything)

            dWolfeFunctionNormLast = dGEMFunctionNorm

            ! Compute the functional norm:
            call CompFunctionNorm

            ! Reiterate:
            cycle LOOP_WOLFE

        end if

    end do LOOP_WOLFE

    return

end subroutine GEMLineSearch

!-------------------------------------------------------------------------------------------------------------
!-------------------------------------------------------------------------------------------------------------
!-------------------------------------------------------------------------------------------------------------

    !---------------------------------------------------------------------------
    !
    ! Purpose:
    ! ========
    !
    ! The purpose of this subroutine is to initialize the line search algorithm.
    ! Specifically, the initial step length needs to be determined before the
    ! line search loop starts.
    !
    !
    ! Revisions:
    ! ==========
    !
    !   Date            Programmer      Description of change
    !   ----            ----------      ---------------------
    !
    !   05/08/2012      M.H.A. Piro     Original Code
    !   07/04/2012      M.H.A. Piro     If the number of moles of any soln
    !                                    phase is to be significantly reduced
    !                                    and the system is not stagnant, then
    !                                    then further dampen the system.
    !   09/29/2012      M.H.A. Piro     Revert the system if the phase assemblage
    !                                    has not changed in 50 iterations and
    !                                    the maximum change to the system
    !                                    variables is extremely large.
    !
    !
    ! Pertinent variables:
    ! ====================
    !
    ! dStepLength       A double real scalar representing the step length.
    ! dMaxIncrease      A double real scalar representing the maximum
    !                    increase in the number of moles of a solution phase.
    ! dMaxDecrease      A double real scalar representing the maximum
    !                    decrease in the number of moles of a solution phase.
    ! lCompEverything   A logical scalar indicating whether everything should
    !                    be computed in a particular subroutine (true) or
    !                    not (false).
    !
    !---------------------------------------------------------------------------

subroutine InitGEMLineSearch(dStepLength,dMolesSpeciesLast,dElementPotentialLast)

    USE ModuleThermo
    USE ModuleGEMSolver

    implicit none

    integer                       :: i, j, k, l, nMisciblePhases
    real(8)                       :: dStepLength, dTemp, dMaxIncrease, dMaxDecrease, dMaxChange, dMaxGamma
    real(8), dimension(nElements) :: dElementPotentialLast
    real(8), dimension(nSpecies)  :: dMolesSpeciesLast
    logical                       :: lCompEverything

    ! Initialize variables:
    lCompEverything  = .FALSE.
    dStepLength      = 1D0
    dMolesPhaseLast  = dMolesPhase
    nMisciblePhases  = 0

    ! Count the number of stable miscible phases:
    do j = 1, nSolnPhases
        k = -iAssemblage(nElements - j + 1)
        if (lMiscibility(k)) nMisciblePhases = nMisciblePhases + 1
    end do

    ! Initialize the maximum increase/decrease of functional variables
    ! (this is dapenned if the number of iterations gets large):
    if (iterGlobal < 1000) then
        dMaxIncrease = 2D0
        dMaxDecrease = 0.5D0
    elseif (iterGlobal < 2000) then
        dMaxIncrease = 1.5D0
        dMaxDecrease = 0.75D0
    else
        dMaxIncrease = 1.25D0
        dMaxDecrease = 0.85D0
    end if

    ! TEMPORARY: I should figure out a better and more permanent solution. My gut tells me that
    ! I need to go away from constraining the maximum change to the element potentials and a more
    ! elegant approach would be leaving it to the Wolfe conditions. For now, do this:
    if (iterGlobal > 1000) then
        dMaxGamma = (5D0 - 1D0) / (100D0 - 25D0) * (dGEMFunctionNorm - 100D0) + 5D0
        dMaxGamma = DMIN1(dMaxGamma,5D0)
        dMaxGamma = DMAX1(dMaxGamma,1D0)
    else
        dMaxGamma = 100D0
    end if

    ! Update the number of moles of pure condensed phases:
    do i = 1, nConPhases
        j = nElements + nSolnPhases + i
        dMolesPhase(i) = dUpdateVar(j)
    end do

    dMaxSpeciesChange = 0D0
    ! Initialize the step length (prevent moles from being negative):
    do l = 1, nSolnPhases
        k = -iAssemblage(nElements - l + 1)     ! Absolute solution phase index.

        do i = nSpeciesPhase(k-1) + 1, nSpeciesPhase(k)
            dTemp = 0D0
            do j = 1, nElements
                dTemp = dTemp + dUpdateVar(j) * dStoichSpecies(i,j)
            end do

            dTemp = dTemp / DFLOAT(iParticlesPerMole(i))

            ! NOTE: The variable dMolFraction is used temporarily to represent
            ! the fractional update to dMolesSpecies, but it does not replace
            ! dMolesSpecies in the event that further dampening is required.
            dMolFraction(i) = (1D0 + dUpdateVar(nElements + l) + dTemp - &
                dChemicalPotential(i))

            if (dMolFraction(i) /= 1D0) dTemp = 1D0 / (1D0 - dMolFraction(i))

            if ((dMolFraction(i) < 0D0).AND.(dTemp < dStepLength)) then
                if (dMolesSpecies(i) < 1D-50) then
                    ! This species is very small and the system is trying to make it negative.
                    ! Reduce the mass of this species by an arbitrary positive value less than
                    ! unity:
                    dMolFraction(i) = 0.5D0

!!!! TEMPORARY   !!!!!
!!!! TEMPORARY   !!!!!
!!!! TEMPORARY   !!!!!

                elseif ((dMolesSpecies(i) / dMolesPhaseLast(nElements - l + 1)) < 1D-10) then

                    dMolFraction(i) = 0.1D0

                else

                    dStepLength = dTemp * 0.99D0
                end if
            end if
            ! This check was getting hung up on tiny tiny mole fractions in some cases,
            ! so I've added the outer if to skip these
            if (dMolesSpecies(i) > 1D-100) then
                if (dMolFraction(i) > 0) then
                    dMaxSpeciesChange = MAX(dMaxSpeciesChange,ABS(LOG(dMolFraction(i)))*dStepLength)
                else
                    ! if set negative or zero do nothing, I guess
                end if
            end if
            dMolesSpecies(i) = dMolesSpecies(i) * dMolFraction(i)

        end do
    end do

    ! Initialize the steplength (constrain the element potentials to only
    ! change by dMaxGamma):
    do i = 1, nElements
        dTemp = DABS(dElementPotential(i) - dUpdateVar(i))

        if (dUpdateVar(i) == 0D0) cycle

        if (dTemp > 1D0) then

            dTemp = dMaxGamma / dTemp
            !dTemp = 1D0 / dTemp
            dStepLength = DMIN1(dTemp, dStepLength)
        end if
    end do

	! Update the element potentials:
    do i = 1, nElements
        dElementPotential(i) = dUpdateVar(i)
    end do

    ! Update the system variables:
    call UpdateSystemVariables(dStepLength,dMolesSpeciesLast,dElementPotentialLast)

    ! Constrain the number of moles of any solution phase to only change by a factor of 2:
    i = 0
    dStepLength = 1D0
    dTemp = 0D0
    do j = 1, nSolnPhases
        k     = nElements - j + 1     ! Solution phase index in iAssemblage and dMolesPhase

        if (dMolesPhase(k) /= dMolesPhaseLast(k)) then
            dTemp = dMolesPhaseLast(k) * (dMaxIncrease - 1D0) / (dMolesPhase(k) - dMolesPhaseLast(k))
        end if

        if (dTemp <= 0D0) dTemp = dStepLength
        dStepLength = DMIN1(dTemp, dStepLength)

        ! TEMPORARY TO AVOID AN INF:
        dMolesPhaseLast(k) = DMAX1(dMolesPhaseLast(k),1D-10)

        ! Do not allow the number of moles of a solution phase to be reduced by less than half if it is
        ! not being pushed to zero:
        dTemp = DABS(dMolesPhase(k) / dMolesPhaseLast(k))

        ! Count the number of times the number of moles of a solution phase is trying to increase
        ! by a large margin:
        if (dTemp > dMaxIncrease) i = i + 1

! THIS FOLLOWING SECTION NEEDS TO BE BETTER FIGURED OUT.  WHY WOULD IT DAMPEN IF A SOLUTOIN PHASE IS
! BECOMING NEGATIVE WITH DTEMP > 0.011 OR DTEMP < 0.009??

        if (dTemp < 1D0) then
        !if (dTemp < 0.1D0) then
            if ((dTemp > 0.05D0).OR.(dTemp < 0.009D0)) then
            !if ((dTemp > 0.011D0).OR.(dTemp < 0.009D0)) then

                ! TEMPORARY:
                if (dMolesPhase(k) < 1D-9) Cycle

                ! Count the number of solution phases that have molar quantities that are tending to zero:
                if (dTemp <= dMaxDecrease) i = i + 1

                ! Determine a step length that would result in halving this solution phase:
                dTemp = DABS(dMaxDecrease * dMolesPhaseLast(k) / (dMolesPhaseLast(k) - dMolesPhase(k)))
                dStepLength = DMIN1(dTemp,dStepLength)
            end if
        end if
    end do

    dTemp = 0.05D0

    ! Count the number of solution phases (i.e., j) that are changing by at least 5%:
    call CheckStagnation(dTemp,dMaxChange,j)
    ! Check if there is at least one solution phase that is tending to zero:
    if (i > 0) then
        ! An issue with the Gibbs Energy Minimization method is that one is effectively attempting to
        ! minimize two objective functions simultaneously: 1) the integral Gibbs energy of the system,
        ! and 2) the residual vector of the mass balance equations.  A common problem is that there
        ! may be a phase change that results in significant changes in both the mass balance residuals
        ! and the Gibbs energy function, making the functional norm a poor indicator of convergence.
        ! Generally, if a single solution phase should be forced out of the system, then the number of
        ! moles of that phase alone will change significantly.  If the number of moles of more than
        ! one solution phase changes, then the system may need to be further dampened.

        ! If there is more than one solution phase that is changing by more than 5%, then dampen some more:
        if (j > 1) dStepLength = dStepLength * 0.5D0

    else
        ! Stagnating? Seen a lot of oscillation in small dMolFractions, try reducing step size
        if (dMaxChange < 1D-6) dStepLength = dStepLength * 0.5D0

    end if

    ! Further dampen the system if a miscible phase was recently added to the system:
    if ((iterGlobal - iterLast < 10).AND.(nMisciblePhases > 0)) dStepLength = dStepLength * 0.5D0

    ! If the number of moles of a solution phase is changing by too large of an amount, then recompute
    if (dStepLength < 1D0) call UpdateSystemVariables(dStepLength,dMolesSpeciesLast,dElementPotentialLast)

    ! Compute the chemical potentials of solution species:
    call CompChemicalPotential(lCompEverything)

    ! Compute the functional norm:
    call CompFunctionNorm

end subroutine InitGEMLineSearch

!-------------------------------------------------------------------------------------------------------------
!-------------------------------------------------------------------------------------------------------------
!-------------------------------------------------------------------------------------------------------------

    !---------------------------------------------------------------------------
    !
    ! Purpose:
    ! ========
    !
    ! The purpose of this subroutine is to update the system variables using
    ! a specified steplength.
    !
    !---------------------------------------------------------------------------

subroutine UpdateSystemVariables(dStepLength,dMolesSpeciesLast,dElementPotentialLast)

    USE ModuleThermo
    USE ModuleGEMSolver

    integer                       :: i, j, k
    real(8)                       :: dTemp, dStepLength
    real(8), dimension(nSpecies)  :: dMolesSpeciesLast
    real(8), dimension(nElements) :: dElementPotentialLast

    ! Loop through all solution phases expected to be stable:
    do j = 1, nSolnPhases
        k     = -iAssemblage(nElements - j + 1)
        dTemp = 0D0
        do i = nSpeciesPhase(k-1) + 1, nSpeciesPhase(k)
            dMolesSpecies(i) = dStepLength * dMolesSpecies(i) + (1D0 - dStepLength) * dMolesSpeciesLast(i)
            dMolesSpecies(i) = DMAX1(dMolesSpecies(i), dTolerance(8))
            dTemp            = dTemp + dMolesSpecies(i)
        end do
        ! Compute the total number of moles of each solution phase:
        dMolesPhase(nElements-j+1) = dTemp
    end do

    ! Dampen the number of moles of pure condensed phases:
    do i = 1, nConPhases
        dMolesPhase(i) = dStepLength * dMolesPhase(i) + (1D0 - dStepLength) * dMolesPhaseLast(i)
    end do

    ! Dampen the element potentials:
    LOOP_Gamma: do i = 1, nElements
        if (dElementPotential(i) == 0D0) then
            dElementPotential(i) = dElementPotentialLast(i)
            cycle LOOP_Gamma
        end if
        dElementPotential(i) = dStepLength * dElementPotential(i) + (1D0 - dStepLength) * dElementPotentialLast(i)
    end do LOOP_Gamma

end subroutine UpdateSystemVariables

    !---------------------------------------------------------------------------
    !                            END - GEMLineSearch.f90
    !---------------------------------------------------------------------------

max_stars_repo_name: MULTIFROZ/thermochimica

max_issues_repo_name: nikolasscuro/thermochimica

max_forks_repo_name: nikolasscuro/thermochimica


text:
!	Find posterior distribution for latent vectors.

!	catobsrange is used to relate observed and underlying categories.
!	Responses are in dat.
!	maxita is used to determine the number of iterations to use to find the maximum posterior density.
!	numcat provides ranges for underlying observations,
!	numcatobs provides ranges of observations.
!	mask indicates items to use for computations.
!	normal determines if a normal or multinomial analysis is used.
!	beta is the parameter vector.
!	changemin is used for minimum change.
!	Predictors are in indvar.
!	lintran is the linear transformation for the latent vector.
!	maxdalpha is used to limit changes in alpha approximations.
!	quadpoint provides quadrature points.
!	quadweight provides quadrature weights.
!	tau controls shrinking of step size.
!	tola is used to determine the desired accuracy of the maximization.
!	tolsing is used to determine the tolerance for the Cholesky decomposition.
!	alpha contains locations of maxima.
!	cholnhess contains modified Cholesky decompositions.
!	postdensity gives weights at quadrature points.
!	theta are quadrature points.
subroutine posterior(catobsrange,dat,maxita,numcat,numcatobs,mask,normal,beta,changemin,indvar,lintran,&
	maxdalpha,quadpoint,quadweight,&
	tau,tola,tolsing,alpha,cholnhess,postdensity,theta)
implicit none
interface

	function chol(sym,tolsing)
!	chol is used to compute a modified Cholesky decomposition of the n by n matrix sym.
!	The singularity tolerance is tolsing.
		implicit none
		real(kind=8),intent(in)::sym(:,:),tolsing
	 	real(kind=8)::chol(size(sym,1),size(sym,1))
	end function chol

!	Get normal density for scaled latent vector.  Place in density.
!	Input involves linear component linth, quadrature points quadpoint, quadratic components quadth, quadrature weights quadweight, and scale
!	factor scaletheta.
	function density(linth,quadpoint,quadth,quadweight,scaletheta)
		implicit none
		real(kind=8),intent(in)::linth(:),quadpoint(:,:),quadth(:,:),quadweight(:),scaletheta
		real(kind=8)::density(size(quadweight))
	end function density
!	Get multinomial probabilities for latent vector.  Place in density.
!	Input involves linear component linth, quadrature points quadpoint, quadratic components quadth, and quadrature weights quadweight.
	function densitym(linth,quadpoint,quadth,quadweight)
		real(kind=8),intent(in)::linth(:),quadpoint(:,:),quadth(:,:),quadweight(:)
		real(kind=8)::densitym(size(quadweight))
	end function densitym

!	convert an n by n symmetric matrix from compact to regular form.
	function expandmat(n,compmat)
		implicit none
		integer,intent(in)::n
		real(kind=8),intent(in)::compmat(:)
		real(kind=8)::expandmat(n,n)
	end function expandmat

!	Extract the linear component lintheta and the quadratic component quadtheta from the parameter vector beta.
	subroutine getlinquad(beta,lintheta,quadtheta)
		implicit none
		real(kind=8),intent(in)::beta(:)
		real(kind=8),intent(out)::lintheta(:,:),quadtheta(:,:,:)
	end subroutine getlinquad

!	Extract the linear component linth and the quadratic component quadth from lintheta and quadtheta for a given predictor thpr.
	subroutine getlinquadth(lintheta,quadtheta,thpr,linth,quadth)
		implicit none
		real(kind=8),intent(in)::lintheta(:,:),quadtheta(:,:,:),thpr(:)
		real(kind=8),intent(out)::linth(:),quadth(:,:)
	end subroutine getlinquadth
!
!	Find location of maximum posterior density of the latent vector.
!	catobsrange is used to relate observed and underlying categories.
!	maxita is used to determine the number of iterations to use to find the maximum posterior density.
!	maxitb is used for subiterations.
!	npred is the number of predictors.
!	numcat provides ranges for underlying observations,
!	numcatobs provides ranges of observations, mask indicates items to use for computations.
!	obsmask is an observation mask that can be used to restrict attention to a portion of the response.
!	beta is the parameter vector.
!	changemin is used for minimum change.
!	linth is the linear component of the log density.
!	lintran is the linear transformation of the latent vector.
!	maxdtheta is the maximum permitted change in the approximation of the maximum in one step.
!	quadth is the quadratic component of the log density.
!	resp is the response.
!	tau is used for minimum step sizes.
!	tol is used to determine the desired accuracy of the maximization.
!	tolsing is used to determine the tolerance for the Cholesky decomposition.
!	nhchol is the modified Cholesky decomposition of the negative hessian.
!	theta is the maximum.
	subroutine maxpost(catobsrange,maxita,maxitb,npred,numcat,numcatobs,resp,obsmask,beta,changemin,&
		linth,lintran,maxdtheta,quadth,tau,tol,tolsing,&
		nhchol,theta)
		implicit none
		integer,intent(in)::catobsrange(:,:),maxita,maxitb,npred,numcat(:),numcatobs(:),resp(:)
		logical,intent(in)::obsmask(:)
		real(kind=8),intent(in)::beta(:),changemin,lintran(:,:),linth(:),maxdtheta,quadth(:,:),tau,tol,tolsing
		real(kind=8),intent(inout)::theta(:)
		real(kind=8),intent(out)::nhchol(:,:)
	end subroutine maxpost

!	probvec is used to compute underlying item probabilities.
!   	numcat indicates the number of underlying categories per item.
!	mask indicates which items were presented.
!	locations provides location parameters and scales provides scale factors.
!	theta is the value of the latent vector.
	function probvec(numcat,mask,locations,scales,theta)
		implicit none
		integer,intent(in)::numcat(:)
		logical,intent(in)::mask(:)
		real(kind=8),intent(in)::locations(:),scales(:,:),theta(:)
		real(kind=8)::probvec(size(locations))
	end function probvec

!	probvecobs is used to compute observed item probabilities given the latent vector.
!	catobsrange is the table of ranges of underlying categories per observed category.
!	numcatobs provides the number of observed categories per item.
!	resp is the response vector.
!	mask indicates which items were presented.
!	probcat is the vector of underlying item probabilities.
	function probvecobs(catobsrange,numcatobs,resp,mask,probcat)
		implicit none
		integer,intent(in)::catobsrange(:,:),numcatobs(:),resp(:)
		logical,intent(in)::mask(:)
		real(kind=8),intent(in)::probcat(:)
		real(kind=8)::probvecobs(size(mask))
	end function probvecobs

!	Convert from adaptive quadrature scale to normalized latent vector.
!	Use nhchol for modified Cholesky decomposition of negative hessian.
!	Use oldalpha for location.
!	Use scalez for scaling.
!	Quadrature point is z.
	function scaleadapt(nhchol,oldalpha,scalez,z)
		implicit none
		real(kind=8),intent(in)::nhchol(:,:),oldalpha(:),scalez(:),z(:)
		real(kind=8)::scaleadapt(size(z))
	end function scaleadapt

!		This program uses the modified Cholesky decomposition
!		in lu to solve the equation lusolve = b.
	function solve(lu,b)
		implicit none

		real(kind=8),intent(in)::lu(:,:)
		real(kind=8),intent(in)::b(:)
		real(kind=8)::solve(size(b))
	end function solve
end interface
integer,intent(in)::catobsrange(:,:),dat(:,:),maxita,numcat(:),numcatobs(:)
logical,intent(in)::mask(:),normal
real(kind=8),intent(in)::beta(:),changemin,indvar(:,:),lintran(:,:),&
	maxdalpha,quadpoint(:,:),quadweight(:),&
	tau,tola,tolsing
real(kind=8),intent(inout)::alpha(:,:),cholnhess(:,:,:)
real(kind=8),intent(out)::postdensity(:,:),theta(:,:,:)

!	dimlatin is the dimension of the latent vector.
!	dimlatout is the dimension of the transformed latent vector.

!	ncat is the total number of underlying categories.
!	nitems is the number of items, and item is an item.
!	nobs is number of observations.
!	npred is the number of predictors.
!	nquad is the number of quadrature points.
!	nscale is the pointer for scale parameters.
!	nscale1 points to the end of the scale parameters.
!	obs is observation number.
!	pred is a predictor counter.
!	quad counts quadrature points.
!	resp is a single response.
!	row is a row counter.
integer::dimlatin,dimlatout,item,ncat,nitems,&
	nobs,npred,nquad,nscale,&
	nscale1,obs,pred,quad,&
	resp(size(mask)),row
!	Data mask.
logical::datamask(size(mask))
!	alpha is used to find the maximum posterior density.
!	locations is for location parameters.
!	linth is a linear component of the log density of a latent vector.
!	lintheta is for linear parameters for the log density.
!	meanth is the mean of the latent vector.
!	newtheta is lintran(theta).
!	nhchol is a Cholesky decomposition of -2*quadth.
!	priordensity provides a prior density.
!	probcat is used for underlying marginal conditional probabilities.
!	probcatobs is used for observed marginal conditional probabilities.
!	probobs is the probability of an observation.
!	quadth is a quadratic component of the log density of a latent vector.
!	quadtheta is for quadratic parameters for the log density.
!	scales is for scale parameters.
!	scaletheta is used for scaling.
!	scalez is used for scaling.
!	thpr is the predictor vector for an observation.
!	z is a transformation of a latent vector value.
real(kind=8)::linth(size(quadpoint,1)),&
	lintheta(size(quadpoint,1),size(indvar,1)),locations(sum(numcat)),&
	newtheta(size(lintran,1)),meanth(size(quadpoint,1)),&
	nhchol(size(quadpoint,1),size(quadpoint,1)),&
	priordensity(size(quadweight)),probcat(sum(numcat)),probcatobs(size(mask)),probobs,&
	quadth(size(quadpoint,1),size(quadpoint,1)),&
	quadtheta(size(quadpoint,1),size(quadpoint,1),size(indvar,1)),&
	scales(size(lintran,1),sum(numcat)),scaletheta,&
	scalez(size(quadpoint,1)),&
	thpr(size(indvar,1)),z(size(quadpoint,1))

!	Set up parameter arrays for simplified processing.
dimlatin=size(quadpoint,1)

dimlatout=size(lintran,1)
ncat=sum(numcat)
nitems=size(mask)
npred=size(indvar,1)

nobs=size(dat,2)
!	Item scale parameters.
nscale=ncat+1

nscale1=ncat*(1+dimlatout)
nquad=size(quadweight)
!	Parameters for quadratic component of log density.

locations=beta(1:ncat)
scales=reshape(beta(nscale:nscale1),(/dimlatout,ncat/))

call getlinquad(beta,lintheta,quadtheta)

!	Obtain results for each observation.
do obs=1,nobs

! Observation obs.
	resp=dat(:,obs)

!	Density specification for observation.
	thpr=indvar(:,obs)

	call getlinquadth(lintheta,quadtheta,thpr,linth,quadth)

!	Mask
	datamask=.false.
	do item=1,nitems
		if(resp(item)>=0.and.resp(item)<numcatobs(item).and.mask(item))datamask(item)=.true.
	end do

!	The normal versus multinomial cases differ here.
	if(normal)then
!	Mean and covariance matrix of latent vector.

		nhchol=chol(-2.0_8*quadth,tolsing)
!	Check that covariance matrix is positive definite.
		if(minval((/(nhchol(row,row),row=1,dimlatin)/))<=0.0_8)stop "Unacceptable covariance matrix"
		meanth=solve(nhchol,linth)
		scaletheta=sqrt(product((/(nhchol(row,row),row=1,dimlatin)/)))
		scaletheta=scaletheta*exp(-0.5_8*dot_product(linth,meanth))

!	See if adaptive approach to be used.
		if(maxita>0) then

			call maxpost(catobsrange,maxita,maxita,npred,numcat,numcatobs,resp,datamask,beta,changemin,&
				linth,lintran,maxdalpha,quadth,tau,tola,tolsing,&
				nhchol,alpha(:,obs))
			cholnhess(:,:,obs)=nhchol
		end if
		nhchol=cholnhess(:,:,obs)
		do row=1,dimlatin
			scalez(row)=1.0_8/sqrt(nhchol(row,row))
		end do
		scaletheta=scaletheta*product(scalez)
		do quad=1,nquad
				theta(:,quad,obs)=scaleadapt(nhchol,alpha(:,obs),scalez,quadpoint(:,quad))
		end do
		priordensity=density(linth,theta(:,:,obs),quadth,quadweight,scaletheta)
	else
		theta(:,:,obs)=quadpoint
		priordensity=densitym(linth,quadpoint,quadth,quadweight)
	end if
!	Cycle through quadrature points.

	do quad=1,nquad

! Linear transformation of latent vector.
		newtheta=matmul(lintran,theta(:,quad,obs))

! Probability computation.

		probcat=probvec(numcat,datamask,locations,scales,newtheta)
		probcatobs=probvecobs(catobsrange,numcatobs,resp,datamask,probcat)
		postdensity(quad,obs)=priordensity(quad)
		if(any(datamask))postdensity(quad,obs)=postdensity(quad,obs)*product(probcatobs,datamask)
	end do
	probobs=sum(postdensity(:,obs))
	postdensity(:,obs)=postdensity(:,obs)/probobs
end do
return
end subroutine posterior

max_stars_repo_name: EducationalTestingService/MIRT

max_issues_repo_name: EducationalTestingService/MIRT

max_forks_repo_name: EducationalTestingService/MIRT


text:
c     prog6
c     c.l.lawson and r.j.hanson, jet propulsion laboratory, 1973 jun 15
c     to appear in 'solving least squares problems', prentice-hall, 1974
c	   demonstrate the use of the subroutine   ldp	for least
c     distance programming by solving the constrained line data fitting
c     problem of chapter 23.
c
      dimension e(4,2), f(4), g(3,2), h(3), g2(3,2), h2(3), x(2), z(2),
     1w(4)
      dimension wldp(21), s(6), t(4)
      integer index(3)
c
      write (6,110)
      mde=4
      mdgh=3
c
      me=4
      mg=3
      n=2
c		 define the least squares and constraint matrices.
      t(1)=0.25
      t(2)=0.5
      t(3)=0.5
      t(4)=0.8
c
      w(1)=0.5
      w(2)=0.6
      w(3)=0.7
      w(4)=1.2
c
	  do 10 i=1,me
	  e(i,1)=t(i)
	  e(i,2)=1.
   10	  f(i)=w(i)
c
      g(1,1)=1.
      g(1,2)=0.
      g(2,1)=0.
      g(2,2)=1.
      g(3,1)=-1.
      g(3,2)=-1.
c
      h(1)=0.
      h(2)=0.
      h(3)=-1.
c
c     compute the singular value decomposition of the matrix, e.
c
      call svdrs (e,mde,me,n,f,1,1,s)
c
      write (6,120) ((e(i,j),j=1,n),i=1,n)
      write (6,130) f,(s(j),j=1,n)
c
c     generally rank determination and levenberg-marquardt
c     stabilization could be inserted here.
c
c	 define the constraint matrix for the z coordinate system.
	  do 30 i=1,mg
	      do 30 j=1,n
	      sm=0.
		  do 20 l=1,n
   20		  sm=sm+g(i,l)*e(l,j)
   30	      g2(i,j)=sm/s(j)
c	  define constraint rt side for the z coordinate system.
	  do 50 i=1,mg
	  sm=0.
	      do 40 j=1,n
   40	      sm=sm+g2(i,j)*f(j)
   50	  h2(i)=h(i)-sm
c
      write (6,140) ((g2(i,j),j=1,n),i=1,mg)
      write (6,150) h2
c
c			 solve the constrained problem in z-coordinates.
c
      call ldp (g2,mdgh,mg,n,h2,z,znorm,wldp,index,mode)
c
      write (6,200) mode,znorm
      write (6,160) z
c
c		     transform back from z-coordinates to x-coordinates.
	  do 60 j=1,n
   60	  z(j)=(z(j)+f(j))/s(j)
	  do 80 i=1,n
	  sm=0.
	      do 70 j=1,n
   70	      sm=sm+e(i,j)*z(j)
   80	  x(i)=sm
      res=znorm**2
      np1=n+1
	  do 90 i=np1,me
   90	  res=res+f(i)**2
      res=sqrt(res)
c				     compute the residuals.
	  do 100 i=1,me
  100	  f(i)=w(i)-x(1)*t(i)-x(2)
      write (6,170) (x(j),j=1,n)
      write (6,180) (i,f(i),i=1,me)
      write (6,190) res
      stop
c
  110 format (46h0prog6..  example of constrained curve fitting,26h usin
     1g the subroutine ldp.,/43h0related intermediate quantities are giv
     2en.)
  120 format (10h0v =	   ,2f10.5/(10x,2f10.5))
  130 format (10h0f tilda =,4f10.5/10h0s =	,2f10.5)
  140 format (10h0g tilda =,2f10.5/(10x,2f10.5))
  150 format (10h0h tilda =,3f10.5)
  160 format (10h0z =	   ,2f10.5)
  170 format (52h0the coeficients of the fitted line f(t)=x(1)*t+x(2),12
     1h are x(1) = ,f10.5,14h	and x(2) = ,f10.5)
  180 format (30h0the consecutive residuals are/1x,4(i10,f10.5))
  190 format (23h0the residuals norm is ,f10.5)
  200 format (18h0mode (from ldp) =,i3,2x,7hznorm =,f10.5)
c
      end

ext: f

max_stars_repo_name: ysBach/irafdocgen

max_issues_repo_name: kirxkirx/iraf

max_forks_repo_name: kirxkirx/iraf


text:
      MODULE lodcslmpi_I
      INTERFACE
!...Translated by Pacific-Sierra Research 77to90  4.3E  14:04:58   1/ 3/07
!...Modified by Charlotte Froese Fischer
!                     Gediminas Gaigalas  10/05/17
      SUBROUTINE lodcslmpi (nfile, ncore, jblock)
      INTEGER :: nfile, ncore, jblock
      END SUBROUTINE
      END INTERFACE
      END MODULE

max_stars_repo_name: sylas/grasp-continuum

max_issues_repo_name: sylas/grasp-continuum

max_forks_repo_name: sylas/grasp-continuum


text:
MODULE Circle
   REAL, PARAMETER :: Pi = 3.1415927
   REAL :: radius
END MODULE Circle

max_stars_repo_path: test cases/fortran/2 modules/stuff.f90

max_issues_repo_path: test cases/fortran/2 modules/stuff.f90

max_forks_repo_path: test cases/fortran/2 modules/stuff.f90


text:
C
C
C --- FREDA --- ORIGINALLY WRITTEN BY SAM MOTHERWELL
C MODIFIED OVER THE YEARS BY ELEANOR DODSON YORK
C
      SUBROUTINE  FREDA(IMIN,IMAX,FL,NCHAR)
      CHARACTER*1 LE,FL(40,NCHAR),COL,INUM,IBLANK,IPLUS,IMINUS,IPOINT
      INTEGER ICHAR(40),IBEGIN(40)
      COMMON /FR1/  COL(80)
      COMMON /FR2/  FN(40),KN,KL
      DIMENSION NUM(10),INUM(10)
      DATA  NUM/1,2,3,4,5,6,7,8,9,0/
      DATA INUM/'1','2','3','4','5','6','7','8','9','0'/
      DATA IBLANK,IPLUS,IMINUS,IPOINT,LE/' ','+','-','.','e'/
C
      DO 1 KN=1,40
      FN(KN)=0.
      DO 1 KL=1,NCHAR
      FL(KN,KL)=IBLANK
1     CONTINUE
      KL=0
      KN=0
C
CXC  GET WORD LIMITS  EACH WORD SEPERATED BY SPACE.
      I1=IMIN-1
      J1=0
      DO 2 I=1,40
      ICHAR(I)=0
      IBEGIN(I)=0
2     CONTINUE
C
10    J1=J1+1
11    I1=I1+1
      IF(I1.GT.IMAX) J1=J1-1
      IF(I1.GT.IMAX) GO TO 20
      IF(COL(I1).EQ.IBLANK) GO TO 11
      ICHAR(J1)=1
      IBEGIN(J1)=I1
12    I1=I1+1
      IF(I1.GT.IMAX) GO TO 20
      IF(COL(I1).EQ.IBLANK) GO TO 10
      ICHAR(J1)=ICHAR(J1)+1
      GO TO 12
CC
CC
CC   UNSCRAMBLE J1 WORDS
20    KN=0
      KL=0
      DO 21 J=1,J1
      IF(ICHAR(J).EQ.0) GO TO 21
      IFIRST=IBEGIN(J)
      ILAST=IFIRST+ICHAR(J)-1
      IDIG=0
      NUMBE=0
      ILIT=0
      DO 22 I=IFIRST,ILAST
      JDIG=IDIG
      DO 23 K=1,10
      IF(COL(I).EQ.INUM(K)) IDIG=IDIG+1
23    CONTINUE
      IF(COL(I).EQ.IPOINT) IDIG=IDIG+1
      IF(COL(I).EQ.IPLUS) IDIG=IDIG+1
      IF(COL(I).EQ.IMINUS) IDIG=IDIG+1
      IF(IDIG.GT.JDIG) GO TO 22
      IF(COL(I).EQ.LE) NUMBE=I
      ILIT=1+ILIT
22    CONTINUE
      IF(NUMBE.GT.0.AND.ILIT.EQ.1.AND.IDIG.GT.0) GO TO 40
      IF(ILIT.GT.0)GO TO 30
C
CC   NUMERIC FILED
      KN=KN+1
      A=0.0
      IP=0
      NP=0
      S=1.
      DO 24 I=IFIRST,ILAST
      DO 25 K=1,10
      IF(COL(I).EQ.INUM(K)) GO TO 26
25    CONTINUE
      IF(COL(I).EQ.IPOINT) IP=1
      IF(COL(I).EQ.IPLUS) S=1.
      IF(COL(I).EQ.IMINUS) S=-1.
      GO TO 24
26    IF(K.EQ.10) K=0
      IF(IP.EQ.1) NP=NP+1
      A=10.*A + K
24    CONTINUE
      FNK=S*A/10.**NP
      IF(IP.EQ.0)THEN
        NFK=NINT(FNK)
        FN(KN)=NFK
      ELSE
        FN(KN)=FNK
      ENDIF
      GO TO 21
CC   NUMERIC E TYPE FIELD
40    KN=KN+1
      ILASE=NUMBE-1
      A=0.0
      IP=0
      NP=0
      S=1.
      DO 44 I=IFIRST,ILASE
      DO 45 K=1,10
      IF(COL(I).EQ.INUM(K)) GO TO 46
45    CONTINUE
      IF(COL(I).EQ.IPOINT) IP=1
      IF(COL(I).EQ.IPLUS) S=1.
      IF(COL(I).EQ.IMINUS) S=-1.
      GO TO 44
46    IF(K.EQ.10) K=0
      IF(IP.EQ.1) NP=NP+1
      A=10.*A + K
44    CONTINUE
      A1=S*A/10.**NP
      A=0.0
      IFIRSE=NUMBE+1
      IP=0
      NP=0
      S=1.
      DO 54 I=IFIRSE,ILAST
      DO 55 K=1,10
      IF(COL(I).EQ.INUM(K)) GO TO 56
55    CONTINUE
      IF(COL(I).EQ.IPLUS) S=1.
      IF(COL(I).EQ.IMINUS) S=-1.
      GO TO 54
56    IF(K.EQ.10) K=0
      A=10.*A + K
54    CONTINUE
      A=S*A
      FN(KN)=A1*10.**A
      GO TO 21
C
C
C  LITERAL FIELD(S)
30    KL=KL+1
      I2=1
      DO 32 I=IFIRST,ILAST
      FL(KL,I2)=COL(I)
      I2=I2+1
      IF(I.EQ.ILAST) GO TO 32
      IF(I2.GT.NCHAR) KL=KL+1
      IF(I2.GT.NCHAR)I2=1
32    CONTINUE
21    CONTINUE
      RETURN
      END


text:
! DART software - Copyright UCAR. This open source software is provided
! by UCAR, "as is", without charge, subject to all terms of use at
! http://www.image.ucar.edu/DAReS/DART/DART_download
!
! $Id$

program test_inv_gamma

! test the gamma distribution random number generator routine
! with different means, standard deviations, and iteration counts.

use      types_mod, only : r8
use  utilities_mod, only : register_module, &
                           open_file, close_file, &
                           initialize_utilities, finalize_utilities, &
                           squeeze_out_blanks
use random_seq_mod, only : random_seq_type, init_random_seq, random_inverse_gamma

implicit none

! version controlled file description for error handling, do not edit
character(len=256), parameter :: source   = &
   "$URL$"
character(len=32 ), parameter :: revision = "$Revision$"
character(len=128), parameter :: revdate  = "$Date$"

type (random_seq_type) :: r
integer :: i, j, n, f
real(r8) :: r1, sum, sumsq, a, b, l
real(r8) :: mean, var, sd, compmean, compvar, compsd, compdiffm, compdiffsd
logical  :: write_this_one
character(len=256) :: fname, temp

character(len=50) :: formf = '(I12,2(F8.3),4(F12.5),2(F16.10),2(F12.5))'
logical :: write_me = .true.       ! if true, write each distribution into a file for later diagnostics
integer :: write_limit = 1000000   ! but only if rep count is not greater than this limit

type t_inputs
 real(r8) :: t_alpha
 real(r8) :: t_beta
 integer  :: t_nreps
end type

! to add more tests or change the parameters, specify a test count
! and update the sets of inputs here:

integer, parameter :: ntests = 27
type(t_inputs) :: t(ntests) = (/ &
     t_inputs(  3.0_r8,  0.20_r8,       100), &
     t_inputs(  3.0_r8,  0.20_r8,     10000), &
     t_inputs(  3.0_r8,  0.20_r8,   1000000), &
     t_inputs(  3.0_r8,  0.20_r8, 100000000), &
     t_inputs(  5.0_r8,  0.40_r8,   4000000), &
     t_inputs(  5.0_r8,  2.00_r8,   5000000), &
     t_inputs(  5.0_r8,  0.60_r8,       100), &
     t_inputs(  5.0_r8,  0.60_r8,     10000), &
     t_inputs(  5.0_r8,  0.60_r8,   1000000), &
     t_inputs(  5.0_r8,  0.60_r8, 100000000), &
     t_inputs(  3.0_r8,  0.50_r8,   1000000), &
     t_inputs(  4.0_r8,  0.50_r8,   1000000), &
     t_inputs(  5.0_r8,  0.50_r8,   1000000), &
     t_inputs(  6.0_r8,  0.50_r8,   1000000), &
     t_inputs(  6.0_r8,  0.01_r8,   1000000), &
     t_inputs( 86.0_r8,  0.01_r8,   1000000), &
     t_inputs(  6.6_r8,  0.01_r8,   1000000), &
     t_inputs(  6.0_r8,  0.01_r8,   1000000), &
     t_inputs(956.0_r8,  0.01_r8,   1000000), &
     t_inputs( 36.0_r8,  0.01_r8,   1000000), &
     t_inputs(  6.0_r8,  0.03_r8,   1000000), &
     t_inputs(  6.0_r8,  0.09_r8,   1000000), &
     t_inputs(  6.0_r8,  9.00_r8,   1000000), &
     t_inputs(  6.0_r8,  0.10_r8,   1000000), &
     t_inputs(  6.0_r8,  0.80_r8,   1000000), &
     t_inputs(  6.0_r8,  0.50_r8,   1000000), &
     t_inputs(  6.0_r8,  3.50_r8,   1000000)  /)

call initialize_utilities('test_inv_gamma')
call register_module(source,revision,revdate)

write(*, *) ''
write(*, *) 'sample size     input a & b       computed mean, sd         actual mean, sd       diff mean         diff sd        % diff mean,  sd'
write(*, *) ''

do j=1, ntests

   call init_random_seq(r, 5)

   a = t(j)%t_alpha
   b = t(j)%t_beta
   n = t(j)%t_nreps

   ! save all values in a file for post-plotting?
   write_this_one = (write_me .and. n <= write_limit)

   if (write_this_one) then
      write(temp, "(A,F8.3,A,F8.3,A,I10)") "invgamma_", a, "_", b, "_", n
      call squeeze_out_blanks(temp, fname)
      f = open_file(fname)
   endif

   ! analytical values:

   ! lambda = scale, beta = rate
   l = 1.0_r8 / b

   if (a > 1.0_r8) then
      mean  = b / (a - 1.0_r8)
   else
      mean  = b  ! ??
   endif
   if (a > 2.0_r8) then
      sd  = sqrt((b*b) / ((a - 1.0)**2 * (a - 2.0_r8)))
   else
      sd  = b    ! ??
   endif

   sum = 0.0_r8
   sumsq = 0.0_r8

   ! inverse gamma wants shape and scale.
   do i = 1, n
      r1 = random_inverse_gamma(r, a, l)

      if (write_this_one) write(f,*) r1

      sum = sum + r1
      sumsq = sumsq + (r1 - mean)**2
   end do

   if (write_this_one) call close_file(f)

   ! computed values:
   compmean = sum / n
   compsd = sqrt(sumsq/n)

   ! differences
   compdiffm  = compmean - mean
   compdiffsd = compsd - sd

   write(*, formf) n, a, b, mean, sd, compmean, compsd, compdiffm, compdiffsd, &
                      abs(compdiffm/mean) * 100._r8, abs(compdiffsd/sd) * 100._r8

enddo

call finalize_utilities()

end program test_inv_gamma

! <next few lines under version control, do not edit>
! $URL$
! $Id$
! $Revision$
! $Date$

max_stars_repo_name: hkershaw-brown/feature-preprocess

max_issues_repo_name: hkershaw-brown/feature-preprocess

max_forks_repo_name: hkershaw-brown/feature-preprocess