text:
module type_backward_euler
#include <messenger.h>
    use messenger,                      only: write_line
    use mod_kinds,                      only: rk,ik
    use mod_constants,                  only: ONE, NO_ID, dQ_DIFF, NO_DIFF
    use mod_spatial,                    only: update_space

    use type_time_integrator_marching,  only: time_integrator_marching_t
    use type_system_assembler,          only: system_assembler_t

    use type_chidg_data,                only: chidg_data_t
    use type_nonlinear_solver,          only: nonlinear_solver_t
    use type_linear_solver,             only: linear_solver_t
    use type_preconditioner,            only: preconditioner_t
    use type_solver_controller,         only: solver_controller_t
    use type_chidg_vector
    use type_chidg_matrix,              only: chidg_matrix_t
    use type_element_info,              only: element_info_t

    implicit none
    private

    !>  Object implementing the backward Euler time integrator
    !!
    !!  @author Mayank Sharma
    !!  @date   3/31/2017
    !!
    !-------------------------------------------------------------------------------
    type, extends(time_integrator_marching_t),  public  :: backward_euler_t

    contains
        procedure   :: init
        procedure   :: step
    end type backward_euler_t
    !*******************************************************************************

    !>  Object for assembling the implicit system
    !!
    !!  @author Mayank Sharma
    !!  @date   3/31/2017
    !!
    !-------------------------------------------------------------------------------
    type, extends(system_assembler_t),  public  :: assemble_backward_euler_t
        type(chidg_vector_t)    :: q_n
    contains
        procedure   :: assemble
    end type assemble_backward_euler_t
    !*******************************************************************************

    !>  Control the lhs update inside the nonlinear solver.
    !!
    !!  Reference:
    !!  Persson, P.-O., "High-Order Navier-Stokes Simulations using a Sparse
    !!  Line-Based Discontinuous Galerkin Method", AIAA-2012-0456
    !!
    !!  @author Nathan A. Wukie (AFRL)
    !!  @date   6/22/2017
    !!
    !--------------------------------------------------------------------------------
    type, extends(solver_controller_t), public :: backwardeuler_solver_controller_t

    contains
        procedure   :: update_lhs
    end type backwardeuler_solver_controller_t
    !********************************************************************************

contains

    !>  Initialize the backward_euler_t time integrator
    !!
    !!  Create the assembler and attach it to the time_integrator object so it can
    !!  be passed to the nonlinear solver
    !!
    !!  @author Mayank Sharma
    !!  @date   3/31/2017
    !!
    !-------------------------------------------------------------------------------
    subroutine init(self)
        class(backward_euler_t),    intent(inout)   :: self

        integer(ik)                     :: ierr
        type(assemble_backward_euler_t) :: assemble_backward_euler

        call self%set_name('Backward Euler')

        if (allocated(self%system)) deallocate(self%system)
        allocate(self%system, source=assemble_backward_euler, stat=ierr)
        if (ierr /= 0) call AllocationError

    end subroutine init
    !*******************************************************************************

    !>  Solution advancement via the backward Euler method
    !!
    !!  Given the system of partial differential equations consisting of the time-derivative of the
    !!  the solution vector and a spatial residual as
    !!
    !!  \f$ M \frac{\partial Q}{\partial t} + R(Q) = 0 \f$
    !!
    !!  the time derivative is discretized with a bacward difference
    !!
    !!  \f$ M \frac{Q^{n + 1} - Q^{n}}{\Delta t} + R(Q^{n + 1}) = 0 \f$
    !!
    !!  which yields a nonlinear equation
    !!
    !!  \f$  \frac{\Delta Q}{\Delta t}M + R(Q^{n + 1}) = 0 \f$
    !!
    !!  which results in the Newton iteration
    !!
    !!  \f$ \left(\frac{M}{\Delta t} + \frac{\partial R(Q^{m})}{\partial Q}\right) \delta Q^{m} =
    !!      -\frac{M}{\Delta t}\Delta Q^{m} - R(Q^{m}) \f$
    !!  \f$ Q^{m} = Q^{n} + \Delta Q^{m} \f$
    !!
    !!  where \f$ \delta Q^{m} = \Delta Q^{m + 1} -  \Delta Q^{m} \f$ for the mth Newton iteration
    !!
    !!  @author Mayank Sharma
    !!  @date   3/31/2017
    !!
    !--------------------------------------------------------------------------------------
    subroutine step(self,data,nonlinear_solver,linear_solver,preconditioner)
        class(backward_euler_t),                intent(inout)   :: self
        type(chidg_data_t),                     intent(inout)   :: data
        class(nonlinear_solver_t),  optional,   intent(inout)   :: nonlinear_solver
        class(linear_solver_t),     optional,   intent(inout)   :: linear_solver
        class(preconditioner_t),    optional,   intent(inout)   :: preconditioner

        class(*),   allocatable     :: assemble_type

        type(backwardeuler_solver_controller_t),    save    :: solver_controller

        select type(associate_name => self%system)
            type is (assemble_backward_euler_t)

                ! Store solution at nth time step to a separate vector
                associate_name%q_n = data%sdata%q

        end select

        !
        ! Solve assembled nonlinear system, the nonlinear update is the step in time
        ! System assembled in subroutine assemble
        !
        data%time_manager%t = data%time_manager%t + data%time_manager%dt
        !call update_grid(data)
        call data%update_grid()
        call nonlinear_solver%solve(data,self%system,linear_solver,preconditioner,solver_controller)

        !
        ! Store end residual from nonlinear solver
        !
        call self%residual_norm%push_back(nonlinear_solver%residual_norm%at(nonlinear_solver%residual_norm%size()))

    end subroutine step
    !**************************************************************************************

    !>  Assemble the system for the backward Euler equations with temporal contributions
    !!
    !!  @author Mayank Sharma
    !!  @date   3/31/2017
    !!
    !!  \f$ lhs = \frac{\partial R(Q)}{\partial Q} \f$
    !!  \f$ rhs = R(Q) \f$
    !!  \f$ M   = element mass matrix \f$
    !!
    !!  For system assembly with temporal contributions
    !!
    !!  \f$ lhs = \frac{M}{dt} + lhs \f$
    !!  \f$ rhs = \frac{M}{dt}*dq + rhs \f$
    !!
    !--------------------------------------------------------------------------------------
    subroutine assemble(self,data,differentiate,timing)
        class(assemble_backward_euler_t),   intent(inout)               :: self
        type(chidg_data_t),                 intent(inout)               :: data
        integer(ik),                        intent(in)                  :: differentiate
        real(rk),                           intent(inout),  optional    :: timing

        integer(ik)                 :: itime, idom, ielem, ifield, ierr
        real(rk)                    :: dt
        type(chidg_vector_t)        :: delta_q
        real(rk),   allocatable     :: temp_1(:), mat(:,:)
        type(element_info_t)        :: elem_info

        associate ( q   => data%sdata%q,   &
                    dq  => data%sdata%dq,  &
                    lhs => data%sdata%lhs, &
                    rhs => data%sdata%rhs)

        ! Clear data containers
        call rhs%clear()
        if (differentiate == dQ_DIFF) call lhs%clear()

        ! Get spatial update
        call update_space(data,differentiate,timing)

        ! Get no. of time levels ( = 1 for time marching) and time step
        dt    = data%time_manager%dt

        ! Compute delta q
        delta_q = q - self%q_n

        do idom = 1,data%mesh%ndomains()
            do ielem = 1,data%mesh%domain(idom)%nelem
                do itime = 1,data%mesh%domain(idom)%ntime
                    elem_info = data%mesh%get_element_info(idom,ielem)

                    do ifield = 1,data%eqnset(elem_info%eqn_ID)%prop%nprimary_fields()

                        ! Add time derivative to left-hand side
                        if (differentiate == dQ_DIFF) then
                            mat = data%mesh%domain(idom)%elems(ielem)%mass * (ONE/dt)
                            call data%sdata%lhs%scale_diagonal(mat,elem_info,ifield,itime)
                        end if

                        ! Add time derivative to right-hand side
                        temp_1 = (ONE/dt)*matmul(data%mesh%domain(idom)%elems(ielem)%mass, delta_q%get_field(elem_info,ifield,itime))
                        call rhs%add_field(temp_1,elem_info,ifield,itime)

                    end do  ! ifield

                end do  ! itime
            end do  ! ielem
        end do  ! idom

        ! Reassemble
        call data%sdata%rhs%assemble()
        if (differentiate == dQ_DIFF) call data%sdata%lhs%assemble()

        end associate

    end subroutine assemble
    !**************************************************************************************

    !>  Control algorithm for selectively updating the lhs matrix in the
    !!  nonlinear solver.
    !!
    !!  Reference:
    !!  Persson, P.-O., "High-Order Navier-Stokes Simulations using a Sparse
    !!  Line-Based Discontinuous Galerkin Method", AIAA-2012-0456
    !!
    !!  @author Nathan A. Wukie
    !!  @date   6/22/2017
    !!
    !!  @param[in]  niter               Number of newton iterations
    !!  @param[in]  residual_ratio      R_{i}/R_{i-1}
    !!
    !----------------------------------------------------------------------------
    function update_lhs(self,A,niter,residual_ratio) result(update)
        class(backwardeuler_solver_controller_t),   intent(inout)   :: self
        type(chidg_matrix_t),                       intent(in)      :: A
        integer(ik),                                intent(in)      :: niter
        real(rk),                                   intent(in)      :: residual_ratio

        integer(ik) :: update

        ! Update lhs if:
        !   1: If matrix(lhs/A) hasn't been updated before
        !   2: number of newton iterations > 10
        !   3: residual norm increases by factor of 10 (divergence)
        !   4: being forced
        if ( all(A%stamp == NO_ID)      .or. &
            (niter > 10)                .or. &
            (residual_ratio > 10._rk)   .or. &
            (self%force_update_lhs) ) then
            !update = .true.
            update = dQ_DIFF
            self%lhs_updated = .true.
        else
            !update = .false.
            update = NO_DIFF
            self%lhs_updated = .false.
        end if

        !! Store action
        !self%lhs_updated = update

        ! Turn off forced update
        self%force_update_lhs = .false.

    end function update_lhs
    !****************************************************************************

end module type_backward_euler

max_stars_repo_name: nwukie/ChiDG

max_issues_repo_name: nwukie/ChiDG

max_forks_repo_name: nwukie/ChiDG


! { dg-do run }
!
! Test the fix for PR55763 comment 9 as part of the overall fix for PR34640.
!
! Contributed by Tobias Burnus  <burnus@gcc.gnu.org>
!
  program change_field_type
    use, intrinsic :: iso_c_binding
    implicit none
    REAL(kind=c_float), POINTER :: vector_comp(:)
    TYPE, BIND(C) :: scalar_vector
       REAL(kind=c_float) :: scalar
       REAL(kind=c_float) :: vec(3)
    END TYPE
    TYPE, BIND(C) :: scalar_vector_matrix
       REAL(kind=c_float) :: scalar
       REAL(kind=c_float) :: vec(3)
       REAL(kind=c_float) :: mat(3,3)
    END TYPE
    CLASS(*), ALLOCATABLE, TARGET :: one_d_field(:)
    real, pointer :: v1(:)

    allocate(one_d_field(3), &
             source = (/ scalar_vector( 1.0, (/ -1.0, 0.0, 1.0 /) ), &
                         scalar_vector( 1.1, (/ -1.2, 0.2, 0.9 /) ), &
                         scalar_vector( 1.2, (/ -1.4, 0.4, 0.8 /) )  /) )

    call extract_vec(one_d_field, 1, 2)
    if (any (abs (vector_comp - [0.0,0.2,0.4]) .gt. 1e-4)) STOP 1
    deallocate(one_d_field)   ! v1 becomes undefined

    allocate(one_d_field(1), &
         source = (/ scalar_vector_matrix( 1.0, (/ -1.0, 0.0, 1.0 /), &
         reshape( (/ 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0 /), &
                 (/3, 3/) ) ) /) )

    call extract_vec(one_d_field, 2, 1)
    if (abs (vector_comp(1) + 1.0) > 1e-4) STOP 2
    call extract_vec(one_d_field, 2, 3)
    if (abs (vector_comp(1) - 1.0) > 1e-4) STOP 3
    deallocate(one_d_field)   ! v1 becomes undefined
  contains
    subroutine extract_vec(field, tag, ic)
        use, intrinsic :: iso_c_binding
        CLASS(*), TARGET :: field(:)
        INTEGER(kind=c_int), value :: tag, ic

        type(scalar_vector), pointer :: sv(:)
        type(scalar_vector_matrix), pointer :: svm(:)

        select type (field)
        type is (real(c_float))
          vector_comp => field
        class default
          select case (tag)
          case (1)
             sv => field
             vector_comp => sv(:)%vec(ic)
          case (2)
             svm => field
             vector_comp => svm(:)%vec(ic)
          end select
        end select
    end subroutine
  end program


text:
!!! Time-stamp: <utils.f90 10:09, 18 Oct 2005 by P Sunthar>

!_______________________________________
!  Contains general purpose utilities
!_______________________________________

!!! $Log: utils.f90,v $
!!! Revision 1.3  2005/10/18 00:08:23  pxs565
!!! With dotlockfile external file locking program
!!!
!!! Revision 1.2  2005/05/24 01:07:10  pxs565
!!! variable NTOPB
!!!
!!! Revision 1.1  2005/03/16 00:54:06  pxs565
!!! Initial revision
!!!

Subroutine Initial_position(SType,N,L0s,R,myseed)
  Use bspglocons
  Use bspglovars
  Use bspintfacs
  !Use iflport ! for portable intel intrinsic functions, here ran()
  Implicit None

  Integer, Intent (in) :: Stype
  Integer(k4b), Intent (in) :: N
  Integer(k4b), Intent (inout) :: myseed
  Real (FPREC) , Intent (in) :: L0s
  Real (FPREC) , Intent (out), Dimension(:,:) :: R

  Integer, Save :: varseed = 201235

  Integer mu, done
  Real (FPREC) , Parameter :: upbound = 2
  Real (FPREC)  modr,b2b(Ndim),rv,eqpr, rvchk, b

  Call GaussRand(Ndim*N,R,myseed)

  R(:,1) = 0

  ! intrinsic ran() modifies the seed for every call
  ! so use a separate varible seed, so that the ran_1 sequence is
  ! not altered by this variation.
  !varseed = seed

  b = L0s*L0s

!!$  do rv = 0,0.99,.01
!!$     eqpr = normeqpr(Stype,rv,b);
!!$     write (*,*) rv,eqpr
!!$  end do

  Do mu = 2,N
     b2b = R(:,mu) ! * sqrt(sigma=1), gaussian dist betwn two adj beads

     If (SType .Ne. HOOK) Then

        done = 0;
        Do While (done .Ne. 1)

           !rv = ran(varseed)
	   call Random_Number(rv) ! standards conforming IBM
           eqpr = normeqpr(Stype,rv,b);

           ! use another uniform random number to decide acceptance
           ! ie, accept only eqpr fraction at this r
           !rvchk = ran(varseed)
	   call Random_Number(rvchk) ! standards conforming IBM
           !write (*,*) rv, rvchk , eqpr
           If (rvchk .Le. eqpr/upbound) done = 1
        End Do

        modr = Sqrt(Sum(b2b*b2b))
        b2b = b2b/modr * (rv * L0s)
     End If

     R(:,mu) = R(:,mu-1) + b2b

     ! in the case of a finitely extensible spring, a gaussian distribution
     ! can lead to incorrect forces some exceptional b2b vectors, we limit
     ! this by altering the distance when preseving the random direction.
     ! and an additional random factor between 0 and 1
     ! A good guess of the factor cud b obtained more rigorously.
  End Do
End Subroutine Initial_position

Subroutine GaussRand(N,GX,seed)
!!! Returns a normalised gaussian random variable with mean zero
!!! and variance 1: f(r) = exp(-r^2/2) / sqrt(2*Pi)
!!! To obtain a dist with a given sigma use:  sqrt(sigma) * Gx

  Use bspglocons
  Use bspintfacs
  Implicit None

  Integer(k4b), Intent (in) :: N
  Integer(k4b), Intent (inout) :: seed
  Real (FPREC) , Intent (out), Dimension(N) :: GX

  Real (FPREC)  rnd(2), x,y, rsq, fac
  Integer i,n2
  Real (FPREC)  :: rtemp(N+1)

  Do i = 1,N,2
     rsq = 0.0
     Do While (rsq .Ge. 1.0  .Or.  rsq .Eq. 0.0)
        Call ran_1(0,2,rnd,seed)
        x = 2*rnd(1) -1
        y = 2*rnd(2) -1
        rsq = x*x + y*y
     End Do
     fac = Sqrt(-2*Log(rsq)/rsq);
     rtemp(i) = x*fac
     rtemp(i+1) = y*fac
  End Do

  GX = rtemp(1:N)

End Subroutine GaussRand

Subroutine ran_1(action, n, X, idum)
  Use bspglocons
  Implicit None
  Integer, Intent(in) :: action
  Integer(k4b), Intent(inout) ::idum
  Integer(k4b), Intent(in) :: n
  Real (FPREC) , Intent(inout) :: X(n)
  !---------------------------------------------------------------------c
  !     Random number generator based on procedure given in             c
  !     Numerical Recipes in Fortran 90, Chapter B7, pp. 1141-1143      c
  !     This generator is supposed to have a period of about 3.1E18.    c
  !                                                                     c
  !     The calling sequence is the same as that of RANL,               c
  !     the BLAS random number generator. The arguments are:            c
  !     n - gives the dimension of vector R                             c
  !     X - 1D array of dimension n; on exit                            c
  !         contains n unifromly distributed random numbers in (0,1)    c
  !     idum |- seeds which are updated on exit                         c
  !---------------------------------------------------------------------c

  Integer(k4b),Parameter :: IA=16807,IM=2147483647,IQ=127773,IR=2836
  Real (FPREC) , Save :: am
  Integer(k4b), Save :: ix=1, iy = -1, k, callcntr=0
  Integer(k4b) count

  Select Case (action)
  Case (RESET)
     ix = 1
     iy = -1
     k = 0
     callcntr = 0

  Case (SHOW)
     !-- pass the value of ix, which can be used to display
     X(1) = ix
     !-- Write (*,*) 'ran_1() status: ', callcntr,ix,iy,idum

  Case DEFAULT
     callcntr=callcntr+1
     If ((idum.Le.0) .Or. (iy.Lt.0)) Then
        am = Nearest(1.0,-1.0)/IM
        iy = Ior(Ieor(888889999,Abs(idum)),1)
        ix = Ieor(777755555,Abs(idum))
        idum = Abs(idum) + 1
     End If

     Do count = 1,n
        ix = Ieor(ix,Ishft(ix,13))
        ix = Ieor(ix, Ishft(ix,-17))
        ix = Ieor(ix,Ishft(ix, 5))
        k = iy/IQ
        iy = IA*(iy-k*IQ)-IR*k
        If (iy.Lt.0) iy = iy + IM
        X(count) = am*Ior(Iand(IM,Ieor(ix,iy)),1)
     End Do

  End Select

End Subroutine ran_1

Subroutine chbyshv (L, d_a, d_b, a)
  Use bspglocons
  Implicit None
  Integer L
  Real (FPREC)  d_a, d_b, a(0:MAXCHEB)
  !---------------------------------------------------------------------c
  !     This routine calculates the Chebyshev coefficients for the      c
  !     Chebyshev polynomial approximation of a square root function.   c
  !     The routine computes L+1 coefficients, which are returned in    c
  !     in the one-dimensional array a, for the approximation of the    c
  !     square root of any variable that lies in the interval           c
  !     (d_a, d_b).                                                     c
  !---------------------------------------------------------------------c

  Integer j, k
  Real (FPREC)  xks(0:L)

  !     Calculate the shift factors
  !	d_a = 2/(lmax-lmin)
  !      -d_b/d_a = (lmax+lmin)/2

  !     Calculate the collocation points
  Do k = 0,L
     xks(k) = Cos(PI*(k+0.5)/(L+1))/d_a -d_b/d_a
  End Do

  !     Calculate the Chebyshev coefficients
  a = 0.0
  Do j = 0,L
     Do k = 0,L
        a(j) = a(j)+(xks(k)**0.5)*Cos(j*(k+0.5)*PI/(L+1))
     End Do
     a(j) = 2.0/(L+1)*a(j)
  End Do
  a(0) = a(0)/2

End Subroutine chbyshv

Subroutine maxminev_fi(n, A, maxev, minev)
  Use bspglocons
  Implicit None
  Integer n
  Real (FPREC)  A(:,:,:,:), maxev, minev
  !---------------------------------------------------------------------c
  !     This routine approximates the maximum and minimum eigen values  c
  !     of a symmetric matrix nxn matrix A using Fixman's suggestion.   c
  !---------------------------------------------------------------------c
  Integer lda, incx, incy, i
  Real (FPREC)  F(n), G(n), alpha, beta
  Real (SNGL) sdot
  Real (DOBL) ddot !! Blas

  F = 1.0
  G = 0.0

  lda = n
  alpha = 1.0
  beta = 0.0
  incx = 1
  incy = 1

  Select Case (FPREC)
  Case(SNGL)
     Call ssymv('U',n,alpha,A,lda,F,incx,beta,G,incy)
     maxev = sdot(n, F, 1, G, 1)
  Case(DOBL)
     Call dsymv('U',n,alpha,A,lda,F,incx,beta,G,incy)
     maxev = ddot(n, F, 1, G, 1)
  End Select

  maxev = 2.0*maxev/n

  !Forall (i = 1:n) F(i) = (-1.0)**i
  Do i = 1,n,2
     F(i)   = -1.
     F(i+1) =  1.
  End Do

  Select Case (FPREC)
  Case(SNGL)
     Call ssymv('U',n,alpha,A,lda,F,incx,beta,G,incy)
     minev = sdot(n, F, 1, G, 1)
  Case(DOBL)
     Call dsymv('U',n,alpha,A,lda,F,incx,beta,G,incy)
     minev = ddot(n, F, 1, G, 1)
  End Select

  minev = minev/2.0/n
  !        write (*,*) maxev, minev

End Subroutine maxminev_fi

Subroutine polish_poly_root(c,xin,atol)
  Use Bspglocons
  Implicit None
  ! use newton raphson to polish the root of a polynomial
  Integer, Parameter :: n=4 ! presently only for cubic
  Real (FPREC) , Intent (in) :: c(n)
  Real (FPREC) , Intent (inout) :: xin
  Real (FPREC) , Intent (in) :: atol

  Real (FPREC)  ::  p, p1,x0,x

  Integer i,iter
  Integer, Parameter :: IMAX = 30

  x = xin

  ! algo from NR: techniques for polishing, roots of polynomials
  Do iter=1,IMAX
     x0 = x
     p =  c(n)*x + c(n-1)
     p1 = c(n)
     Do i=n-2,1,-1
        p1 = p + p1*x
        p  = c(i) + p*x
     End Do

     !if (abs(p1) < atol) then ! f' = 0 is not solvable in newt-raph
     If (Abs(p1) .Eq. 0.0) Then ! f' = 0 is not solvable in newt-raph
        !! for WLC, in this case f'' also = 0, therefore
        p1 = 6 * c(4) ! f'''
        p1 = 6*p/p1
        ! note sign(a,b) returns sgn(b) * mod(a)
        x = x - Sign(1._FPREC,p1) * Abs(p1)**(1./3.)
        !! we omit considering cases for ILC and FENE, as they
        !! donot seem to have this singularity
     Else
        x = x - p/p1
     End If

     If (Abs(p) .Le. atol) Exit
  End Do
  If (iter>IMAX) Then
     !write (5,*) 'poly-root: loop exceeded'
  End If

  xin = x
End Subroutine polish_poly_root

Subroutine numint(f,dx,nmin,nmax,nord,sumf)
  Use Bspglocons
  Implicit None

  Real (FPREC) , Intent(in), Dimension(:) :: f
  Real (FPREC) , Intent(in) :: dx
  Real (FPREC) , Intent(out) :: sumf
  Integer, Intent (in) :: nmin, nmax,nord

  Integer, Parameter :: stdout=5, stderr=6
  Integer nbound,nint,j

  nbound = Ubound(f,1)
  nint = nmax-nmin

  !write (*,*) 'nbound = ', nbound

  If (nint > nbound .Or. nmin < 1 .Or. nmax > nbound ) Then
     Write(stderr,*) 'Array out of bounds: ', nmin,nmax,nbound
     Return
  End If

  sumf = 0.

  Select Case (nord)
  Case (1) ! trapezoidal rule
     sumf = 0.5 * (f(nmax) + f(nmin))

  Case(2)
     sumf = 5./12 * (f(nmin) + f(nmax)) + &
          13./12 *(f(nmin+1) + f(nmax-1))
  Case (3)
     sumf = 3./8 * (f(nmin) + f(nmax)) + &
          7./6 *(f(nmin+1) + f(nmax-1)) + &
          23./24 *(f(nmin+2) + f(nmax-2))
  End Select

  Do j = nmin+nord, nmax-nord
     sumf = sumf + f(j)
  End Do
  sumf =  dx*sumf

End Subroutine numint

Subroutine meanerr(vec,mean,err)
  Use Bspglocons
  Real (FPREC), Intent(in) :: vec(:)
  Real (FPREC), Intent(out) :: mean
  Real (FPREC), Intent(out) :: err ! standard error of mean

  Integer n

  n = Ubound(vec,1)
  mean = Sum(vec)/n

  ! though the following is correct,
  ! err = sqrt((sum(vec*vec) - mean*mean*n)/(n-1))
  ! it does not always numerically evaluate to a quantity > 0
  ! for very small errors, such as in the case of diffusivity.
  ! so use the sure shot formula.
  err = 2*Sqrt(Sum((vec-mean)*(vec-mean))/(n-1)/n)
  ! the standard error of mean is approximately in the 95% confidence
  ! interval, giving the factor 2 (from the t-distribution)

End Subroutine meanerr

Function normeqpr(Stype, r, b)
  Use bspglocons
  Implicit None
  Integer, Intent (in) :: Stype
  Real (FPREC) , Intent (in) :: b, r
  Real (FPREC)   normeqpr

  Real (FPREC)  expphi
  External expphi

  Integer, Save :: norm_computed = 0
  Real (FPREC)    , Save :: norm_b, bsav

  Integer , Parameter :: Ln = 21 ! 21 always see below for X,W data
  Real (DOBL) X(Ln), W(Ln), scratch(Ln)
  Real (DOBL)  endpts(2)

  Integer n
  Real (FPREC)  xfact,rg, fb , M

  If ( norm_computed .Eq. 0  .And. bsav .Ne. b) Then

     !! generate the gauss abscissa and weights for legendre
     !call gaussq(1, Ln, 0, 0, 0, endpts, scratch, X, W) ;

     !! generated from c code
     X(01) = -0.993752170620389;   W(01) = 0.016017228257774
     X(02) = -0.967226838566306;   W(02) = 0.03695378977085309
     X(03) = -0.920099334150401;   W(03) = 0.05713442542685689
     X(04) = -0.853363364583318;   W(04) = 0.0761001136283793
     X(05) = -0.768439963475678;   W(05) = 0.09344442345603385
     X(06) = -0.667138804197413;   W(06) = 0.1087972991671478
     X(07) = -0.55161883588722;   W(07) = 0.1218314160537286
     X(08) = -0.424342120207439;   W(08) = 0.1322689386333376
     X(09) = -0.288021316802401;   W(09) = 0.1398873947910733
     X(10) = -0.145561854160895;   W(10) = 0.14452440398997
     X(11) = -2.4782829604619e-16;   W(11) = 0.1460811336496907
     X(12) = 0.145561854160895;   W(12) = 0.1445244039899697
     X(13) = 0.288021316802401;   W(13) = 0.1398873947910732
     X(14) = 0.424342120207439;   W(14) = 0.1322689386333371
     X(15) = 0.55161883588722;   W(15) = 0.1218314160537288
     X(16) = 0.667138804197413;   W(16) = 0.1087972991671492
     X(17) = 0.768439963475678;   W(17) = 0.09344442345603389
     X(18) = 0.853363364583317;   W(18) = 0.07610011362837897
     X(19) = 0.920099334150401;   W(19) = 0.05713442542685797
     X(20) = 0.967226838566306;   W(20) = 0.03695378977085323
     X(21) = 0.993752170620389;   W(21) = 0.01601722825777395

     M = 18  ! some large number, obtained by trial and error

     ! limit the upper bound of M
     If (M > b/2) M = b/2

     xfact = Sqrt(2*M/b)

     norm_b = 0
     fb = 0
     Do n=1,Ln
        rg = (X(n) + 1)/2 * xfact
        norm_b = norm_b +  W(n) * expphi(Stype,rg,b) * rg * rg
        fb = fb +  W(n) * expphi(Stype,rg,b) * rg**4
     End Do

     !Write (*,*) '# b = ', b
     !write (*,*) '# computed normb = ', norm_b
     !Write (*,*) '# f(b) = ', b*fb/norm_b

     bsav = b
     norm_computed = 1
  End If

  normeqpr = r*r*expphi(Stype,r,b)/norm_b

End Function normeqpr

Function expphi(Stype, r, b)
!!! return the distribution function (without normalisation)
  Use bspglocons
  Implicit None
  Integer, Intent (in) :: Stype
  Real (FPREC) , Intent (in) :: r, b
  Real (FPREC)   expphi

  Select Case (Stype)
  Case (HOOK)
     Write (*,*) 'Hookean case called in expb: Check the code'
  Case (FENE)
     expphi = (1-r*r)**(b/2)

  Case (WLC)
     expphi = Exp( -b/6.0 * ( 2*r*r + 1/(1-r) -r -1 ));
  End Select
End Function expphi

Function  lam1_th(hs, N)
  Use bspglocons
  Implicit None
  Real (FPREC)  lam1_th
  Real (FPREC) ,  Intent (in) :: hs
  Integer,  Intent (in) :: N

  Real (FPREC) s,b,aj

  b   = 1 - 1.66*hs**0.78;
  s   =   - 1.40*hs**0.78;

  aj = 4*Sin(Pi/2./N)*Sin(Pi/2./N);
  lam1_th = 2./( aj * b / N**s);

End Function lam1_th

!!!===================================================================!!!
!!! This module contains utilities to lock files to prevent other programs
!!! to open and write data into it.  This is done by creating another file in
!!! the directory called <file>.lock, where <file> is the file to be locked.
!!! The module contains routines for locking, unlocking (removing the lockfile
!!! The locking and unlocking is achieved by an external executable
!!!   dotlockfile
!!! Since an external program is called, make sure the PATH variable
!!! contains the location of dotlockfile
!!! Attempts to lock files from within this program was not entirely
!!! successfull, and difficult to find where the error occurred.
!!!===================================================================!!!

Module flock_utils
Implicit none

!--a constant number to be added to the unit of the file to be locked
Integer, parameter :: lunit = 3141
Character(40) :: lfile,syscmd

Contains

  Function lockfilename(file)
    Character (*) file
    Character (40) lockfilename
    !--string suffix for the lock file
    Character(5), parameter :: sufix = '.lock'

    lockfilename = trim(file) // sufix
  end Function lockfilename

  Subroutine islocked(file,lstat)

    Character(*), intent(in) :: file
    Logical, intent(out) :: lstat

    lfile=lockfilename(file)
    Inquire (file=lfile, exist=lstat)

  end Subroutine islocked

  Subroutine getlock(file)
  !! Obtain a lock
  !  Either a new one, or wait till the present one is released.
  !  If a stale file (not used for 5 mins), it is deleted
  !  and a fresh lock is obtained (see details in manual page
  !  of dotlockfile

    Character(*), intent(in) :: file

    Character(8) date
    Character(10) time

    lfile = lockfilename(file)

    !Call Date_and_Time(date,time)
    !write(*,*) time,': Waiting to obtain lock'
    !-- 12 retries (390 secs) or stale file (300 secs) which ever is earlier
    syscmd= 'dotlockfile -r 12 ' // lfile
    call system(syscmd)
    !Call Date_and_Time(date,time)
    !write(*,*) time,': Obtained lock'

  end Subroutine getlock

  Subroutine unlockfile(file)

    Character(*), intent(in) :: file

    lfile = lockfilename(file)
    syscmd= 'dotlockfile -u ' // lfile
    call system(syscmd)

  end Subroutine unlockfile
end Module flock_utils

max_stars_repo_name: psunthar/bds

max_issues_repo_name: psunthar/bds

max_forks_repo_name: psunthar/bds


! Copyright 2014 Lorenz Hüdepohl
!
! This file is part of ftimings.
!
! ftimings is free software: you can redistribute it and/or modify
! it under the terms of the GNU Lesser General Public License as published by
! the Free Software Foundation, either version 3 of the License, or
! (at your option) any later version.
!
! ftimings 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
! GNU Lesser General Public License for more details.
!
! You should have received a copy of the GNU Lesser General Public License
! along with ftimings.  If not, see <http://www.gnu.org/licenses/>.

!> \mainpage Ftimings
!>
!> An almost pure-fortran attempt to play with tree structures, which evolved
!> into the timing library used e.g. by the VERTEX supernova code.
!>
!> All you need to know is contained in the \ref ftimings::timer_t derived type.
module ftimings
   use ftimings_type
   use ftimings_value
   use, intrinsic :: iso_fortran_env, only : error_unit, output_unit
   implicit none
   save

   private

   ! this is mainly needed for Doxygen, they are
   ! by implicitly reachable as type-bound procedures
   ! of timer_t, however Doxygen does not document them
   ! if they are not also public
   public timer_start, timer_stop, timer_free, timer_print, &
      timer_enable, timer_disable, timer_is_enabled, &
      timer_in_entries, timer_get, timer_since, timer_sort, &
      timer_set_print_options, &
      timer_measure_flops, &
      timer_measure_allocated_memory, &
      timer_measure_virtual_memory,   &
      timer_measure_max_allocated_memory,   &
      timer_measure_memory_bandwidth

   character(len=name_length), private, parameter :: own = "(own)"
   character(len=name_length), private, parameter :: below = "(below threshold)"

   !> Type for a timer instance.
   !>
   !> Typical usage:
   !> \code{.f90}
   !>   type(timer_t) :: timer
   !>
   !>   call timer%enable()
   !>
   !>   call timer%start("section")
   !>     ...
   !>   call timer%start("subsection")
   !>     ...
   !>   call timer%stop("subsection")
   !>     ...
   !>   call timer%stop("section")
   !>
   !>   call timer%print()
   !> \endcode
   !>
   !> Every first call to timer%start() at a certain point in the graph
   !> allocates a small amount of memory. If the timer is no longer needed,
   !> all that memory can be freed again with
   !>
   !> \code{.f90}
   !>   call timer%free()
   !> \endcode
   type, public :: timer_t
      logical, private :: active = .false.                         !< If set to .false., most operations return immediately without any action
      logical, private :: record_allocated_memory = .false.        !< IF set to .true., record also the current resident set size
      logical, private :: record_virtual_memory = .false.          !< IF set to .true., record also the virtual memory
      logical, private :: record_max_allocated_memory = .false.    !< IF set to .true., record also the max resident set size ("high water mark")
      logical, private :: record_flop_counts = .false.             !< If set to .true., record also FLOP counts via PAPI calls
      logical, private :: record_memory_bandwidth = .false.        !< If set to .true., record also FLOP counts via PAPI calls

      logical, private :: print_allocated_memory = .false.
      logical, private :: print_max_allocated_memory = .false.
      logical, private :: print_virtual_memory = .false.
      logical, private :: print_flop_count = .false.
      logical, private :: print_flop_rate = .false.
      logical, private :: print_ldst = .false.
      logical, private :: print_memory_bandwidth = .false.
      logical, private :: print_ai = .false.
      integer, private :: bytes_per_ldst = 8

      type(node_t), private, pointer :: root => NULL()             !< Start of graph
      type(node_t), private, pointer :: current_node => NULL()     !< Current position in the graph
   contains
      procedure, pass :: start => timer_start
      procedure, pass :: stop => timer_stop
      procedure, pass :: free => timer_free
      procedure, pass :: print => timer_print
      procedure, pass :: enable => timer_enable
      procedure, pass :: disable => timer_disable
      procedure, pass :: is_enabled => timer_is_enabled
      procedure, pass :: measure_flops => timer_measure_flops
      procedure, pass :: measure_allocated_memory => timer_measure_allocated_memory
      procedure, pass :: measure_virtual_memory => timer_measure_virtual_memory
      procedure, pass :: measure_max_allocated_memory => timer_measure_max_allocated_memory
      procedure, pass :: measure_memory_bandwidth => timer_measure_memory_bandwidth
      procedure, pass :: set_print_options => timer_set_print_options
      procedure, pass :: in_entries => timer_in_entries
      procedure, pass :: get => timer_get
      procedure, pass :: since => timer_since
      procedure, pass :: sort => timer_sort
   end type

   ! Private type node_t, representing a graph node
   !
   type :: node_t
      character(len=name_length) :: name             ! Descriptive name, used when printing the timings
      integer :: count = 0                           ! Number of node_stop calls
      type(value_t) :: value                         ! The actual counter data, see ftimings_values.F90
      logical :: is_running = .false.                ! .true. if still running
      type(node_t), pointer :: firstChild => NULL()
      type(node_t), pointer :: lastChild => NULL()
      type(node_t), pointer :: parent => NULL()
      type(node_t), pointer :: nextSibling => NULL()
      class(timer_t), pointer :: timer
   contains
      procedure, pass :: now => node_now
      procedure, pass :: start => node_start
      procedure, pass :: stop => node_stop
      procedure, pass :: get_value => node_get_value
      procedure, pass :: new_child => node_new_child
      procedure, pass :: get_child => node_get_child
      procedure, pass :: sum_of_children => node_sum_of_children
      procedure, pass :: sum_of_children_with_name => node_sum_of_children_with_name
      procedure, pass :: sum_of_children_below => node_sum_of_children_below
      procedure, pass :: print => node_print
      procedure, pass :: print_graph => node_print_graph
      procedure, pass :: sort_children => node_sort_children
   end type

   interface
      function microseconds_since_epoch() result(us) bind(C, name="ftimings_microseconds_since_epoch")
         use, intrinsic :: iso_c_binding
         implicit none
         integer(kind=C_INT64_T) :: us
      end function
   end interface

   interface
      function resident_set_size() result(rsssize) bind(C, name="ftimings_resident_set_size")
         use, intrinsic :: iso_c_binding
         implicit none
         integer(kind=C_LONG) :: rsssize
      end function
   end interface

   interface
      function virtual_memory() result(virtualmem) bind(C, name="ftimings_virtual_memory")
         use, intrinsic :: iso_c_binding
         implicit none
         integer(kind=C_LONG) :: virtualmem
      end function
   end interface

   interface
      function max_resident_set_size() result(maxrsssize) bind(C, name="ftimings_highwater_mark")
         use, intrinsic :: iso_c_binding
         implicit none
         integer(kind=C_LONG) :: maxrsssize
      end function
   end interface

contains

   !> Activate the timer, without this, most methods are non-ops.
   !>
   subroutine timer_enable(self)
      class(timer_t), intent(inout), target :: self

      self%active = .true.
   end subroutine

   !> Call with enabled = .true. to also record amount of newly allocated memory.
   !> By default, memory usage is not recored. Call with .false. to deactivate again.
   !>
   !> This opens /proc/self/statm, parses it, and closes it agagain and is thus
   !> quite costly, use when appropriate.
   !>
   subroutine timer_measure_allocated_memory(self, enabled)
      class(timer_t), intent(inout) :: self
      logical, intent(in) :: enabled

      self%record_allocated_memory = enabled
   end subroutine

   !> Call with enabled = .true. to also record amount of newly created virtual memory.
   !> By default, memory usage is not recored. Call with .false. to deactivate again.
   !>
   !> This opens /proc/self/statm, parses it, and closes it agagain and is thus
   !> quite costly, use when appropriate.
   !>
   subroutine timer_measure_virtual_memory(self, enabled)
      class(timer_t), intent(inout) :: self
      logical, intent(in) :: enabled

      self%record_virtual_memory = enabled
   end subroutine

   !> Call with enabled = .true. to also record amount of newly increase of max.
   !> resident memory
   !> By default, memory usage is not recored. Call with .false. to deactivate again.
   !>
   !> This opens /proc/self/status, parses it, and closes it agagain and is thus
   !> quite costly, use when appropriate.
   !>
   subroutine timer_measure_max_allocated_memory(self, enabled)
      class(timer_t), intent(inout) :: self
      logical, intent(in) :: enabled

      self%record_max_allocated_memory = enabled
   end subroutine

   !> Call with enabled = .true. to also record the memory bandwidth with PAPI
   !> By default, this is not recorded. Call with .false. to deactivate again.
   !>
   subroutine timer_measure_memory_bandwidth(self, enabled)
      class(timer_t), intent(inout) :: self
      logical, intent(in) :: enabled

      if (enabled) then
         write(0,'(a)') "ftimings: not compiled with PAPI support, disabling memory bandwidth counter"
         self%record_memory_bandwidth = .false.
      else
         ! explicitly set to .false. by caller
         self%record_memory_bandwidth = .false.
      endif
   end subroutine

   !> Call with enabled = .true. to also record FLOP counts via PAPI calls.
   !> By default no FLOPS are recored. Call with .false. to deactivate again.
   !>
   subroutine timer_measure_flops(self, enabled)
      class(timer_t), intent(inout) :: self
      logical, intent(in) :: enabled

      if (enabled) then
         write(0,'(a)') "ftimings: not compiled with PAPI support, disabling FLOP counter"
         self%record_flop_counts = .false.
      else
         ! Explicitly set to .false. by caller
         self%record_flop_counts = .false.
      endif
   end subroutine

   !> Deactivate the timer
   !>
   subroutine timer_disable(self)
      class(timer_t), intent(inout), target :: self
      self%active = .false.
   end subroutine

   !> Return whether the timer is currently running
   !>
   function timer_is_enabled(self) result(is)
      class(timer_t), intent(inout), target :: self
      logical :: is
      is = self%active
   end function

   !> Control what to print on following %print calls
   !>
   !> \param     print_allocated_memory       Amount of newly allocated,
   !>                                         resident memory
   !> \param     print_virtual_memory         Amount of newly created virtual
   !>                                         memory
   !> \param     print_max_allocated_memory   Amount of new increase of max.
   !>                                         resident memory ("high water mark")
   !> \param     print_flop_count             Number of floating point operations
   !> \param     print_flop_rate              Rate of floating point operations per second
   !> \param     print_ldst                   Number of loads+stores
   !> \param     print_memory_bandwidth       Rate of loads+stores per second
   !> \param     print_ai                     Arithmetic intensity, that is number of
   !>                                         floating point operations per
   !>                                         number of load and store
   !>                                         operations (currently untested)
   !> \param     bytes_per_ldst               For calculating the AI, assume this number
   !>                                         of bytes per load or store (default: 8)
   subroutine timer_set_print_options(self, &
      print_allocated_memory, &
      print_virtual_memory, &
      print_max_allocated_memory, &
      print_flop_count, &
      print_flop_rate, &
      print_ldst, &
      print_memory_bandwidth, &
      print_ai, &
      bytes_per_ldst)
      class(timer_t), intent(inout) :: self
      logical, intent(in), optional :: &
         print_allocated_memory, &
         print_virtual_memory, &
         print_max_allocated_memory, &
         print_flop_count, &
         print_flop_rate, &
         print_ldst, &
         print_memory_bandwidth, &
         print_ai
      integer, intent(in), optional :: bytes_per_ldst

      if (present(print_allocated_memory)) then
         self%print_allocated_memory = print_allocated_memory
         if ((.not. self%record_allocated_memory) .and. self%print_allocated_memory) then
            write(0,'(a)') "ftimings: Warning: RSS size recording was disabled, expect zeros!"
         endif
      endif

      if (present(print_virtual_memory)) then
         self%print_virtual_memory = print_virtual_memory
         if ((.not. self%record_virtual_memory) .and. self%print_virtual_memory) then
            write(0,'(a)') "ftimings: Warning: Virtual memory recording was disabled, expect zeros!"
         endif
      endif

      if (present(print_max_allocated_memory)) then
         self%print_max_allocated_memory = print_max_allocated_memory
         if ((.not. self%record_max_allocated_memory) .and. self%print_max_allocated_memory) then
            write(0,'(a)') "ftimings: Warning: HWM recording was disabled, expect zeros!"
         endif
      endif

      if (present(print_flop_count)) then
         self%print_flop_count = print_flop_count
         if ((.not. self%record_flop_counts) .and. self%print_flop_count) then
            write(0,'(a)') "ftimings: Warning: FLOP counter was disabled, expect zeros!"
         endif
      endif

      if (present(print_flop_rate)) then
         self%print_flop_rate = print_flop_rate
         if ((.not. self%record_flop_counts) .and. self%print_flop_rate) then
            write(0,'(a)') "ftimings: Warning: FLOP counter was disabled, expect zeros!"
         endif
      endif

      if (present(print_ldst)) then
         self%print_ldst = print_ldst
         if ((.not. self%record_memory_bandwidth) .and. self%print_ldst) then
            write(0,'(a)') "ftimings: Warning: Load+Store counters were disabled, expect zeros!"
         endif
      endif
      if (present(print_memory_bandwidth)) then
         self%print_memory_bandwidth = print_memory_bandwidth
         if ((.not. self%record_memory_bandwidth) .and. self%print_memory_bandwidth) then
            write(0,'(a)') "ftimings: Warning: Load+Store counters were disabled, expect zeros for memory bandwidth!"
         endif
      endif

      if (present(print_ai)) then
         self%print_ai = print_ai
         if (.not. (self%record_memory_bandwidth .and. self%record_flop_counts)) then
            write(0,'(a)') "ftimings: Warning: Memory bandwidth or FLOP counters were disabled, expect invalid values for AI"
         endif
      endif

      if (present(bytes_per_ldst)) then
         self%bytes_per_ldst = bytes_per_ldst
      endif
   end subroutine

   !> Start a timing section
   !>
   !> \param name        A descriptive name
   !> \param replace     If .true. (default .false.), replace any entries at the
   !>                    current position with the same name. If .false., add the
   !>                    time to a possibly existing entry
   !>
   !> Care must be taken to balance any invocations of %start() and %stop(), e.g.
   !> the following is valid
   !>
   !> \code{.f90}
   !>   call timer%start("A")
   !>   call timer%start("B")
   !>   call timer%stop("B")
   !>   call timer%stop("A")
   !> \endcode
   !>
   !> while the following is not
   !>
   !> \code{.f90}
   !>   call timer%start("A")
   !>   call timer%start("B")
   !>   call timer%stop("A")
   !>   call timer%stop("B")
   !> \endcode
   !>
   subroutine timer_start(self, name, replace)
      class(timer_t), intent(inout), target :: self
      character(len=*), intent(in)  :: name
      logical, intent(in), optional  :: replace
      type(node_t), pointer :: node
!$    integer :: omp_get_thread_num, omp_get_num_threads, omp_get_level, omp_get_ancestor_thread_num
!$    integer :: i

      if (.not. self%active) then
         return
      endif

      ! Deal with nested parallelization
!$    do i = 0, omp_get_level()
!$       if (omp_get_ancestor_thread_num(i) > 0) then
!$          return
!$       endif
!$    end do

      !$omp master

      if (.not. associated(self%current_node)) then
         ! First call to timer_start()
         allocate(self%root)
         self%root%name = "[Root]"
         self%root%timer => self
         call self%root%start()
         nullify(self%root%firstChild)
         nullify(self%root%lastChild)
         nullify(self%root%parent)
         nullify(self%root%nextSibling)
         self%current_node => self%root
      endif

      if (string_eq(self%current_node%name, name)) then
         !$omp critical
         write(error_unit,*) "Recursion error! Printing tree so far.."
         write(error_unit,*) "Got %start(""" // trim(name) // """), while %start(""" // trim(name) // """) was still active"
!$       write(*,*) "omp_get_thread_num() = ", omp_get_thread_num()
!$       write(*,*) "omp_get_num_threads() = ", omp_get_num_threads()
!$       write(*,*) "omp_get_level() = ", omp_get_level()
!$       do i = 0, omp_get_level()
!$          write(*,*) "omp_get_ancestor_thread_num(", i, ") = ", omp_get_ancestor_thread_num(i)
!$       end do
         call self%root%print_graph(0)
         !$omp end critical
         stop "timer_start() while same timer was active"
      endif
      node => self%current_node%get_child(name)
      if (.not. associated(node)) then
         node => self%current_node%new_child(name)
      else if (present(replace)) then
         if (replace) then
            node%value = null_value
            node%count = 0
            if (associated(node%firstChild)) then
               call deallocate_node(node%firstChild)
               nullify(node%firstChild)
               nullify(node%lastChild)
            endif
         endif
      endif

      call node%start()

      self%current_node => node

      !$omp end master

   end subroutine

   !> End a timing segment, \sa timer_start
   !>
   !> \param name        The exact same name as was used for %start().
   !>                    If not provided, close the currently active region.
   !>                    If given, warns if it does not match the last %start()
   !>                    call on stderr and disables the current timer instance.
   !>
   subroutine timer_stop(self, name)
      class(timer_t), intent(inout), target :: self
      character(len=*), intent(in), optional :: name
      logical :: error
!$    integer :: omp_get_level, omp_get_ancestor_thread_num
!$    integer :: i

      if (.not. self%active) then
         return
      endif

      ! Deal with nested parallelization
!$    do i = 0, omp_get_level()
!$       if (omp_get_ancestor_thread_num(i) > 0) then
!$          return
!$       endif
!$    end do

      !$omp master
      error = .false.

      if (.not. associated(self%current_node)) then
         write(error_unit,'(a)') "Called timer_stop() without first calling any timer_start(), disabling timings"
         call self%free()
         self%active = .false.
         error = .true.
      else if (present(name)) then
         if (.not. string_eq(self%current_node%name, name)) then
            write(error_unit,'(a)') "Expected %stop(""" // trim(self%current_node%name)  // """),&
            & but got %stop(""" // trim(name) //  """), disabling timings"
            call self%free()
            self%active = .false.
            error = .true.
         endif
      endif

      if (.not. error) then
         call self%current_node%stop()

         ! climb up to parent
         if (.not. associated(self%current_node%parent)) then
            write(error_unit,'(a)') "Error: No valid parent node found for node '" // trim(self%current_node%name) // "'"
            call self%free()
            self%active = .false.
         endif
         self%current_node => self%current_node%parent

      endif
      !$omp end master

   end subroutine

   !> Deallocate all objects associated with (but not including) self
   !>
   subroutine timer_free(self)
      class(timer_t), intent(inout), target :: self
      if (associated(self%root)) then
         call deallocate_node(self%root)
      endif
      nullify(self%root)
      nullify(self%current_node)
   end subroutine

   !> Print a timing graph
   !>
   !> \param name1       If given, first descend one level to the node with name name1
   !> \param name2       If given, also descend another level to the node with name2 there
   !> \param name3       etc.
   !> \param name4       etc.
   !> \param threshold   If given, subsume any entries with a value of threshold
   !>                    seconds in a single node "(below threshold)"
   !> \param is_sorted   Assume a sorted graph for inserting "(own)" and "(below threshold)"
   !> \param unit        The unit number on which to print, default stdout
   !>
   subroutine timer_print(self, name1, name2, name3, name4, threshold, is_sorted, unit)
      class(timer_t), intent(in), target :: self
      character(len=*), intent(in), optional :: name1, name2, name3, name4
      real(kind=rk), intent(in), optional :: threshold
      logical, intent(in), optional :: is_sorted
      integer, intent(in), optional :: unit

      integer :: unit_act

      type(node_t), pointer :: node
      character(len=64) :: format_spec

      ! I hate fortran's string handling
      character(len=name_length), parameter :: group = "Group"
      character(len=12), parameter :: seconds    = "         [s]"
      character(len=12), parameter :: fract      = "    fraction"
      character(len=12), parameter :: ram        = "  alloc. RAM"
      character(len=12), parameter :: vmem       = "   alloc. VM"
      character(len=12), parameter :: hwm        = "  alloc. HWM"
      character(len=12), parameter :: flop_rate  = "     Mflop/s"
      character(len=12), parameter :: flop_count = "       Mflop"
      character(len=12), parameter :: ldst       = "loads+stores"
      character(len=12), parameter :: bandwidth  = "  mem bandw."
      character(len=12), parameter :: ai         = "arithm. Int."
      character(len=12), parameter :: dash       = "============"

      if (.not. self%active) then
         return
      endif

      if (present(unit)) then
         unit_act = unit
      else
         unit_act = output_unit
      endif

      node => self%root
      if (present(name1)) then
         node => node%get_child(name1)
         if (.not. associated(node)) then
            write(error_unit,'(a)') "Could not descend to """ // trim(name1)  // """"
            return
         endif
      end if
      if (present(name2)) then
         node => node%get_child(name2)
         if (.not. associated(node)) then
            write(error_unit,'(a)') "Could not descend to """ // trim(name2)  // """"
            return
         endif
      end if
      if (present(name3)) then
         node => node%get_child(name3)
         if (.not. associated(node)) then
            write(error_unit,'(a)') "Could not descend to """ // trim(name3)  // """"
            return
         endif
      end if
      if (present(name4)) then
         node => node%get_child(name4)
         if (.not. associated(node)) then
            write(error_unit,'(a)') "Could not descend to """ // trim(name4)  // """"
            return
         endif
      end if

      ! I really do hate it ..
      write(format_spec,'("("" /= "",a",i0,",2x,a12,1x,a12)")') name_length
      write(unit_act, format_spec, advance='no') adjustl(group), seconds, fract

      if (self%print_allocated_memory) then
         write(unit_act,'(1x,a12)',advance='no') ram
      endif

      if (self%print_virtual_memory) then
         write(unit_act,'(1x,a12)',advance='no') vmem
      endif

      if (self%print_max_allocated_memory) then
         write(unit_act,'(1x,a12)',advance='no') hwm
      endif

      if (self%print_flop_count) then
         write(unit_act,'(1x,a12)',advance='no') flop_count
      endif
      if (self%print_flop_rate) then
         write(unit_act,'(1x,a12)',advance='no') flop_rate
      endif
      if (self%print_ldst) then
         write(unit_act,'(1x,a12)',advance='no') ldst
      endif
      if (self%print_memory_bandwidth) then
         write(unit_act,'(1x,a12)',advance='no') bandwidth
      endif
      if (self%print_ai) then
         write(unit_act,'(1x,a12)',advance='no') ai
      endif

      write(unit_act,'(a)') ""

      write(format_spec,'("("" |  "",a",i0,",1x,2(1x,a12))")') name_length
      write(unit_act, format_spec, advance='no') "", dash, dash

      if (self%print_allocated_memory) then
         write(unit_act,'(1x,a12)',advance='no') dash
      endif

      if (self%print_virtual_memory) then
         write(unit_act,'(1x,a12)',advance='no') dash
      endif

      if (self%print_max_allocated_memory) then
         write(unit_act,'(1x,a12)',advance='no') dash
      endif

      if (self%print_flop_count) then
         write(unit_act,'(1x,a12)',advance='no') dash
      endif
      if (self%print_flop_rate) then
         write(unit_act,'(1x,a12)',advance='no') dash
      endif
      if (self%print_ldst) then
         write(unit_act,'(1x,a12)',advance='no') dash
      endif
      if (self%print_memory_bandwidth) then
         write(unit_act,'(1x,a12)',advance='no') dash
      endif
      if (self%print_ai) then
         write(unit_act,'(1x,a12)',advance='no') dash
      endif

      write(unit_act,'(a)') ""

      call node%print_graph(0, threshold, is_sorted, unit=unit)

   end subroutine

   !> Return the sum of all entries with a certain name below
   !> a given node. Specify the name with the last argument, the
   !> path to the starting point with the first few parameters
   !>
   !> \param name1, .., namei-1  The path to the starting node
   !> \param namei               The name of all sub-entries below this
   !>                            node which should be summed together
   !>
   !> For example timer%in_entries("foo", "bar", "parallel") returns
   !> the sum of all entries named "parallel" below the foo->bar node
   !>
   function timer_in_entries(self, name1, name2, name3, name4) result(s)
      use, intrinsic :: iso_fortran_env, only : error_unit
      class(timer_t), intent(in), target :: self
      character(len=*), intent(in) :: name1
      character(len=*), intent(in), optional :: name2, name3, name4
      real(kind=rk) :: s
      type(node_t), pointer :: node ! the starting node
      type(value_t) :: val
      character(len=name_length) :: name ! the name of the sections

      s = 0._rk

      if (.not. self%active) then
         return
      endif

      node => self%root
      name = name1

      if (present(name2)) then
         node => node%get_child(name1)
         if (.not. associated(node)) then
            write(error_unit,'(a)') "Could not descend to """ // trim(name1)  // """"
            return
         endif
         name = name2
      end if
      if (present(name3)) then
         node => node%get_child(name2)
         if (.not. associated(node)) then
            write(error_unit,'(a)') "Could not descend to """ // trim(name2)  // """"
            return
         endif
         name = name3
      end if
      if (present(name4)) then
         node => node%get_child(name3)
         if (.not. associated(node)) then
            write(error_unit,'(a)') "Could not descend to """ // trim(name3)  // """"
            return
         endif
         name = name4
      end if

      val = node%sum_of_children_with_name(name)
      s = real(val%micros, kind=rk) * 1e-6_rk
   end function

   !> Access a specific, already stopped entry of the graph by specifying the
   !> names of the nodes along the graph from the root node
   !>
   !> The result is only meaningfull if the entry was never appended by
   !> additional %start() calls.
   !>
   function timer_get(self, name1, name2, name3, name4, name5, name6) result(s)
      class(timer_t), intent(in), target :: self
      ! this is clunky, but what can you do..
      character(len=*), intent(in), optional :: name1, name2, name3, name4, name5, name6
      real(kind=rk) :: s
      type(node_t), pointer :: node

      s = 0._rk

      if (.not. self%active) then
         return
      endif

      node => self%root
      if (present(name1)) then
         node => node%get_child(name1)
         if (.not. associated(node)) then
            write(error_unit,'(a)') "Could not descend to """ // trim(name1)  // """"
            return
         endif
      end if
      if (present(name2)) then
         node => node%get_child(name2)
         if (.not. associated(node)) then
            write(error_unit,'(a)') "Could not descend to """ // trim(name2)  // """"
            return
         endif
      end if
      if (present(name3)) then
         node => node%get_child(name3)
         if (.not. associated(node)) then
            write(error_unit,'(a)') "Could not descend to """ // trim(name3)  // """"
            return
         endif
      end if
      if (present(name4)) then
         node => node%get_child(name4)
         if (.not. associated(node)) then
            write(error_unit,'(a)') "Could not descend to """ // trim(name4)  // """"
            return
         endif
      end if
      if (present(name5)) then
         node => node%get_child(name5)
         if (.not. associated(node)) then
            write(error_unit,'(a)') "Could not descend to """ // trim(name5)  // """"
            return
         endif
      end if
      if (present(name6)) then
         node => node%get_child(name6)
         if (.not. associated(node)) then
            write(error_unit,'(a)') "Could not descend to """ // trim(name6)  // """"
            return
         endif
      end if
      if (node%is_running) then
         write(error_unit,'(a)') "Timer """ // trim(node%name) // """ not yet stopped"
         return
      endif
      s = real(node%value%micros, kind=rk) * 1e-6_rk
   end function

   !> Access a specific, not yet stopped entry of the graph by specifying the
   !> names of the nodes along the graph from the root node and return the
   !> seconds that have passed since the entry was created.
   !>
   !> The result is only meaningfull if the entry was never appended by
   !> additional %start() calls.
   !>
   function timer_since(self, name1, name2, name3, name4) result(s)
      class(timer_t), intent(in), target :: self
      character(len=*), intent(in), optional :: name1, name2, name3, name4
      real(kind=rk) :: s
      type(value_t) :: val
      type(node_t), pointer :: node

      s = 0._rk

      node => self%root
      if (present(name1)) then
         node => node%get_child(name1)
         if (.not. associated(node)) then
            write(error_unit,'(a)') "Could not descend to """ // trim(name1)  // """"
            return
         endif
      end if
      if (present(name2)) then
         node => node%get_child(name2)
         if (.not. associated(node)) then
            write(error_unit,'(a)') "Could not descend to """ // trim(name2)  // """"
            return
         endif
      end if
      if (present(name3)) then
         node => node%get_child(name3)
         if (.not. associated(node)) then
            write(error_unit,'(a)') "Could not descend to """ // trim(name3)  // """"
            return
         endif
      end if
      if (present(name4)) then
         node => node%get_child(name4)
         if (.not. associated(node)) then
            write(error_unit,'(a)') "Could not descend to """ // trim(name4)  // """"
            return
         endif
      end if
      if (node%is_running .neqv. .true.) then
         write(error_unit,'(a)') "Timer """ // trim(node%name) // """ already stopped"
         return
      endif
      val = node%value + node%now()
      s = real(val%micros, kind=rk) * 1e-6_rk
   end function

   !> Sort the graph on each level.
   !> Warning: This irrevocable destroys the old ordering.
   !>
   subroutine timer_sort(self)
      class(timer_t), intent(inout), target :: self
      type(node_t), pointer :: node

      call sort_nodes(self%root, node)

      node => self%root
      do while (associated(node))
         call node%sort_children()
         node => node%nextSibling
      enddo
   end subroutine

   ! Now methods of node_t:

   ! This is the function that actually returns the current timestamp and all other counters
   function node_now(self) result(val)
      use, intrinsic :: iso_c_binding
      class(node_t), intent(in) :: self
      type(value_t) :: val

      ! current time
      val%micros = microseconds_since_epoch()

      if (self%timer%record_allocated_memory) then
         val%rsssize = resident_set_size()
      endif

      if (self%timer%record_virtual_memory) then
         val%virtualmem = virtual_memory()
      endif

      if (self%timer%record_max_allocated_memory) then
         val%maxrsssize = max_resident_set_size()
      endif

   end function

   subroutine node_start(self)
      class(node_t), intent(inout) :: self

      ! take the time
      self%value = self%value - self%now()
      self%is_running = .true.
   end subroutine

   subroutine node_stop(self)
      class(node_t), intent(inout) :: self

      self%count = self%count + 1

      ! take the time
      self%value = self%value + self%now()
      self%is_running = .false.
   end subroutine

   function node_get_value(self) result(val)
      class(node_t), intent(in) :: self
      type(value_t) :: val
      val = self%value
      if (self%is_running) then
         ! we have not finished, give time up to NOW
         val = val + self%now()
      endif
   end function

   function node_new_child(self, name) result(new)
      class(node_t), intent(inout), target :: self
      character(len=*), intent(in) :: name
      type(node_t), pointer :: new

      if (.not. associated(self%lastChild)) then
         allocate(self%lastChild)
         new => self%lastChild
         self%firstChild => new
      else
         allocate(self%lastChild%nextSibling)
         new => self%lastChild%nextSibling
         self%lastChild => new
      endif

      select type (self)
       type is (node_t)
         new%parent => self
       class default
         stop "node_new_child(): This should not happen"
      end select

      new%name = name
      new%count = 0
      new%timer => self%timer

      nullify(new%firstChild)
      nullify(new%lastChild)
      nullify(new%nextSibling)
   end function

   function string_eq(str1, str2) result(eq)
      character(len=name_length), intent(in) :: str1
      character(len=*), intent(in) :: str2
      logical :: eq
      eq = trim(str1) .eq. str2(1:min(len(trim(str2)), name_length))
   end function

   function node_get_child(self, name) result(child)
      class(node_t), intent(in) :: self
      character(len=*), intent(in) :: name
      type(node_t), pointer :: child

      child => self%firstChild
      do while (associated(child))
         if (string_eq(child%name, name)) then
            return
         endif
         child => child%nextSibling
      enddo
      nullify(child)
   end function

   recursive subroutine deallocate_node(entry)
      type(node_t), intent(inout), pointer :: entry
      type(node_t), pointer :: nextSibling

      if (associated(entry%firstChild)) then
         call deallocate_node(entry%firstChild)
      endif
      nextSibling => entry%nextSibling
      deallocate(entry)
      nullify(entry)
      if (associated(nextSibling)) then
         call deallocate_node(nextSibling)
      endif
   end subroutine

   function node_sum_of_children(self) result(sum_time)
      class(node_t), intent(in) :: self
      type(node_t), pointer :: cur_entry
      type(value_t) :: sum_time

      cur_entry => self%firstChild
      do while (associated(cur_entry))
         sum_time = sum_time + cur_entry%get_value()
         cur_entry => cur_entry%nextSibling
      enddo
   end function

   recursive function node_sum_of_children_with_name(self, name) result(sum_time)
      class(node_t), intent(in) :: self
      character(len=*), intent(in) :: name
      type(node_t), pointer :: cur_entry
      type(value_t) :: sum_time

      cur_entry => self%firstChild
      do while (associated(cur_entry))
         if (string_eq(cur_entry%name, name)) then
            sum_time = sum_time + cur_entry%value
         else
            sum_time = sum_time + cur_entry%sum_of_children_with_name(name)
         endif
         cur_entry => cur_entry%nextSibling
      enddo
   end function

   function node_sum_of_children_below(self, threshold) result(sum_time)
      class(node_t), intent(in) :: self
      real(kind=rk), intent(in), optional :: threshold
      type(node_t), pointer :: cur_entry
      type(value_t) :: sum_time, cur_value

      if (.not. present(threshold)) then
         return
      endif

      cur_entry => self%firstChild

      do while (associated(cur_entry))
         cur_value = cur_entry%get_value()
         if (cur_value%micros * 1e-6_rk < threshold) then
            sum_time = sum_time + cur_value
         endif
         cur_entry => cur_entry%nextSibling
      enddo
   end function

   subroutine insert_into_sorted_list(head, node)
      type(node_t), pointer, intent(inout) :: head
      type(node_t), target, intent(inout)    :: node
      type(node_t), pointer :: cur

      if (node%value%micros >= head%value%micros) then
         node%nextSibling => head
         head => node
         return
      endif

      cur => head
      do while (associated(cur%nextSibling))
         if (cur%value%micros > node%value%micros .and. node%value%micros >= cur%nextSibling%value%micros) then
            node%nextSibling => cur%nextSibling
            cur%nextSibling => node
            return
         endif
         cur => cur%nextSibling
      end do

      ! node has to be appended at the end
      cur%nextSibling => node
      node%nextSibling => NULL()
   end subroutine

   subroutine remove_from_list(head, node)
      type(node_t), pointer, intent(inout) :: head
      type(node_t), pointer, intent(in)    :: node
      type(node_t), pointer :: cur

      if (associated(head,node)) then
         head => head%nextSibling
         return
      endif

      cur => head
      do while (associated(cur%nextSibling))
         if (associated(cur%nextSibling,node)) then
            cur%nextSibling => cur%nextSibling%nextSibling
            return
         endif
         cur => cur%nextSibling
      end do
   end subroutine

   subroutine node_print(self, indent_level, total, unit)
      class(node_t), intent(inout) :: self
      integer, intent(in) :: indent_level
      type(value_t), intent(in) :: total
      type(value_t) :: val
      integer, intent(in) :: unit
      character(len=name_length) :: name, suffix

      if (self%is_running) then
         name = trim(self%name) // " (running)"
      else
         name = self%name
      endif

      if (self%count > 1) then
         write(suffix, '(" (",i0,"x)")') self%count
         name = trim(name) // " " // trim(suffix)
      endif

      if (self%is_running) then
         val = self%value + self%now()
      else
         val = self%value
      endif
      call print_value(val, self%timer, indent_level, name, total, unit)
   end subroutine

   recursive subroutine node_print_graph(self, indent_level, threshold, is_sorted, total, unit)
      use, intrinsic :: iso_fortran_env, only : output_unit
      class(node_t), intent(inout) :: self
      integer, intent(in) :: indent_level
      real(kind=rk), intent(in), optional :: threshold
      logical, intent(in), optional :: is_sorted
      type(value_t), intent(in), optional :: total
      integer, intent(in), optional :: unit

      type(node_t), pointer :: node
      integer :: i
      type(value_t) :: cur_value, node_value, own_value, below_threshold_value, total_act
      type(node_t), pointer :: own_node, threshold_node
      real(kind=rk) :: threshold_act
      logical :: is_sorted_act, print_own, print_threshold
      integer :: unit_act

      nullify(own_node)
      nullify(threshold_node)

      if (present(threshold)) then
         threshold_act = threshold
      else
         threshold_act = 0
      endif

      if (present(is_sorted)) then
         is_sorted_act = is_sorted
      else
         is_sorted_act = .false.
      endif

      cur_value = self%get_value()

      if (present(total)) then
         total_act = total
      else
         total_act = cur_value
      endif

      if (present(unit)) then
         unit_act = unit
      else
         unit_act = output_unit
      endif

      call self%print(indent_level, total_act, unit_act)

      own_value = cur_value - self%sum_of_children()
      below_threshold_value = self%sum_of_children_below(threshold)

      print_own = associated(self%firstChild)
      print_threshold = below_threshold_value%micros > 0

      ! Deal with "(own)" and "(below threshold)" entries
      if (is_sorted_act) then
         ! sort them in
         if (print_own) then
            ! insert an "(own)" node
            allocate(own_node)
            own_node%value = own_value
            own_node%name = own
            own_node%timer => self%timer
            call insert_into_sorted_list(self%firstChild, own_node)
         endif

         if (print_threshold) then
            ! insert a "(below threshold)" node
            allocate(threshold_node)
            threshold_node%value = below_threshold_value
            threshold_node%name = below
            threshold_node%timer => self%timer
            call insert_into_sorted_list(self%firstChild, threshold_node)
         endif

      else
         ! print them first
         if (print_own) then
            call print_value(own_value, self%timer, indent_level + 1, own, cur_value, unit_act)
         endif
         if (print_threshold) then
            call print_value(below_threshold_value, self%timer, indent_level + 1, below, cur_value, unit_act)
         endif
      endif

      ! print children
      node => self%firstChild
      do while (associated(node))
         node_value = node%get_value()
         if (node_value%micros * 1e-6_rk >= threshold_act &
            .or. associated(node, threshold_node) &
            .or. associated(node, own_node)) then
            call node%print_graph(indent_level + 1, threshold, is_sorted, cur_value, unit_act)
         endif
         node => node%nextSibling
      end do

      if (is_sorted_act) then
         ! remove inserted dummy nodes again
         if (print_own) then
            call remove_from_list(self%firstChild, own_node)
            deallocate(own_node)
         endif
         if (print_threshold) then
            call remove_from_list(self%firstChild, threshold_node)
            deallocate(threshold_node)
         endif
      endif

   end subroutine

   ! In-place sort a node_t linked list and return the first and last element,
   subroutine sort_nodes(head, tail)
      type(node_t), pointer, intent(inout) :: head, tail

      type(node_t), pointer :: p, q, e
      type(value_t) :: p_val, q_val
      integer :: insize, nmerges, psize, qsize, i

      if (.not. associated(head)) then
         nullify(tail)
         return
      endif

      insize = 1

      do while (.true.)
         p => head
         nullify(head)
         nullify(tail)
         nmerges = 0

         do while(associated(p))
            nmerges = nmerges + 1
            q => p
            psize = 0
            do i = 1, insize
               psize = psize + 1
               q => q%nextSibling
               if (.not. associated(q)) then
                  exit
               endif
            end do

            qsize = insize

            do while (psize > 0 .or. (qsize > 0 .and. associated(q)))
               if (psize == 0) then
                  e => q
                  q => q%nextSibling
                  qsize = qsize - 1

               else if (qsize == 0 .or. (.not. associated(q))) then
                  e => p;
                  p => p%nextSibling
                  psize = psize - 1
               else
                  p_val = p%get_value()
                  q_val = q%get_value()
                  if (p_val%micros >= q_val%micros) then
                     e => p
                     p => p%nextSibling
                     psize = psize - 1

                  else
                     e => q
                     q => q%nextSibling
                     qsize = qsize - 1

                  end if
               end if

               if (associated(tail)) then
                  tail%nextSibling => e
               else
                  head => e
               endif
               tail => e

            end do

            p => q

         end do

         nullify(tail%nextSibling)

         if (nmerges <= 1) then
            return
         endif

         insize = insize * 2

      end do
   end subroutine

   recursive subroutine node_sort_children(self)
      class(node_t), intent(inout) :: self
      type(node_t), pointer :: node

      call sort_nodes(self%firstChild, self%lastChild)

      node => self%firstChild
      do while (associated(node))
         call node%sort_children()
         node => node%nextSibling
      enddo
   end subroutine

   subroutine print_value(value, timer, indent_level, label, total, unit)
      type(value_t), intent(in) :: value
      type(timer_t), intent(in) :: timer
      integer, intent(in) :: indent_level
      character(len=name_length), intent(in) :: label
      type(value_t), intent(in) :: total
      integer, intent(in) :: unit

      character(len=64) :: format_spec

      write(format_spec,'("(",i0,"x,""|_ "",a",i0,",2x,f12.6,1x,f12.3)")') indent_level * 2 + 1, name_length
      write(unit,format_spec,advance='no') &
         label, &
         real(value%micros, kind=rk) * 1e-6_rk, &
         real(value%micros, kind=rk) / real(total%micros, kind=rk)

      if (timer%print_allocated_memory) then
         write(unit,'(1x,a12)',advance='no') &
            nice_format(real(value%rsssize, kind=C_DOUBLE))
      endif

      if (timer%print_virtual_memory) then
         write(unit,'(1x,a12)',advance='no') &
            nice_format(real(value%virtualmem, kind=C_DOUBLE))
      endif

      if (timer%print_max_allocated_memory) then
         write(unit,'(1x,a12)',advance='no') &
            nice_format(real(value%maxrsssize, kind=C_DOUBLE))
      endif

      if (timer%print_flop_count) then
         write(unit,'(1x,f12.2)',advance='no') real(value%flop_count, kind=rk) / 1e6_rk
      endif
      if (timer%print_flop_rate) then
         write(unit,'(1x,f12.2)',advance='no') real(value%flop_count, kind=rk) / value%micros
      endif
      if (timer%print_ldst) then
         write(unit,'(1x,a12)',advance='no') nice_format(real(value%ldst, kind=rk))
      endif
      if (timer%print_memory_bandwidth) then
         write(unit,'(1x,a12)',advance='no') nice_format(real(value%ldst*timer%bytes_per_ldst, kind=rk) / (value%micros * 1e-6_rk))
      endif
      if (timer%print_ai) then
         write(unit,'(1x,f12.4)',advance='no') real(value%flop_count, kind=rk) / value%ldst / timer%bytes_per_ldst
      endif

      write(unit,'(a)') ""
   end subroutine

   pure elemental function nice_format(number) result(string)
      real(kind=C_DOUBLE), intent(in) :: number
      character(len=12) :: string
      real(kind=C_DOUBLE), parameter :: &
         kibi = 2.0_C_DOUBLE**10, &
         mebi = 2.0_C_DOUBLE**20, &
         gibi = 2.0_C_DOUBLE**30, &
         tebi = 2.0_C_DOUBLE**40, &
         pebi = 2.0_C_DOUBLE**50

      if (abs(number) >= pebi) then
         write(string,'(es12.2)') number
      else if (abs(number) >= tebi) then
         write(string,'(f9.2,'' Ti'')') number / tebi
      else if (abs(number) >= gibi) then
         write(string,'(f9.2,'' Gi'')') number / gibi
      else if (abs(number) >= mebi) then
         write(string,'(f9.2,'' Mi'')') number / mebi
      else if (abs(number) >= kibi) then
         write(string,'(f9.2,'' ki'')') number / kibi
      else
         write(string,'(f12.2)') number
      endif
   end function

end module