!*****************************************************************************************
!> author: Jacob Williams
!
!  Gooding's Kepler and universal elements conversion routines.
!
!# Notes
!  The Gooding universal elements are:
!
!   * `alpha` - mu/a [km^2/s^2]
!   * `rp`    - periapsis radius [km]
!   * `inc`   - inclination [rad]
!   * `raan`  - right ascension of the ascending node [rad]
!   * `w`     - argument of periapsis [rad]
!   * `tau`   - time since last periapsis passage [sec]
!
!# References
!  1. A. W. Odell, R. H. Gooding, "Procedures for solving Kepler's equation"
!     Celestial Mechanics 38 (1986), 307-334.
!  2. R. H. Gooding, "On universal elements, and conversion procedures
!     to and from position and velocity"
!     Celestial Mechanics 44 (1988), 283-298.
!  3. R. H. Gooding, A. W. Odell. "The hyperbolic Kepler equation
!     (and the elliptic equation revisited)"
!     Celestial Mechanics 44 (1988), 267-282.

    module gooding_module

    use kind_module,    only: wp
    use numbers_module

    implicit none

    private

    !constants:
    real(wp),parameter :: ntwo   = -two
    real(wp),parameter :: pineg  = -pi
    real(wp),parameter :: half   = 0.5_wp
    real(wp),parameter :: halfpi = pi/two
    real(wp),parameter :: athird = one/three
    real(wp),parameter :: asixth = one/six

    public :: els3pv,pv3els
    public :: ekepl,ekepl1,ekepl2,emkepl,emkep,shkepl,shmkep
    public :: propagate

    contains
!*****************************************************************************************

!*****************************************************************************************
!> author: Jacob Williams
!
!  Basic two-body propagator using the Gooding universal element routines.

    pure subroutine propagate(mu, rv0, dt, rvf)

    implicit none

    real(wp),intent(in)               :: mu    !! grav. parameter [km^3/s^2]
    real(wp),dimension(6),intent(in)  :: rv0   !! initial state [km, km/s]
    real(wp),intent(in)               :: dt    !! time step [sec]
    real(wp),dimension(6),intent(out) :: rvf   !! final state [km, km/s]

    real(wp),dimension(6) :: e

    !convert to elements, increment time,
    ! then convert back to cartesian:

    call pv3els(mu, rv0, e)
    e(6) = e(6) + dt
    call els3pv(mu, e, rvf)

    end subroutine propagate
!*****************************************************************************************

!*****************************************************************************************
!>
!  Kepler's equation, `em = ekepl - (1 - e1)*sin(ekepl)`,
!  with `e1` in range 1 to 0 inclusive, solved accurately
!  (based on ekepl3, but entering `e1`, not `e`)

    pure function ekepl(em, e1)

    implicit none

    real(wp) :: ekepl
    real(wp),intent(in) :: em
    real(wp),intent(in) :: e1

    real(wp) :: emr,ee,e,w,fdd,fddd,f,fd,dee
    integer :: iter

    real(wp),parameter :: sw = 0.25_wp

    !range-reduce em to lie in range -pi to pi
    emr = mod(em,twopi)
    if (emr<pineg) emr = emr + twopi
    if (emr>pi) emr = emr - twopi
    ee = emr

    if (ee/=zero) then

        if (ee<zero) ee = -ee

        !(emr is range-reduced em & ee is absolute value of emr)
        !starter by first solving cubic equation
        e = one - e1
         w = dcbsol(e,two*e1, three*ee)

        !effectively interpolate in emr (absolute value)
        ee = (ee*ee + (pi - ee)*w)/pi
        if (emr<zero) ee = -ee

        !do two iterations of halley, each followed by newton
        do iter=1,2
            fdd = e*sin(ee)
            fddd = e*cos(ee)
            if (ee*ee/six + e1 >= sw) then
                f = (ee - fdd) - emr
                fd = one - fddd
            else
                f = emkep(e1,ee) - emr
                fd = two*e*sin(half*ee)**2 + e1
            end if
            dee = f*fd/(half*f*fdd - fd*fd)
            f = f + dee*(fd + half*dee*(fdd + athird*dee*fddd))
            !to reduce the danger of underflow replace the last line by
            !    w = fd + half*dee*(fdd + athird*dee*fddd)
            fd = fd + dee*(fdd + half*dee*fddd)
            ee = ee + dee - f/fd
            !if replacing as above, then also replace the last line by
            !ee = ee - (f - dee*(fd - w))/fd
        end do

    end if

    !range-expand
    ekepl = ee + (em - emr)

    end function ekepl
!*****************************************************************************************

!*****************************************************************************************
!>
!  Solve kepler's equation, `em = ekepl - e*sin(ekepl)`,
!  with legendre-based starter and halley iterator
!  (function has also been used under the name eafkep)

    pure function ekepl1(em, e)

    implicit none

    real(wp) :: ekepl1
    real(wp),intent(in) :: em
    real(wp),intent(in) :: e

    real(wp) :: c,s,psi,xi,eta,fd,fdd,f

    real(wp),parameter :: testsq = 1.0e-8_wp

    c = e*cos(em)
    s = e*sin(em)
    psi = s/sqrt(one - c - c + e*e)

    do

         xi = cos(psi)
         eta = sin(psi)
         fd = (one - c*xi) + s*eta
         fdd = c*eta + s*xi
         f = psi - fdd
         psi = psi - f*fd/(fd*fd - half*f*fdd)
         if (f*f < testsq) exit

     end do

     ekepl1 = em + psi

    end function ekepl1
!*****************************************************************************************

!*****************************************************************************************
!>
!  Kepler's equation, `em = ekepl - e*sin(ekepl)` with
!  e in range 0 to 1 inclusive, solved accurately

    pure function ekepl2(em, e)

    implicit none

    real(wp)            :: ekepl2
    real(wp),intent(in) :: em
    real(wp),intent(in) :: e

    real(wp) :: emr, ee, w, e1, fdd, fddd, f, fd, dee
    logical :: l
    integer :: iter

    real(wp),parameter :: sw = 0.1_wp
    real(wp),parameter :: a  = (pi-one)**2/(pi+two/three)
    real(wp),parameter :: b  = two*(pi-asixth)**2/(pi+two/three)

    !range-reduce em to line in range -pi to pi
    emr = mod(em,twopi)
    if (emr<pineg) emr = emr + twopi
    if (emr>pi) emr = emr - twopi
    ee = emr

    if (ee/=zero) then

        if (ee<zero) ee = -ee

        !(emr is range-reduced em & ee is absolute value of emr)
        !started for e = 1 by cube root of bilinear function
        if (ee<asixth) then
            ee = (six*ee)**athird
        else
            w = pi - ee
            ee = pi - a*w/(b - w)
        end if
        if (emr<zero) ee = -ee

        !interpolate for e
        ee = emr + (ee - emr)*e

        !do two iterations of halley, each followed by newton
        e1 = one - e
        l = (e1 + ee*ee/six) >= sw
        do iter=1,2
            fdd = e*sin(ee)
            fddd = e*cos(ee)
            if (l) then
                f = (ee - fdd) - emr
                fd = one - fddd
            else
                f = emkepl(e,ee) - emr
                fd = e1 + two*e*sin(half*ee)**2
            end if
            dee = f*fd/(half*f*fdd - fd*fd)
            f = f + dee*(fd + half*dee*(fdd + athird*dee*fddd))
            !to reduce the danger of underflow replace the last line by
            !w = fd + half*dee*(fdd + athird*dee*fddd)
            fd = fd + dee*(fdd + half*dee*fddd)
            ee = ee + dee - f/fd
            !if replacing as above, then also replace the last line by
            !ee = ee - (f - dee*(fd - w))/fd
        end do

    end if

    !range-expand
    ekepl2 = ee + (em - emr)

    end function ekepl2
!*****************************************************************************************

!*****************************************************************************************
!>
!  Accurate computation of `ee - e*sin(ee)`
!  when (e, ee) is close to (1, 0)
!
!@note must not be used for large ee (absolute)
!      as then rounding worse not better

    pure function emkepl(e, ee)

    implicit none

    real(wp)            :: emkepl
    real(wp),intent(in) :: e
    real(wp),intent(in) :: ee

    real(wp) :: x, ee2, term, d, x0

    x    = (one - e)*sin(ee)
    ee2  = -ee*ee
    term = ee
    d    = zero

    do

        d = d + two
        term = term*ee2/(d*(d + one))
        x0 = x
        x = x - term
        if (x==x0) exit

    end do

    emkepl = x

    end function emkepl
!*****************************************************************************************

!*****************************************************************************************
!>
!  Similar to emkepl, except input is `1-e`.

    pure function emkep(e1,ee)

    implicit none

    real(wp)            :: emkep
    real(wp),intent(in) :: e1
    real(wp),intent(in) :: ee

    real(wp) :: x, ee2, term, d, x0

    x    = e1*sin(ee)
    ee2  = -ee*ee
    term = ee
    d    = zero

    do

        d = d + two
        term = term*ee2/(d*(d + one))
        x0 = x
        x = x - term
        if (x==x0) exit

    end do

    emkep = x

    end function emkep
!*****************************************************************************************

!*****************************************************************************************
!>
!  Equation `el = shkepl + (g1 - 1)*asinh(shkepl)`,
!  with g1 in range 0 to 1 inclusive, solved accurately.

    pure function shkepl (el, g1)

    implicit none

    real(wp)            :: shkepl
    real(wp),intent(in) :: el
    real(wp),intent(in) :: g1

    real(wp) :: s,g,cl,al,w,s0,s1,s2,s3,fdd,fddd,f,fd,ds,stemp
    integer  :: iter

    real(wp),parameter :: sw=half

    s = el

    if (el/=zero) then

        !started based on lagrange's theorem
        g = one - g1
        cl = sqrt(one + el**2)
        al = asinh(el)
        w = g**2*al/cl**3
        s = one - g/cl
        s = el + g*al/dcubrt(s**3 + w*el*(1.5_wp - g/0.75_wp))

        !two iterations (at most) of halley-then-newton process
        do iter=1,2
            s0 = s*s
            s1 = s0 + one
            s2 = sqrt(s1)
            s3 = s1*s2
            fdd = g*s/s3
            fddd = g*(one - two*s0)/(s1*s3)
            if (asixth*s0 + g1 >= sw) then
                f = (s - g*asinh(s)) - el
                fd = one - g/s2
            else
                f = shmkep(g1, s) - el
                fd = (s0/(s2 + one) + g1)/s2
            end if
            ds = f*fd/(half*f*fdd - fd*fd)
            stemp = s + ds
            if (stemp==s) exit
            f = f + ds*(fd + half*ds*(fdd + athird*ds*fddd))
            fd = fd + ds*(fdd + half*ds*fddd)
             s = stemp - f/fd
        end do

    end if

    shkepl = s

     end function shkepl
!*****************************************************************************************

!*****************************************************************************************
!>
!  Accurate computation of `s - (1 - g1)*asinh(s)`
!  when (g1, s) is close to (0, 0)

    pure function shmkep (g1, s)

    implicit none

    real(wp)            :: shmkep
    real(wp),intent(in) :: g1
    real(wp),intent(in) :: s

    real(wp) :: g,t,tsq,x,term,twoi1,x0

    g = one - g1
    t = s/(one + sqrt(one + s*s))
    tsq = t*t
    x = s*(g1 + g*tsq)
    term = two*g*t
    twoi1 = one

    do

        twoi1 = twoi1 + two
        term = term*tsq
        x0 = x
        x = x - term/twoi1
        if (x==x0) exit

     end do

     shmkep = x

     end function shmkep
!*****************************************************************************************

!*****************************************************************************************
!>
!  Algorithm for two-dimensional conversion
!  from orbital elements to position and velocity.

    pure subroutine els2pv (gm, al, q, om, tau, r, u, vr, vt)

    implicit none

    real(wp),intent(in)  :: gm    !! grav. parameter [km^3/s^2]
    real(wp),intent(in)  :: al    !! alpha [km^2/s^2]
    real(wp),intent(in)  :: q     !! periapsis distance [km]
    real(wp),intent(in)  :: om    !! argument of periapsis relative to assumed reference direction [rad]
    real(wp),intent(in)  :: tau   !! time from periapsis [sec]
    real(wp),intent(out) :: r     !! radial distance [km]
    real(wp),intent(out) :: u     !! angle from reference direction [rad]
    real(wp),intent(out) :: vr    !! radial velocity [km/2]
    real(wp),intent(out) :: vt    !! transverse velocity >=0 [km/s]

    real(wp) :: d,h,v,e1,e,ep1,alp,rtal,em,ee2,s2,c2,emv,s,c

    if (al==zero) then

        !(parabola - gm cannot be zero)

        d = dcbsol(half/gm, q, 1.5_wp*gm*tau)
        r = q + half*d*d/gm
        h = sqrt(two*gm*q)
        v = two*atan2(d,h)

    else

        !(ellipse or hyperbola)

        e1 = al*q
        e = gm - e1
        ep1 = gm + e
        h = sqrt(q*ep1)
        alp = abs(al)
        rtal = sqrt(alp)
        !(last 6 items could be saved if repeating gm, al & q)

        em = tau*alp*rtal
        if (al>zero) then

            !(ellipse - gm cannot be zero)
            ! make sure e1 argument to ekepl is between [0,1]
            ee2 = half*ekepl(em/gm, max(zero,min(one,e1/gm)))
            s2 = sin(ee2)
            c2 = cos(ee2)
            r = q + two*e*s2*s2/al
            d = two*e*s2*c2/rtal
            v = two*atan2(ep1*s2, h*rtal*c2)
            emv = em/gm - v
            v = v + fourpi*sign(real(int(abs(emv/fourpi) + half),wp), emv)

        else

            !(hyperbola)
            s = shkepl(em/e, -e1/e)
            s2 = s*s
            c = sqrt(one + s2)
            s2 = s2/(c+one)
            r = q - e*s2/al
            d = e*s/rtal
            v = atan2(s*h*rtal, -gm*s2 - e1)

        end if

    end if

    !(all orbits)
    u = om + v
    vr = d/r
    vt = h/r

    end subroutine els2pv
!*****************************************************************************************

!*****************************************************************************************
!>
!  Algorithm for three-dimensional conversion
!  from orbital elements to position and velocity

    pure subroutine els3pv (gm, e, pv)

    implicit none

    real(wp),intent(in)               :: gm   !! grav. parameter [km^3/sec^2]
    real(wp),dimension(6),intent(in)  :: e    !! [al, q, ei, bom, om, tau]
    real(wp),dimension(6),intent(out) :: pv   !! [x, y, z, xdot, ydot, zdot]

    real(wp) :: x,y,z,xdot,ydot,zdot,al,q,ei,bom,om,tau
    real(wp) :: r,u,vr,vt,c,s,x1,x2,y1,y2

    if (all(e==zero)) then

        pv = zero

    else

        al  = e(1)
        q   = e(2)
        ei  = e(3)
        bom = e(4)
        om  = e(5)
        tau = e(6)

        call els2pv (gm, al, q, om, tau, r, u, vr, vt)

        c = cos(u)
        s = sin(u)
        x1 = r*c
        y1 = r*s
        x2 = vr*c - vt*s
        y2 = vr*s + vt*c
        c = cos(ei)
        s = sin(ei)
        z = y1*s
        y1 = y1*c
        zdot = y2*s
        y2 = y2*c
        c = cos(bom)
        s = sin(bom)
        x = x1*c - y1*s
        y = x1*s + y1*c
        xdot = x2*c - y2*s
        ydot = x2*s + y2*c

        pv(1) = x
        pv(2) = y
        pv(3) = z
        pv(4) = xdot
        pv(5) = ydot
        pv(6) = zdot

    end if

    end subroutine els3pv
!*****************************************************************************************

!*****************************************************************************************
!>
!  Algorithm for two-dimensional conversion
!  from position and velocity to orbital elements.

    pure subroutine pv2els (gm, r, u, vr, vt, al, q, om, tau)

    implicit none

    real(wp),intent(in)  :: gm    !! grav. parameter [km^3/s^2]
    real(wp),intent(in)  :: r     !! radial distance [km]
    real(wp),intent(in)  :: u     !! angle from assumed reference direction [rad]
    real(wp),intent(in)  :: vr    !! radial velocity [km/2]
    real(wp),intent(in)  :: vt    !! transverse velocity >=0 [km/s]
    real(wp),intent(out) :: al    !! alpha: gm/a [km^2/s^2]
    real(wp),intent(out) :: q     !! periapsis distance [km]
    real(wp),intent(out) :: om    !! argument of periapsis relative to reference direction [rad]
    real(wp),intent(out) :: tau   !! time from periapsis [sec]

    real(wp) :: esq1,es,eses,ec,ecec,esq,e,v,e1
    real(wp) :: eh,em,ecesq,en,adj,vsq,rtal,d,h,p,alp

    real(wp),parameter   :: sw = 0.25_wp
    logical,parameter    :: l = .false.

    !(all orbits)
    vsq = vr*vr + vt*vt
    al = two*gm/r - vsq
    alp = abs(al)
    rtal = sqrt(alp)
    d = r*vr
    h = r*vt
    p = h*h
    esq1 = p*al
    es = d*rtal
    eses = es*es
    ec = r*vsq - gm
    ecec = ec*ec
    if (al>zero) then
        !(one esq formula superior for the ellipse)
        esq = ecec + eses
    else
        !(different formula superior for the hyperbola)
        esq = gm*gm - esq1
    end if
    e = sqrt(esq)
    q = p/(gm + e)
    if (al==zero) then
        !(parabola)
        tau = d*(two*q + r)/(three*gm)
        v = two*atan2(vr, vt)
    else if (e==zero) then
        !(circle)
        tau = zero
        v = zero
    else
        !(ellipse or hyperbola)
        e1 = al*q
        if (al>zero) then
            !(ellipse)
            eh = atan2(es, ec)
            if (gm*eh*eh/six + e1 >= gm*sw) then
                !(general case)
                em = gm*eh - es
                ecesq = gm*ec - esq
            else
                !(for e1 & eh both near zero)
                em = gm*emkep (e1/gm, eh)
                ecesq = (esq1*ecec - esq*eses)/(esq + gm*ec)
            end if
        else
            !(hyperbola)
            eh = asinh(es/e)
            if (gm*eh*eh/six - e1 >= gm*sw) then
                !(general case)
                em = es - gm*eh
                ecesq = esq - gm*ec
            else
                !(for e1 & eh both near zero)
                em = e*shmkep(-e1/e, es/e)
                ecesq = -(esq1*ecec + esq*eses)/(esq + gm*ec)
            end if
        end if
            !(ellipse or hyperbola still)
            en = alp*rtal
            tau = em/en
            v = atan2(es*h*rtal, ecesq)
    end if

    !(all orbits)
    om = u - v

    !
    !  note: the following is never executed... set l=true and test...
    !

    if (l .and. al>zero) then
        !(for ellipse, adjust revolutions if required (using l))
        adj = twopi*sign(real(int(abs(om/twopi) + half),wp), om)
        om = om - adj
        tau = tau + adj/en
    end if

    end subroutine pv2els
!*****************************************************************************************

!*****************************************************************************************
!>
!  Algorithm for three-dimensional conversion
!  from position and velocity to orbital elements.

    pure subroutine pv3els (gm, pv, e)

    implicit none

    real(wp),intent(in)               :: gm   !! grav. parameter [km^3/s^2]
    real(wp),dimension(6),intent(in)  :: pv   !! [x, y, z, xdot, ydot, zdot]
    real(wp),dimension(6),intent(out) :: e    !! [al, q, ei, bom, om, tau]

    real(wp) :: x,y,z,xdot,ydot,zdot,al,q,ei,bom,om,tau,xsqysq,&
                rsq,r,vr,hx,hy,hz,hsq,u,vt,bx,by,bz,w,h

    if (all(pv==zero)) then

        e = zero

    else

        x    = pv(1)
        y    = pv(2)
        z    = pv(3)
        xdot = pv(4)
        ydot = pv(5)
        zdot = pv(6)

        xsqysq = x*x + y*y
        rsq = xsqysq + z*z
        r = sqrt(rsq)
        vr = (x*xdot + y*ydot + z*zdot)/r
        hx = y*zdot - z*ydot
        hy = z*xdot - x*zdot
        hz = x*ydot - y*xdot
        hsq = hx*hx + hy*hy + hz*hz
        if (hsq==zero) then
            !(rectilinear orbit)
            ei = halfpi
            if (xsqysq==zero) then
                !(axial orbit)
                bom = zero
            else
                !(general rectilinear orbit)
                bom = atan2(y, x)
            end if
            u = atan2(z, sqrt(xsqysq))
            vt = zero
        else
            !(non-degenerate orbit)
            bx = hy*z - hz*y
            by = hz*x - hx*z
            bz = hx*y - hy*x
            hx = y*bz - z*by
            hy = z*bx - x*bz
            hz = x*by - y*bx
            w = hx*hx + hy*hy
            h = sqrt(w + hz*hz)
            ei = atan2(sqrt(w), hz)
            if (w==zero) then
                !(orbit in reference plane)
                bom = zero
                u = atan2(y*sign(one,hz), x)
            else
                !(general orbit)
                bom = atan2(hx, -hy)
                u = atan2(h*z, rsq*bz)
            end if
            vt = h/(r*rsq)
        end if
        call pv2els (gm, r, u, vr, vt, al, q, om, tau)

        e(1) = al
        e(2) = q
        e(3) = ei
        e(4) = bom
        e(5) = om
        e(6) = tau

    end if

    end subroutine pv3els
!*****************************************************************************************

!*****************************************************************************************
!>
!  Solution to `a*x**3 + 3*b*x - 2c = 0`, where
!  `a` and `b**3 + a*c**2` are both non-negative
!  (zero generated, in lieu of infinity, if `a = b = 0`)

    pure function dcbsol (a, b, c) result(x)

    implicit none

    real(wp)               :: x
    real(wp),intent(in)    :: a
    real(wp),intent(in)    :: b
    real(wp),intent(in)    :: c

    real(wp) :: bsq,d

    if (a==zero .and. b==zero .or. c==zero) then
        x = zero
    else
        bsq = b*b
        d = sqrt(a) * abs(c)
        d = dcubrt(d + sqrt(b*bsq + d*d))**2
        x = two * c / (d + b + bsq / d)
    end if

    end function dcbsol
!*****************************************************************************************

!*****************************************************************************************
!>
!  Cube root computed accurately, by incorporating
!  one Newton-Raphson iteration.

    pure function dcubrt(x) result(c)

    implicit none

    real(wp)            :: c
    real(wp),intent(in) :: x

    real(wp) :: y

    if (x==zero) then
        c = zero
    else
        y = abs(x)
        c = y**athird
        c = c - athird*(c - y/c**2)
        c = sign(c,x)
    end if

    end function dcubrt
!*****************************************************************************************

!*****************************************************************************************
    end module gooding_module
!*****************************************************************************************


!=============================================================================!
        subroutine csolve1p( nsys, np, nblk, mat, rhs, per, per2, mult, fact)
!=============================================================================!
!
!  Solves a block tridiagonal system of equations without pivoting.
!  Vectorized over the second index.
!
!  This routine solves a PERIODIC block tridiagonal system.
!
!  WARNING: Since pivoting is not done, this algorithm could fail
!
!  input:
!    nsys       : number of systems to vectorize over
!    np         : number of rows in the systems including periodic node
!    nblk       : dimension of the blocks in each row
!
!  inout:
!    mat        : the block tridiagonal matrix to be solved
!    rhs        : the right hand side
!=============================================================================!
        implicit none

        integer, intent(in) :: np, nblk, nsys
        complex, intent(inout) :: mat(3,nblk,nblk,nsys,np), rhs(nblk,nsys,np)
        complex, intent(inout) :: per(nblk,nblk,nsys,np)
        complex, intent(inout) :: per2(2,nblk,nblk,nsys)
        complex, intent(inout) :: mult(nsys), fact(nsys)

!.... useful constants

        real, parameter :: zero = 0.0, one = 1.0

!.... local variables

        integer :: nc, nr, i, l, m, n, p, q, r, iv
!=============================================================================!

c!.... don't solve for the redundant nodes

        n = np - 1

        !$omp parallel do private(nc,i,l,m,p,q,r,iv)
        do iv = 1, nsys

!.... setup the periodic fillin vector

        per(:,1:nblk,iv,1)  = mat(1,:,1:nblk,iv,1)
        per2(1,:,1:nblk,iv) = mat(3,:,1:nblk,iv,n)

!=============================================================================!
!                     F O R W A R D   E L I M I N A T I O N
!=============================================================================!
        do i = 1, n - 3

!.... do the first four columns

          do nc = 1, nblk - 1

            do m = nc, nblk - 1
              l = m + 1
              mult(iv) = -mat(2,l,nc,iv,i) / mat(2,nc,nc,iv,i)
              do q = nc + 1, nblk
                mat(2,l,q,iv,i) = mat(2,l,q,iv,i) + mult(iv) * mat(2,nc,q,iv,i)
              end do
              do q = 1, nblk
                mat(3,l,q,iv,i) = mat(3,l,q,iv,i) + mult(iv) * mat(3,nc,q,iv,i)
                per(l,q,iv,i) = per(l,q,iv,i) + mult(iv) * per(nc,q,iv,i)
              end do
              rhs(l,iv,i) = rhs(l,iv,i) + mult(iv) * rhs(nc,iv,i)
            end do

          end do        ! end of loop on nc

!.... for all columns

          nc = 1

          p = i + 1

          fact(iv) = one / mat(2,nc,nc,iv,i)
          do m = 1, nblk
            mult(iv) = -mat(1,m,nc,iv,p) * fact(iv)
            do q = nc + 1, nblk
              mat(1,m,q,iv,p) = mat(1,m,q,iv,p) + mult(iv) * mat(2,nc,q,iv,i)
            end do
            do q = 1, nblk
              mat(2,m,q,iv,p) = mat(2,m,q,iv,p) + mult(iv) * mat(3,nc,q,iv,i)
              per(m,q,iv,p) = mult(iv) * per(nc,q,iv,i)
            end do
            rhs(m,iv,p) = rhs(m,iv,p) + mult(iv) * rhs(nc,iv,i)

!.... the periodic stuff on the last row

            mult(iv) = -per2(1,m,nc,iv) * fact(iv)
            do q = nc + 1, nblk
              per2(1,m,q,iv) = per2(1,m,q,iv) + mult(iv) * mat(2,nc,q,iv,i)
            end do
            do q = 1, nblk
              per2(2,m,q,iv) = mult(iv) * mat(3,nc,q,iv,i)
              mat(2,m,q,iv,n) = mat(2,m,q,iv,n) + mult(iv) * per(nc,q,iv,i)
            end do
            rhs(m,iv,n) = rhs(m,iv,n) + mult(iv) * rhs(nc,iv,i)
          end do

          do nc = 2, nblk

            fact(iv) = one / mat(2,nc,nc,iv,i)
            do m = 1, nblk
              mult(iv) = -mat(1,m,nc,iv,p) * fact(iv)
              do q = nc + 1, nblk
                mat(1,m,q,iv,p) = mat(1,m,q,iv,p) + mult(iv) * mat(2,nc,q,iv,i)
              end do
              do q = 1, nblk
                mat(2,m,q,iv,p) = mat(2,m,q,iv,p) + mult(iv) * mat(3,nc,q,iv,i)
                per(m,q,iv,p) = per(m,q,iv,p) + mult(iv) * per(nc,q,iv,i)
              end do
              rhs(m,iv,p) = rhs(m,iv,p) + mult(iv) * rhs(nc,iv,i)

!.... the periodic stuff on the last row

              mult(iv) = -per2(1,m,nc,iv) * fact(iv)
              do q = nc + 1, nblk
                per2(1,m,q,iv) = per2(1,m,q,iv) + mult(iv) * mat(2,nc,q,iv,i)
              end do
              do q = 1, nblk
                per2(2,m,q,iv) = per2(2,m,q,iv) + mult(iv) * mat(3,nc,q,iv,i)
                mat(2,m,q,iv,n) = mat(2,m,q,iv,n) + mult(iv) * per(nc,q,iv,i)
              end do
              rhs(m,iv,n) = rhs(m,iv,n) + mult(iv) * rhs(nc,iv,i)
            end do

          end do        ! end of loop on nc

          per2(1,:,:,iv) = per2(2,:,:,iv)

        end do          ! end of loop on i

!.... do the third to the last row

        i = n - 2

!.... do the first four columns

        do nc = 1, nblk - 1

          do m = nc, nblk - 1
            l = m + 1
            mult(iv) = -mat(2,l,nc,iv,i) / mat(2,nc,nc,iv,i)
            do q = nc + 1, nblk
              mat(2,l,q,iv,i) = mat(2,l,q,iv,i) + mult(iv) * mat(2,nc,q,iv,i)
            end do
            do q = 1, nblk
              mat(3,l,q,iv,i) = mat(3,l,q,iv,i) + mult(iv) * mat(3,nc,q,iv,i)
              per(l,q,iv,i) = per(l,q,iv,i) + mult(iv) * per(nc,q,iv,i)
            end do
            rhs(l,iv,i) = rhs(l,iv,i) + mult(iv) * rhs(nc,iv,i)
          end do

        end do  ! end of loop on nc

!.... for all columns

        do nc = 1, nblk

          p  = i + 1

          fact(iv) = one / mat(2,nc,nc,iv,i)
          do m = 1, nblk
            mult(iv) = -mat(1,m,nc,iv,p) * fact(iv)
            do q = nc + 1, nblk
              mat(1,m,q,iv,p) = mat(1,m,q,iv,p) + mult(iv) * mat(2,nc,q,iv,i)
            end do
            do q = 1, nblk
              mat(2,m,q,iv,p) = mat(2,m,q,iv,p) + mult(iv) * mat(3,nc,q,iv,i)
              mat(3,m,q,iv,p) = mat(3,m,q,iv,p) + mult(iv) * per(nc,q,iv,i)
            end do
            rhs(m,iv,p) = rhs(m,iv,p) + mult(iv) * rhs(nc,iv,i)

!.... the periodic stuff on the last row

            mult(iv) = -per2(1,m,nc,iv) * fact(iv)
            do q = nc + 1, nblk
              per2(1,m,q,iv) = per2(1,m,q,iv) + mult(iv) * mat(2,nc,q,iv,i)
            end do
            do q = 1, nblk
              mat(1,m,q,iv,n) = mat(1,m,q,iv,n) + mult(iv) * mat(3,nc,q,iv,i)
              mat(2,m,q,iv,n) = mat(2,m,q,iv,n) + mult(iv) * per(nc,q,iv,i)
            end do
            rhs(m,iv,n) = rhs(m,iv,n) + mult(iv) * rhs(nc,iv,i)
          end do

        end do  ! end of loop on nc

!.... do the second to the last row

        i = n - 1

!.... do the first four columns

          do nc = 1, nblk - 1

            do m = nc, nblk - 1
              l = m + 1
              mult(iv) = -mat(2,l,nc,iv,i) / mat(2,nc,nc,iv,i)
              do q = nc + 1, nblk
                mat(2,l,q,iv,i) = mat(2,l,q,iv,i) + mult(iv) * mat(2,nc,q,iv,i)
              end do
              do q = 1, nblk
                mat(3,l,q,iv,i) = mat(3,l,q,iv,i) + mult(iv) * mat(3,nc,q,iv,i)
              end do
              rhs(l,iv,i) = rhs(l,iv,i) + mult(iv) * rhs(nc,iv,i)
            end do

          end do        ! end of loop on nc

!.... for all columns

          do nc = 1, nblk

            p  = i + 1

            fact(iv) = one / mat(2,nc,nc,iv,i)
            do m = 1, nblk
              mult(iv) = -mat(1,m,nc,iv,p) * fact(iv)
              do q = nc + 1, nblk
                mat(1,m,q,iv,p) = mat(1,m,q,iv,p) + mult(iv) * mat(2,nc,q,iv,i)
              end do
              do q = 1, nblk
                mat(2,m,q,iv,p) = mat(2,m,q,iv,p) + mult(iv) * mat(3,nc,q,iv,i)
              end do
              rhs(m,iv,p) = rhs(m,iv,p) + mult(iv) * rhs(nc,iv,i)
            end do

          end do        ! end of loop on nc

!.... do the last row

        i = n

        do nc = 1, nblk - 1

          do m = nc, nblk - 1
            l = m + 1
            mult(iv) = -mat(2,l,nc,iv,i) / mat(2,nc,nc,iv,i)
            do q = nc + 1, nblk
              mat(2,l,q,iv,i) = mat(2,l,q,iv,i) + mult(iv) * mat(2,nc,q,iv,i)
            end do
            rhs(l,iv,i) = rhs(l,iv,i) + mult(iv) * rhs(nc,iv,i)
          end do

        end do  ! end of loop on nc

!=============================================================================!
!                   B A C K W A R D   S U B S T I T U T I O N
!=============================================================================!

!.... do the last row first

        i = n

        do m = nblk, 1, -1
          do q = m+1, nblk
            rhs(m,iv,i) = rhs(m,iv,i) - rhs(q,iv,i) * mat(2,m,q,iv,i)
          end do
          rhs(m,iv,i) = rhs(m,iv,i) / mat(2,m,m,iv,i)
        end do

!.... do the second to the last row

        i = n - 1

        p = i + 1

        do m = nblk, 1, -1
          do q = m+1, nblk
            rhs(m,iv,i) = rhs(m,iv,i) - rhs(q,iv,i) * mat(2,m,q,iv,i)
          end do
          do q= 1, nblk
            rhs(m,iv,i) = rhs(m,iv,i) - rhs(q,iv,p) * mat(3,m,q,iv,i)
          end do
          rhs(m,iv,i) = rhs(m,iv,i) / mat(2,m,m,iv,i)
        end do

!.... now do the rest of the rows in reverse order

        do i = n - 2, 1, -1

          p = i + 1

          do m = nblk, 1, -1
            do q = m+1, nblk
              rhs(m,iv,i) = rhs(m,iv,i) - rhs(q,iv,i) * mat(2,m,q,iv,i)
            end do
            do q = 1, nblk
              rhs(m,iv,i) = rhs(m,iv,i) - rhs(q,iv,p) * mat(3,m,q,iv,i)
              rhs(m,iv,i) = rhs(m,iv,i) - rhs(q,iv,n) * per(m,q,iv,i)
            end do
            rhs(m,iv,i) = rhs(m,iv,i) / mat(2,m,m,iv,i)
          end do

        end do

!.... account for periodicity

        rhs(1,iv,np) = rhs(1,iv,1)
        rhs(2,iv,np) = rhs(2,iv,1)
        rhs(3,iv,np) = rhs(3,iv,1)
        rhs(4,iv,np) = rhs(4,iv,1)
        rhs(5,iv,np) = rhs(5,iv,1)

        end do

        return
        end


text:
C $Id$
************************************************************************
*                                                                      *
      subroutine ecp_contract (n_ab,n_c,m_count,Qk,coeff_c,Q)
*                                                                      *
*   Contract ECP radial integrals over potential expansion             *
*                                                                      *
*   Argument (status) - description                                    *
*                                                                      *
*   n_ab (inp) - product of numbers of functions on centres a and b    *
*   n_c (inp) - number of functions on centre c (ecp centre)           *
*   m_count (inp) - number of m values                                 *
*   Qk (inp) - uncontracted ECP integrals                              *
*   coeff_c (inp) - coefficients of functions on centre c              *
*   Q (out) - contracted ECP integrals                                 *
*                                                                      *
*   Written by K. G. Dyall                                             *
*                                                                      *
************************************************************************
      implicit none
      integer i_ab,i_c,m,m_count,n_ab,n_c
      double precision Qk(n_ab,n_c,m_count),Q(n_ab,m_count),
     &    coeff_c(n_c)
*
      do m = 1,m_count
        do i_ab = 1,n_ab
          Q(i_ab,m) = 0.0d00
        end do
        do i_c = 1,n_c
C          write (6,'(1p4e20.10)') (Qk(i_ab,i_c,m),i_ab = 1,n_ab),
C     &        coeff_c(i_c)
          do i_ab = 1,n_ab
            Q(i_ab,m) = Q(i_ab,m)+Qk(i_ab,i_c,m)*coeff_c(i_c)
          end do
        end do
      end do
*
      return
      end

ext: f

max_stars_repo_name: dinisAbranches/nwchem

max_issues_repo_name: dinisAbranches/nwchem

max_forks_repo_name: dinisAbranches/nwchem


text:
!*******************************************************************************
!>  @file extcurz.f
!>  Contains module @ref extcurz.
!
!  Note separating the Doxygen comment block here so detailed decription is
!  found in the Module not the file.
!
!>  Represents a signal obtained by integrating around the magnetic field to get
!>  the current enclosed in the Z direction using Ampere's Law. Defines the base
!>  class of the type @ref extcurz_class.
!>  @par Super Class:
!>  @ref diagnostic
!*******************************************************************************

      MODULE extcurz

      USE stel_kinds, only: rprec
      USE mpi_inc
      USE profiler
      USE signal

      IMPLICIT NONE

!*******************************************************************************
!  DERIVED-TYPE DECLARATIONS
!  1) extcurz base class
!
!*******************************************************************************
!-------------------------------------------------------------------------------
!>  Base class representing a external Z currents enclosed by an inegration
!>  loop.
!>  @par Super Class:
!>  @ref diagnostic
!-------------------------------------------------------------------------------
      TYPE, EXTENDS(signal_class) :: extcurz_class
!>  Radial surface to compute external current from.
         REAL (rprec) :: r
!>  Poloidal angle to compute current from.
         REAL (rprec) :: theta
      CONTAINS
         PROCEDURE    ::                                                       &
     &      get_modeled_signal_last => extcurz_get_modeled_signal
         PROCEDURE    :: get_type => extcurz_get_type
         FINAL        :: extcurz_destruct
      END TYPE extcurz_class

!*******************************************************************************
!  INTERFACE BLOCKS
!*******************************************************************************
!-------------------------------------------------------------------------------
!>  Interface for the construction of @ref extcurz_class types using
!>  @ref extcurz_construct.
!-------------------------------------------------------------------------------
      INTERFACE extcurz_class
         MODULE PROCEDURE extcurz_construct
      END INTERFACE

      CONTAINS
!*******************************************************************************
!  CONSTRUCTION SUBROUTINES
!*******************************************************************************
!-------------------------------------------------------------------------------
!>  @brief Construct a @ref extcurz_class object measureing temperature.
!>
!>  Allocates memory and initializes a @ref extcurz_class object.
!>
!>  @param[in] r     The radial position to integrate at.
!>  @param[in] theta The poloidal angle to integrate at.
!>  @returns A pointer to a constructed @ref extcurz_class object.
!-------------------------------------------------------------------------------
      FUNCTION extcurz_construct(r, theta)

      IMPLICIT NONE

!  Declare Arguments
      CLASS (extcurz_class), POINTER :: extcurz_construct
      REAL(rprec), INTENT(in)        :: r
      REAL(rprec), INTENT(in)        :: theta

!  local variables
      REAL (rprec)                   :: start_time

!  Start of executable code
      start_time = profiler_get_start_time()

      ALLOCATE(extcurz_construct)
      extcurz_construct%r = r
      extcurz_construct%theta = theta

      CALL profiler_set_stop_time('extcurz_construct', start_time)

      END FUNCTION

!*******************************************************************************
!  DESTRUCTION SUBROUTINES
!*******************************************************************************
!-------------------------------------------------------------------------------
!>  @brief Deconstruct a @ref extcurz_class object.
!>
!>  Deallocates memory and uninitializes a @ref extcurz_class object.
!>
!>  @param[inout] this A @ref extcurz_class instance.
!-------------------------------------------------------------------------------
      SUBROUTINE extcurz_destruct(this)

      IMPLICIT NONE

!  Declare Arguments
      TYPE (extcurz_class), INTENT(inout) :: this

!  Start of executable code
      this%r = 0.0
      this%theta = 0.0

      END SUBROUTINE

!*******************************************************************************
!  GETTER SUBROUTINES
!*******************************************************************************
!-------------------------------------------------------------------------------
!>  @brief Calculates the modeled signal.
!>
!>  Calculates the inclused Z current from Ampere's Law from an loop integration
!>  of the magnetic field. The loop integrated magnetic field is provided by
!>  @ref equilibrium::equilibrium_get_Int_B_dphi
!>
!>  @param[inout] this       A @ref extcurz_class instance.
!>  @param[in]    a_model    A @ref model instance.
!>  @param[out]   sigma      The modeled sigma.
!>  @param[in]    last_value Last good value in case the signal did not change.
!>  @returns The model value.
!-------------------------------------------------------------------------------
      FUNCTION extcurz_get_modeled_signal(this, a_model, sigma,                &
     &                                    last_value)
      USE model
      USE stel_constants, only: mu0

      IMPLICIT NONE

!  Declare Arguments
      REAL (rprec), DIMENSION(4) :: extcurz_get_modeled_signal
      CLASS (extcurz_class), INTENT(inout)    :: this
      TYPE (model_class), POINTER             :: a_model
      REAL (rprec), DIMENSION(4), INTENT(out) :: sigma
      REAL (rprec), DIMENSION(4), INTENT(in)  :: last_value

!  local variables
      REAL (rprec)                            :: start_time

!  Start of executable code
      start_time = profiler_get_start_time()

      sigma = 0

      IF (BTEST(a_model%state_flags, model_state_vmec_flag)   .or.             &
     &    BTEST(a_model%state_flags, model_state_siesta_flag) .or.             &
     &    BTEST(a_model%state_flags, model_state_signal_flag)) THEN
         extcurz_get_modeled_signal(1) =                                       &
     &      a_model%equilibrium%get_Int_B_dphi(this%r, this%theta)/mu0

         CALL this%scale_and_offset(a_model,                                   &
     &                              extcurz_get_modeled_signal(1))
      ELSE
         extcurz_get_modeled_signal = last_value
      END IF

      CALL profiler_set_stop_time('extcurz_get_modeled_signal',                &
     &                            start_time)

      END FUNCTION

!-------------------------------------------------------------------------------
!>  @brief Gets a discription of the extcurz type.
!>
!>  Returns a description of the extcurz type for use when writting output files.
!>
!>  @param[in] this A @ref extcurz_class instance.
!>  @returns A string describing the extcurz type.
!-------------------------------------------------------------------------------
      FUNCTION extcurz_get_type(this)
      USE data_parameters

      IMPLICIT NONE

!  Declare Arguments
      CHARACTER (len=data_name_length)  :: extcurz_get_type
      CLASS (extcurz_class), INTENT(in) :: this

!  local variables
      REAL (rprec)                     :: start_time

!  Start of executable code
      start_time = profiler_get_start_time()

      extcurz_get_type = 'extcur'

      CALL profiler_set_stop_time('extcurz_get_type', start_time)

      END FUNCTION

      END MODULE


text:
      PROGRAM FM357

C***********************************************************************00010357
C*****  FORTRAN 77                                                      00020357
C*****   FM357               XAMOD - (159)                              00030357
C*****                                                                  00040357
C***********************************************************************00050357
C*****  GENERAL PURPOSE                                         ANS REF 00060357
C*****    TEST INTRINSIC FUNCTIONS AMOD AND MOD - REMAINDERING,  15.3   00070357
C*****    WHICH IS DEFINED AS A1-(A1/A2)A2 WHERE (X) IS AN     (TABLE 5)00080357
C*****    INTEGER WHOSE MAGNITUDE IS LE ABS(X) AND WHOSE SIGN           00090357
C*****    IS THE SAME AS X.                                             00100357
C*****                                                                  00110357
CBB** ********************** BBCCOMNT **********************************00120357
C****                                                                   00130357
C****            1978 FORTRAN COMPILER VALIDATION SYSTEM                00140357
C****                          VERSION 2.1                              00150357
C****                                                                   00160357
C****                                                                   00170357
C****           SUGGESTIONS AND COMMENTS SHOULD BE FORWARDED TO         00180357
C****          NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY           00190357
C****               SOFTWARE STANDARDS VALIDATION GROUP                 00200357
C****                      BUILDING 225  RM A266                        00210357
C****                     GAITHERSBURG, MD  20899                       00220357
C****                                                                   00230357
C****                                                                   00240357
C****                                                                   00250357
CBE** ********************** BBCCOMNT **********************************00260357
CBB** ********************** BBCINITA **********************************00270357
C**** SPECIFICATION STATEMENTS                                          00280357
C****                                                                   00290357
      CHARACTER ZVERS*13, ZVERSD*17, ZDATE*17, ZPROG*5, ZCOMPL*20,      00300357
     1          ZNAME*20, ZTAPE*10, ZPROJ*13, REMRKS*31, ZTAPED*13      00310357
CBE** ********************** BBCINITA **********************************00320357
CBB** ********************** BBCINITB **********************************00330357
C**** INITIALIZE SECTION                                                00340357
      DATA  ZVERS,                  ZVERSD,             ZDATE           00350357
     1      /'VERSION 2.1  ',  '93/10/21*21.02.00',  '*NO DATE*TIME'/   00360357
      DATA       ZCOMPL,             ZNAME,             ZTAPE           00370357
     1      /'*NONE SPECIFIED*', '*NO COMPANY NAME*', '*NO TAPE*'/      00380357
      DATA       ZPROJ,           ZTAPED,         ZPROG                 00390357
     1      /'*NO PROJECT*',   '*NO TAPE DATE',  'XXXXX'/               00400357
      DATA   REMRKS /'                               '/                 00410357
C**** THE FOLLOWING 9 COMMENT LINES (CZ01, CZ02, ...) CAN BE REPLACED   00420357
C**** FOR IDENTIFYING THE TEST ENVIRONMENT                              00430357
C****                                                                   00440357
CZ01  ZVERS  = 'VERSION OF THE COMPILER VALIDATION SYSTEM'              00450357
CZ02  ZVERSD = 'CREATION DATE/TIME OF THE COMPILER VALIDATION SYSTEM'   00460357
CZ03  ZPROG  = 'PROGRAM NAME'                                           00470357
CZ04  ZDATE  = 'DATE OF TEST'                                           00480357
CZ05  ZCOMPL = 'COMPILER IDENTIFICATION'                                00490357
CZ06  ZPROJ  = 'PROJECT NUMBER/IDENTIFICATION'                          00500357
CZ07  ZNAME  = 'NAME OF USER'                                           00510357
CZ08  ZTAPE  = 'TAPE OWNER/ID'                                          00520357
CZ09  ZTAPED = 'DATE TAPE COPIED'                                       00530357
C                                                                       00540357
      IVPASS = 0                                                        00550357
      IVFAIL = 0                                                        00560357
      IVDELE = 0                                                        00570357
      IVINSP = 0                                                        00580357
      IVTOTL = 0                                                        00590357
      IVTOTN = 0                                                        00600357
      ICZERO = 0                                                        00610357
C                                                                       00620357
C     I01 CONTAINS THE LOGICAL UNIT NUMBER FOR THE CARD READER.         00630357
      I01 = 05                                                          00640357
C     I02 CONTAINS THE LOGICAL UNIT NUMBER FOR THE PRINTER.             00650357
      I02 = 06                                                          00660357
C                                                                       00670357
CX010   REPLACED BY FEXEC X-010 CONTROL CARD (CARD-READER UNIT NUMBER). 00680357
C     THE CX010 CARD IS FOR OVERRIDING THE PROGRAM DEFAULT I01 = 5      00690357
CX011   REPLACED BY FEXEC X-011 CONTROL CARD.  CX011 IS FOR SYSTEMS     00700357
C     REQUIRING ADDITIONAL STATEMENTS FOR FILES ASSOCIATED WITH CX010.  00710357
C                                                                       00720357
CX020   REPLACED BY FEXEC X-020 CONTROL CARD (PRINTER UNIT NUMBER).     00730357
C     THE CX020 CARD IS FOR OVERRIDING THE PROGRAM DEFAULT I02= 6       00740357
CX021   REPLACED BY FEXEC X-021 CONTROL CARD.  CX021 IS FOR SYSTEMS     00750357
C     REQUIRING ADDITIONAL STATEMENTS FOR FILES ASSOCIATED WITH CX020.  00760357
C                                                                       00770357
CBE** ********************** BBCINITB **********************************00780357
      NUVI = I02                                                        00790357
      IVTOTL = 22                                                       00800357
      ZPROG = 'FM357'                                                   00810357
CBB** ********************** BBCHED0A **********************************00820357
C****                                                                   00830357
C**** WRITE REPORT TITLE                                                00840357
C****                                                                   00850357
      WRITE (I02, 90002)                                                00860357
      WRITE (I02, 90006)                                                00870357
      WRITE (I02, 90007)                                                00880357
      WRITE (I02, 90008)  ZVERS, ZVERSD                                 00890357
      WRITE (I02, 90009)  ZPROG, ZPROG                                  00900357
      WRITE (I02, 90010)  ZDATE, ZCOMPL                                 00910357
CBE** ********************** BBCHED0A **********************************00920357
C*****                                                                  00930357
C*****    HEADER FOR SEGMENT 159 WRITTEN                                00940357
        WRITE (NUVI,15901)                                              00950357
15901   FORMAT (" ", //, 2X,"XAMOD - (159) INTRINSIC FUNCTION-- " //16X,00960357
     1          "AMOD, MOD (REMAINDERING)" //" SUBSET REF. - 15.3" )    00970357
CBB** ********************** BBCHED0B **********************************00980357
C**** WRITE DETAIL REPORT HEADERS                                       00990357
C****                                                                   01000357
      WRITE (I02,90004)                                                 01010357
      WRITE (I02,90004)                                                 01020357
      WRITE (I02,90013)                                                 01030357
      WRITE (I02,90014)                                                 01040357
      WRITE (I02,90015) IVTOTL                                          01050357
CBE** ********************** BBCHED0B **********************************01060357
C*****                                                                  01070357
C*****    TEST OF AMOD                                                  01080357
C*****                                                                  01090357
        WRITE(NUVI, 15902)                                              01100357
15902   FORMAT (/ 8X, "TEST OF AMOD" )                                  01110357
C*****                                                                  01120357
CT001*  TEST 1                      FIRST VALUE ZERO, SECOND NON-ZERO   01130357
           IVTNUM = 1                                                   01140357
        REBVS = 0.0                                                     01150357
        REDVS = 4.5                                                     01160357
        REAVS = AMOD(REBVS, REDVS)                                      01170357
           IF (REAVS + 0.00005) 20010, 10010, 40010                     01180357
40010      IF (REAVS - 0.00005) 10010, 10010, 20010                     01190357
10010      IVPASS = IVPASS + 1                                          01200357
           WRITE (NUVI, 80002) IVTNUM                                   01210357
           GO TO 0011                                                   01220357
20010      IVFAIL = IVFAIL + 1                                          01230357
           RVCORR = 0.0                                                 01240357
           WRITE (NUVI,80012) IVTNUM, REAVS, RVCORR                     01250357
 0011      CONTINUE                                                     01260357
CT002*  TEST 2                                      BOTH VALUES EQUAL   01270357
           IVTNUM = 2                                                   01280357
        REBVS = 3.5                                                     01290357
        REDVS = 3.5                                                     01300357
        REAVS = AMOD(REBVS, REDVS)                                      01310357
           IF (REAVS + 0.00005) 20020, 10020, 40020                     01320357
40020      IF (REAVS - 0.00005) 10020, 10020, 20020                     01330357
10020      IVPASS = IVPASS + 1                                          01340357
           WRITE (NUVI, 80002) IVTNUM                                   01350357
           GO TO 0021                                                   01360357
20020      IVFAIL = IVFAIL + 1                                          01370357
           RVCORR = 0.0                                                 01380357
           WRITE(NUVI, 80012) IVTNUM, REAVS, RVCORR                     01390357
 0021      CONTINUE                                                     01400357
CT003*  TEST 3         FIRST MAGNITUDE LARGER, NOT MULTIPLE OF SECOND   01410357
           IVTNUM = 3                                                   01420357
        REBVS = -10.9                                                   01430357
        REDVS = -3.3                                                    01440357
        REAVS = AMOD(REBVS, REDVS)                                      01450357
           IF (REAVS + 1.0001) 20030, 10030, 40030                      01460357
40030      IF (REAVS + 0.99995) 10030, 10030, 20030                     01470357
10030      IVPASS = IVPASS + 1                                          01480357
           WRITE (NUVI, 80002) IVTNUM                                   01490357
           GO TO 0031                                                   01500357
20030      IVFAIL = IVFAIL + 1                                          01510357
           RVCORR = -1.0                                                01520357
           WRITE(NUVI, 80012) IVTNUM, REAVS, RVCORR                     01530357
 0031      CONTINUE                                                     01540357
CT004*  TEST 4             FIRST MAGNITUDE LARGER, MULTIPLE OF SECOND   01550357
           IVTNUM = 4                                                   01560357
        REDVS = 1.5                                                     01570357
        REBVS = 1.5 + REDVS + 1.5                                       01580357
        REAVS = AMOD(REBVS, REDVS)                                      01590357
           IF (REAVS + 0.00005) 20040, 10040, 40040                     01600357
40040      IF (REAVS - 0.00005) 10040, 10040, 20040                     01610357
10040      IVPASS = IVPASS + 1                                          01620357
           WRITE (NUVI, 80002) IVTNUM                                   01630357
           GO TO 0041                                                   01640357
20040      IVFAIL = IVFAIL + 1                                          01650357
           RVCORR = 0.0                                                 01660357
           WRITE(NUVI, 80012) IVTNUM, REAVS, RVCORR                     01670357
 0041      CONTINUE                                                     01680357
CT005*  TEST 5         FIRST MAGNITUDE LARGER, NOT MULTIPLE OF SECOND   01690357
           IVTNUM = 5                                                   01700357
        REBVS = 7.625                                                   01710357
        REDVS = 2.125                                                   01720357
        REAVS = AMOD(REBVS, REDVS)                                      01730357
           IF (REAVS - 1.2499) 20050, 10050, 40050                      01740357
40050      IF (REAVS - 1.2501) 10050, 10050, 20050                      01750357
10050      IVPASS = IVPASS + 1                                          01760357
           WRITE (NUVI, 80002) IVTNUM                                   01770357
           GO TO 0051                                                   01780357
20050      IVFAIL = IVFAIL + 1                                          01790357
           RVCORR = 1.25                                                01800357
           WRITE(NUVI, 80012) IVTNUM, REAVS, RVCORR                     01810357
 0051      CONTINUE                                                     01820357
CT006*  TEST 6                      FIRST VALUE ZERO, SECOND NEGATIVE   01830357
           IVTNUM = 6                                                   01840357
        REBVS = 0.0                                                     01850357
        REDVS = -4.5                                                    01860357
        REAVS = AMOD(REBVS, REDVS)                                      01870357
           IF (REAVS + 0.00005) 20060, 10060, 40060                     01880357
40060      IF (REAVS - 0.00005) 10060, 10060, 20060                     01890357
10060      IVPASS = IVPASS + 1                                          01900357
           WRITE (NUVI, 80002) IVTNUM                                   01910357
           GO TO 0061                                                   01920357
20060      IVFAIL = IVFAIL + 1                                          01930357
           RVCORR = 0.0                                                 01940357
           WRITE(NUVI, 80012) IVTNUM, REAVS, RVCORR                     01950357
 0061      CONTINUE                                                     01960357
CT007*  TEST 7                       BOTH VALUES EQUAL, BOTH NEGATIVE   01970357
           IVTNUM = 7                                                   01980357
        REBVS = -3.5                                                    01990357
        REDVS = -3.5                                                    02000357
        REAVS = AMOD(REBVS, REDVS)                                      02010357
           IF (REAVS + 0.00005) 20070, 10070, 40070                     02020357
40070      IF (REAVS - 0.00005) 10070, 10070, 20070                     02030357
10070      IVPASS = IVPASS + 1                                          02040357
           WRITE (NUVI, 80002) IVTNUM                                   02050357
           GO TO 0071                                                   02060357
20070      IVFAIL = IVFAIL + 1                                          02070357
           RVCORR = 0.0                                                 02080357
           WRITE(NUVI, 80012) IVTNUM, REAVS, RVCORR                     02090357
 0071      CONTINUE                                                     02100357
CT008*  TEST 8        FIRST VALUE NEGATIVE, SECOND POSITIVE, MULTIPLE   02110357
           IVTNUM = 8                                                   02120357
        REBVS = 1.5                                                     02130357
        REDVS = -(1.5 + REDVS + 1.5)                                    02140357
        REAVS = AMOD(-REBVS, -REDVS)                                    02150357
           IF (REAVS + 0.00005) 20080, 10080, 40080                     02160357
40080      IF (REAVS - 0.00005) 10080, 10080, 20080                     02170357
10080      IVPASS = IVPASS + 1                                          02180357
           WRITE (NUVI, 80002) IVTNUM                                   02190357
           GO TO 0081                                                   02200357
20080      IVFAIL = IVFAIL + 1                                          02210357
           RVCORR = 0.0                                                 02220357
           WRITE(NUVI, 80012) IVTNUM, REAVS, RVCORR                     02230357
 0081      CONTINUE                                                     02240357
CT009*  TEST 9         FIRST MAGNITUDE LARGER, NOT MULTIPLE OF SECOND   02250357
           IVTNUM = 9                                                   02260357
        REBVS = 10.5                                                    02270357
        REDVS = -3.3                                                    02280357
        REAVS = AMOD(REBVS, REDVS)                                      02290357
           IF (REAVS - 0.59997) 20090, 10090, 40090                     02300357
40090      IF (REAVS - 0.60003) 10090, 10090, 20090                     02310357
10090      IVPASS = IVPASS + 1                                          02320357
           WRITE (NUVI, 80002) IVTNUM                                   02330357
           GO TO 0091                                                   02340357
20090      IVFAIL = IVFAIL + 1                                          02350357
           RVCORR = 0.6                                                 02360357
           WRITE(NUVI, 80012) IVTNUM, REAVS, RVCORR                     02370357
 0091      CONTINUE                                                     02380357
CT010*  TEST 10       PAIR OF ARITHMETIC EXPRESSIONS USED AS ARGUMENT   02390357
           IVTNUM = 10                                                  02400357
        RECVS = 7.625                                                   02410357
        REDVS = 2.125                                                   02420357
        REFVS = 2.0                                                     02430357
        REAVS = AMOD(RECVS - REFVS, REDVS + REFVS)                      02440357
           IF (REAVS - 1.4999) 20100, 10100, 40100                      02450357
40100      IF (REAVS - 1.5001) 10100, 10100, 20100                      02460357
10100      IVPASS = IVPASS + 1                                          02470357
           WRITE (NUVI, 80002) IVTNUM                                   02480357
           GO TO 0101                                                   02490357
20100      IVFAIL = IVFAIL + 1                                          02500357
           RVCORR = 1.5                                                 02510357
           WRITE(NUVI, 80012) IVTNUM, REAVS, RVCORR                     02520357
 0101      CONTINUE                                                     02530357
CT011*  TEST 11                TEST LOW AND HIGH MAGNITUDE ARGUMENTS    02540357
           IVTNUM = 11                                                  02550357
        RECVS = 1.0E-16                                                 02560357
        REDVS = 1.0E+16                                                 02570357
        REAVS = AMOD(RECVS, REDVS)                                      02580357
           IF (REAVS - 0.99995E-16) 20110, 10110, 40110                 02590357
40110      IF (REAVS - 1.0001E-16) 10110, 10110, 20110                  02600357
10110      IVPASS = IVPASS + 1                                          02610357
           WRITE (NUVI, 80002) IVTNUM                                   02620357
           GO TO 0111                                                   02630357
20110      IVFAIL = IVFAIL + 1                                          02640357
           RVCORR = 1.0E-16                                             02650357
           WRITE(NUVI, 80012) IVTNUM, REAVS, RVCORR                     02660357
 0111      CONTINUE                                                     02670357
C*****                                                                  02680357
C*****    TEST OF MOD                                                   02690357
C*****                                                                  02700357
        WRITE(NUVI, 15904)                                              02710357
15904   FORMAT (/ 8X, "TEST OF MOD" )                                   02720357
C*****                                                                  02730357
CT012*  TEST 12                     FIRST VALUE ZERO, SECOND NON-ZERO   02740357
           IVTNUM = 12                                                  02750357
        IEBVI = 0                                                       02760357
        IEDVI = 4                                                       02770357
        IEAVI = MOD(IEBVI, IEDVI)                                       02780357
           IF (IEAVI - 0) 20120, 10120, 20120                           02790357
10120      IVPASS = IVPASS + 1                                          02800357
           WRITE (NUVI, 80002) IVTNUM                                   02810357
           GO TO 0121                                                   02820357
20120      IVFAIL = IVFAIL + 1                                          02830357
           IVCORR = 0                                                   02840357
           WRITE (NUVI, 80010) IVTNUM, IEAVI, IVCORR                    02850357
 0121      CONTINUE                                                     02860357
CT013*  TEST 13                                     BOTH VALUES EQUAL   02870357
           IVTNUM = 13                                                  02880357
        IEBVI = 3                                                       02890357
        IEDVI = 3                                                       02900357
        IEAVI = MOD(IEBVI, IEDVI)                                       02910357
           IF (IEAVI - 0) 20130, 10130, 20130                           02920357
10130      IVPASS = IVPASS + 1                                          02930357
           WRITE (NUVI, 80002) IVTNUM                                   02940357
           GO TO 0131                                                   02950357
20130      IVFAIL = IVFAIL + 1                                          02960357
           IVCORR = 0                                                   02970357
           WRITE (NUVI, 80010) IVTNUM, IEAVI, IVCORR                    02980357
 0131      CONTINUE                                                     02990357
CT014*  TEST 14        FIRST MAGNITUDE LARGER, NOT MULTIPLE OF SECOND   03000357
           IVTNUM = 14                                                  03010357
        IEBVI = -10                                                     03020357
        IEDVI = -3                                                      03030357
        IEAVI = MOD(IEBVI, IEDVI)                                       03040357
           IF (IEAVI + 1) 20140, 10140, 20140                           03050357
10140      IVPASS = IVPASS + 1                                          03060357
           WRITE (NUVI, 80002) IVTNUM                                   03070357
           GO TO 0141                                                   03080357
20140      IVFAIL = IVFAIL + 1                                          03090357
           IVCORR = -1                                                  03100357
           WRITE (NUVI, 80010) IVTNUM, IEAVI, IVCORR                    03110357
 0141      CONTINUE                                                     03120357
CT015*  TEST 15            FIRST MAGNITUDE LARGER, MULTIPLE OF SECOND   03130357
           IVTNUM = 15                                                  03140357
        IEBVI = 9                                                       03150357
        IEDVI = 3                                                       03160357
        IEAVI = MOD(IEBVI, IEDVI)                                       03170357
           IF (IEAVI - 0) 20150, 10150, 20150                           03180357
10150      IVPASS = IVPASS + 1                                          03190357
           WRITE (NUVI, 80002) IVTNUM                                   03200357
           GO TO 0151                                                   03210357
20150      IVFAIL = IVFAIL + 1                                          03220357
           IVCORR = 0                                                   03230357
           WRITE (NUVI, 80010) IVTNUM, IEAVI, IVCORR                    03240357
 0151      CONTINUE                                                     03250357
CT016*  TEST 16        FIRST MAGNITUDE LARGER, NOT MULTIPLE OF SECOND   03260357
           IVTNUM = 16                                                  03270357
        IEBVI = 7                                                       03280357
        IEDVI = 2                                                       03290357
        IEAVI = MOD(IEBVI, IEDVI)                                       03300357
           IF (IEAVI - 1) 20160, 10160, 20160                           03310357
10160      IVPASS = IVPASS + 1                                          03320357
           WRITE (NUVI, 80002) IVTNUM                                   03330357
           GO TO 0161                                                   03340357
20160      IVFAIL = IVFAIL + 1                                          03350357
           IVCORR = 1                                                   03360357
           WRITE (NUVI, 80010) IVTNUM, IEAVI, IVCORR                    03370357
 0161      CONTINUE                                                     03380357
CT017*  TEST 17                     FIRST VALUE ZERO, SECOND NEGATIVE   03390357
           IVTNUM = 17                                                  03400357
        IEBVI = 0                                                       03410357
        IEDVI = -4                                                      03420357
        IEAVI = MOD(IEBVI, IEDVI)                                       03430357
           IF (IEAVI - 0) 20170, 10170, 20170                           03440357
10170      IVPASS = IVPASS + 1                                          03450357
           WRITE (NUVI, 80002) IVTNUM                                   03460357
           GO TO 0171                                                   03470357
20170      IVFAIL = IVFAIL + 1                                          03480357
           IVCORR = 0                                                   03490357
           WRITE (NUVI, 80010) IVTNUM, IEAVI, IVCORR                    03500357
 0171      CONTINUE                                                     03510357
CT018*  TEST 18                      BOTH VALUES EQUAL, BOTH NEGATIVE   03520357
           IVTNUM = 18                                                  03530357
        IEBVI = -3                                                      03540357
        IEDVI = -3                                                      03550357
        IEAVI = MOD(IEBVI, IEDVI)                                       03560357
           IF (IEAVI - 0) 20180, 10180, 20180                           03570357
10180      IVPASS = IVPASS + 1                                          03580357
           WRITE (NUVI, 80002) IVTNUM                                   03590357
           GO TO 0181                                                   03600357
20180      IVFAIL = IVFAIL + 1                                          03610357
           IVCORR = 0                                                   03620357
           WRITE (NUVI, 80010) IVTNUM, IEAVI, IVCORR                    03630357
 0181      CONTINUE                                                     03640357
CT019*  TEST 19       FIRST MAGNITUDE LARGER, MULTIPLE, BOTH NEGATIVE   03650357
           IVTNUM = 19                                                  03660357
        IEBVI = -9                                                      03670357
        IEDVI = -3                                                      03680357
        IEAVI = MOD(IEBVI, IEDVI)                                       03690357
           IF (IEAVI - 0) 20190, 10190, 20190                           03700357
10190      IVPASS = IVPASS + 1                                          03710357
           WRITE (NUVI, 80002) IVTNUM                                   03720357
           GO TO 0191                                                   03730357
20190      IVFAIL = IVFAIL + 1                                          03740357
           IVCORR = 0                                                   03750357
           WRITE (NUVI, 80010) IVTNUM, IEAVI, IVCORR                    03760357
 0191      CONTINUE                                                     03770357
CT020*  TEST 20      FIRST NUMBER NEGATIVE, SECOND POSITIVE, MULTIPLE   03780357
           IVTNUM = 20                                                  03790357
        IEBVI = -9                                                      03800357
        IEDVI = 3                                                       03810357
        IEAVI = MOD(IEBVI, IEDVI)                                       03820357
           IF (IEAVI - 0) 20200, 10200, 20200                           03830357
10200      IVPASS = IVPASS + 1                                          03840357
           WRITE (NUVI, 80002) IVTNUM                                   03850357
           GO TO 0201                                                   03860357
20200      IVFAIL = IVFAIL + 1                                          03870357
           IVCORR = 0                                                   03880357
           WRITE (NUVI, 80010) IVTNUM, IEAVI, IVCORR                    03890357
 0201      CONTINUE                                                     03900357
CT021*  TEST 21               FIRST VALUE ZERO PRECEDED BY MINUS SIGN   03910357
           IVTNUM = 21                                                  03920357
        IEBVI = 0                                                       03930357
        IEDVI = 4                                                       03940357
        IEAVI = MOD(-IEBVI, IEDVI)                                      03950357
           IF (IEAVI - 0) 20210, 10210, 20210                           03960357
10210      IVPASS = IVPASS + 1                                          03970357
           WRITE (NUVI, 80002) IVTNUM                                   03980357
           GO TO 0211                                                   03990357
20210      IVFAIL = IVFAIL + 1                                          04000357
           IVCORR = 0                                                   04010357
           WRITE (NUVI, 80010) IVTNUM, IEAVI, IVCORR                    04020357
 0211      CONTINUE                                                     04030357
CT022*  TEST 22       PAIR OF ARITHMETIC EXPRESSIONS USED AS ARGUMENT   04040357
           IVTNUM = 22                                                  04050357
        IEDVI = 10                                                      04060357
        IEEVI = 3                                                       04070357
        IEFVI = 2                                                       04080357
        IEAVI = MOD(IEDVI - IEFVI, IEEVI + IEFVI)                       04090357
           IF (IEAVI - 3) 20220, 10220, 20220                           04100357
10220      IVPASS = IVPASS + 1                                          04110357
           WRITE (NUVI, 80002) IVTNUM                                   04120357
           GO TO 0221                                                   04130357
20220      IVFAIL = IVFAIL + 1                                          04140357
           IVCORR = 3                                                   04150357
           WRITE (NUVI, 80010) IVTNUM, IEAVI, IVCORR                    04160357
 0221      CONTINUE                                                     04170357
C*****                                                                  04180357
CBB** ********************** BBCSUM0  **********************************04190357
C**** WRITE OUT TEST SUMMARY                                            04200357
C****                                                                   04210357
      IVTOTN = IVPASS + IVFAIL + IVDELE + IVINSP                        04220357
      WRITE (I02, 90004)                                                04230357
      WRITE (I02, 90014)                                                04240357
      WRITE (I02, 90004)                                                04250357
      WRITE (I02, 90020) IVPASS                                         04260357
      WRITE (I02, 90022) IVFAIL                                         04270357
      WRITE (I02, 90024) IVDELE                                         04280357
      WRITE (I02, 90026) IVINSP                                         04290357
      WRITE (I02, 90028) IVTOTN, IVTOTL                                 04300357
CBE** ********************** BBCSUM0  **********************************04310357
CBB** ********************** BBCFOOT0 **********************************04320357
C**** WRITE OUT REPORT FOOTINGS                                         04330357
C****                                                                   04340357
      WRITE (I02,90016) ZPROG, ZPROG                                    04350357
      WRITE (I02,90018) ZPROJ, ZNAME, ZTAPE, ZTAPED                     04360357
      WRITE (I02,90019)                                                 04370357
CBE** ********************** BBCFOOT0 **********************************04380357
CBB** ********************** BBCFMT0A **********************************04390357
C**** FORMATS FOR TEST DETAIL LINES                                     04400357
C****                                                                   04410357
80000 FORMAT (" ",2X,I3,4X,"DELETED",32X,A31)                           04420357
80002 FORMAT (" ",2X,I3,4X," PASS  ",32X,A31)                           04430357
80004 FORMAT (" ",2X,I3,4X,"INSPECT",32X,A31)                           04440357
80008 FORMAT (" ",2X,I3,4X," FAIL  ",32X,A31)                           04450357
80010 FORMAT (" ",2X,I3,4X," FAIL  ",/," ",15X,"COMPUTED= " ,           04460357
     1I6,/," ",15X,"CORRECT=  " ,I6)                                    04470357
80012 FORMAT (" ",2X,I3,4X," FAIL  ",/," ",16X,"COMPUTED= " ,           04480357
     1E12.5,/," ",16X,"CORRECT=  " ,E12.5)                              04490357
80018 FORMAT (" ",2X,I3,4X," FAIL  ",/," ",16X,"COMPUTED= " ,           04500357
     1A21,/," ",16X,"CORRECT=  " ,A21)                                  04510357
80020 FORMAT (" ",16X,"COMPUTED= " ,A21,1X,A31)                         04520357
80022 FORMAT (" ",16X,"CORRECT=  " ,A21,1X,A31)                         04530357
80024 FORMAT (" ",16X,"COMPUTED= " ,I6,16X,A31)                         04540357
80026 FORMAT (" ",16X,"CORRECT=  " ,I6,16X,A31)                         04550357
80028 FORMAT (" ",16X,"COMPUTED= " ,E12.5,10X,A31)                      04560357
80030 FORMAT (" ",16X,"CORRECT=  " ,E12.5,10X,A31)                      04570357
80050 FORMAT (" ",48X,A31)                                              04580357
CBE** ********************** BBCFMT0A **********************************04590357
CBB** ********************** BBCFMT0B **********************************04600357
C**** FORMAT STATEMENTS FOR PAGE HEADERS                                04610357
C****                                                                   04620357
90002 FORMAT ("1")                                                      04630357
90004 FORMAT (" ")                                                      04640357
90006 FORMAT (" ",20X,"NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY" )04650357
90007 FORMAT (" ",19X,"FORTRAN COMPILER VALIDATION SYSTEM" )            04660357
90008 FORMAT (" ",21X,A13,A17)                                          04670357
90009 FORMAT (" ",/," *",A5,"BEGIN*",12X,"TEST RESULTS - " ,A5,/)       04680357
90010 FORMAT (" ",8X,"TEST DATE*TIME= " ,A17,"  -  COMPILER= " ,A20)    04690357
90013 FORMAT (" "," TEST   ","PASS/FAIL " ,6X,"DISPLAYED RESULTS" ,     04700357
     1       7X,"REMARKS",24X)                                          04710357
90014 FORMAT (" ","----------------------------------------------" ,    04720357
     1        "---------------------------------" )                     04730357
90015 FORMAT (" ",48X,"THIS PROGRAM HAS " ,I3," TESTS",/)               04740357
C****                                                                   04750357
C**** FORMAT STATEMENTS FOR REPORT FOOTINGS                             04760357
C****                                                                   04770357
90016 FORMAT (" ",/," *",A5,"END*",14X,"END OF TEST - " ,A5,/)          04780357
90018 FORMAT (" ",A13,13X,A20,"   *   ",A10,"/",                        04790357
     1        A13)                                                      04800357
90019 FORMAT (" ","FOR OFFICIAL USE ONLY     " ,35X,"COPYRIGHT  1982" ) 04810357
C****                                                                   04820357
C**** FORMAT STATEMENTS FOR RUN SUMMARY                                 04830357
C****                                                                   04840357
90020 FORMAT (" ",21X,I5," TESTS PASSED" )                              04850357
90022 FORMAT (" ",21X,I5," TESTS FAILED" )                              04860357
90024 FORMAT (" ",21X,I5," TESTS DELETED" )                             04870357
90026 FORMAT (" ",21X,I5," TESTS REQUIRE INSPECTION" )                  04880357
90028 FORMAT (" ",21X,I5," OF ",I3," TESTS EXECUTED" )                  04890357
CBE** ********************** BBCFMT0B **********************************04900357
C*****                                                                  04910357
C*****    END OF TEST SEGMENT 159                                       04920357
      STOP                                                              04930357
      END                                                               04940357
                                                                        04950357

max_stars_repo_name: MrSidims/llvm-test-suite

max_issues_repo_name: MrSidims/llvm-test-suite

max_forks_repo_name: MrSidims/llvm-test-suite


text:
*DECK ORTBAK
      SUBROUTINE ORTBAK (NM, LOW, IGH, A, ORT, M, Z)
C***BEGIN PROLOGUE  ORTBAK
C***PURPOSE  Form the eigenvectors of a general real matrix from the
C            eigenvectors of the upper Hessenberg matrix output from
C            ORTHES.
C***LIBRARY   SLATEC (EISPACK)
C***CATEGORY  D4C4
C***TYPE      SINGLE PRECISION (ORTBAK-S, CORTB-C)
C***KEYWORDS  EIGENVALUES, EIGENVECTORS, EISPACK
C***AUTHOR  Smith, B. T., et al.
C***DESCRIPTION
C
C     This subroutine is a translation of the ALGOL procedure ORTBAK,
C     NUM. MATH. 12, 349-368(1968) by Martin and Wilkinson.
C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 339-358(1971).
C
C     This subroutine forms the eigenvectors of a REAL GENERAL
C     matrix by back transforming those of the corresponding
C     upper Hessenberg matrix determined by  ORTHES.
C
C     On INPUT
C
C        NM must be set to the row dimension of the two-dimensional
C          array parameters, A and Z, as declared in the calling
C          program dimension statement.  NM is an INTEGER variable.
C
C        LOW and IGH are two INTEGER variables determined by the
C          balancing subroutine  BALANC.  If  BALANC  has not been
C          used, set LOW=1 and IGH equal to the order of the matrix.
C
C        A contains some information about the orthogonal trans-
C          formations used in the reduction to Hessenberg form by
C          ORTHES  in its strict lower triangle.  A is a two-dimensional
C          REAL array, dimensioned A(NM,IGH).
C
C        ORT contains further information about the orthogonal trans-
C          formations used in the reduction by  ORTHES.  Only elements
C          LOW through IGH are used.  ORT is a one-dimensional REAL
C          array, dimensioned ORT(IGH).
C
C        M is the number of columns of Z to be back transformed.
C          M is an INTEGER variable.
C
C        Z contains the real and imaginary parts of the eigenvectors to
C          be back transformed in its first M columns.  Z is a two-
C          dimensional REAL array, dimensioned Z(NM,M).
C
C     On OUTPUT
C
C        Z contains the real and imaginary parts of the transformed
C          eigenvectors in its first M columns.
C
C        ORT has been used for temporary storage as is not restored.
C
C     NOTE that ORTBAK preserves vector Euclidean norms.
C
C     Questions and comments should be directed to B. S. Garbow,
C     APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
C     ------------------------------------------------------------------
C
C***REFERENCES  B. T. Smith, J. M. Boyle, J. J. Dongarra, B. S. Garbow,
C                 Y. Ikebe, V. C. Klema and C. B. Moler, Matrix Eigen-
C                 system Routines - EISPACK Guide, Springer-Verlag,
C                 1976.
C***ROUTINES CALLED  (NONE)
C***REVISION HISTORY  (YYMMDD)
C   760101  DATE WRITTEN
C   890831  Modified array declarations.  (WRB)
C   890831  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   920501  Reformatted the REFERENCES section.  (WRB)
C***END PROLOGUE  ORTBAK
C
      INTEGER I,J,M,LA,MM,MP,NM,IGH,KP1,LOW,MP1
      REAL A(NM,*),ORT(*),Z(NM,*)
      REAL G
C
C***FIRST EXECUTABLE STATEMENT  ORTBAK
      IF (M .EQ. 0) GO TO 200
      LA = IGH - 1
      KP1 = LOW + 1
      IF (LA .LT. KP1) GO TO 200
C     .......... FOR MP=IGH-1 STEP -1 UNTIL LOW+1 DO -- ..........
      DO 140 MM = KP1, LA
         MP = LOW + IGH - MM
         IF (A(MP,MP-1) .EQ. 0.0E0) GO TO 140
         MP1 = MP + 1
C
         DO 100 I = MP1, IGH
  100    ORT(I) = A(I,MP-1)
C
         DO 130 J = 1, M
            G = 0.0E0
C
            DO 110 I = MP, IGH
  110       G = G + ORT(I) * Z(I,J)
C     .......... DIVISOR BELOW IS NEGATIVE OF H FORMED IN ORTHES.
C                DOUBLE DIVISION AVOIDS POSSIBLE UNDERFLOW ..........
            G = (G / ORT(MP)) / A(MP,MP-1)
C
            DO 120 I = MP, IGH
  120       Z(I,J) = Z(I,J) + G * ORT(I)
C
  130    CONTINUE
C
  140 CONTINUE
C
  200 RETURN
      END


text:
C Copyright(C) 2011-2017 National Technology & Engineering Solutions of
C Sandia, LLC (NTESS).  Under the terms of Contract DE-NA0003525 with
C NTESS, the U.S. Government retains certain rights in this software.
C
C Redistribution and use in source and binary forms, with or without
C modification, are permitted provided that the following conditions are
C met:
C
C * Redistributions of source code must retain the above copyright
C    notice, this list of conditions and the following disclaimer.
C
C * Redistributions in binary form must reproduce the above
C   copyright notice, this list of conditions and the following
C   disclaimer in the documentation and/or other materials provided
C   with the distribution.
C
C * Neither the name of NTESS nor the names of its
C   contributors may be used to endorse or promote products derived
C   from this software without specific prior written permission.
C
C THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
C "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
C LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
C A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
C OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
C SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
C LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
C DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
C THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
C (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
C OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

C=======================================================================
      SUBROUTINE DBISTE (NDB, OPTION, ISTEP,
     &  NVARGL,
     *  NVARNP, NUMNP,
     *  NVAREL, NELBLK, NUMELB, ISEVOK, IDELB,
     *  NVARNS, NUMNPS, NNNPS,  ISNSOK, IDNPS,
     *  NVARSS, NUMESS, NEESS,  ISSSOK, IDESS,
     *  TIME, VARGL, VARNP, VAREL, VARNS, VARSS, *)
C=======================================================================
C     --*** DBISTE *** (EXOLIB) Read database variables for one time step
C     --   Written by Amy Gilkey - revised 10/14/87
C     --
C     --DBISTE reads the database global, nodal, and element variables
C     --for one time step.
C     --
C     --Parameters:
C     --   NDB - IN - the database number
C     --   OPTION - IN - ' ' to not store, '*' to store all, else store options:
C     --      'G' to store global variables
C     --      'E' to store element variables
C     --      'N' to store nodal variables
C     --      'M' to store nodeset variables
C     --      'S' to store sideset variables
C     --   ISTEP - IN - the time step number
C     --   NVARGL - IN - the number of global variables
C     --   NVARNP - IN - the number of nodal variables
C     --   NUMNP - IN - the number of nodes
C     --   NVAREL - IN - the number of element variables
C     --   NELBLK - IN - the number of element blocks
C     --   NUMELB - IN - the number of elements per block
C     --   ISEVOK - IN - the element block variable truth table;
C     --      variable i of block j exists iff ISEVOK(j,i)
C     --   TIME - OUT - the time step time
C     --   VARGL - OUT - the global variables for the time step (if OPTION)
C     --   VARNP - OUT - the nodal variables for the time step (if OPTION)
C     --   VAREL - OUT - the element variables for the time step (if OPTION)
C     --   * - return statement if error encountered, including end-of-file;
C     --      message is printed
C     --
C     --Database must be positioned in front of time step upon entry;
C     --upon exit positioned after time step.

      include 'exodusII.inc'
      INTEGER NDB
      CHARACTER*(*) OPTION
      INTEGER ISTEP
      INTEGER NVARGL
      INTEGER NVARNP, NUMNP
      INTEGER NVAREL, NELBLK, NUMELB(*), IDELB(*), ISEVOK(*)
      INTEGER NVARNS, NUMNPS, NNNPS(*),  IDNPS(*), ISNSOK(*)
      INTEGER NVARSS, NUMESS, NEESS(*),  IDESS(*), ISSSOK(*)
      REAL TIME
      REAL VARGL(*)
      REAL VARNP(*)
      REAL VAREL(*)
      REAL VARNS(*)
      REAL VARSS(*)

      CHARACTER*80 ERRMSG

C     --Read step time
      CALL EXGTIM (NDB, ISTEP, TIME, IERR)

C     --Read global variables

      IF ((OPTION .EQ. '*') .OR. (INDEX (OPTION, 'G') .GT. 0)) THEN
         if (nvargl .gt. 0) then
           call exggv (ndb, istep, nvargl, vargl, ierr)
           if (ierr .lt. 0) goto 180
         end if
      END IF

C     --Read nodal variables

      IF ((OPTION .EQ. '*') .OR. (INDEX (OPTION, 'N') .GT. 0)) THEN
         ioff = 1
         do 150 ivar = 1, nvarnp
           call exgnv (ndb, istep, ivar, numnp, varnp(ioff), ierr)
           if (ierr .lt. 0) goto 190
           ioff = ioff + numnp
 150     continue
      END IF

C     --Read element variables
      IF ((OPTION .EQ. '*') .OR. (INDEX (OPTION, 'E') .GT. 0)) THEN
         CALL DBIST2 (NDB, 'E', ISTEP, NVAREL, MAX(NVAREL,1), NELBLK,
     &        MAX(NELBLK,1), ISEVOK, VAREL, NUMELB,
     &        IDELB, IVAR, IELB, *200)
      END IF

C     --Read nodeset variables
      IF ((OPTION .EQ. '*') .OR. (INDEX (OPTION, 'M') .GT. 0)) THEN
         CALL DBIST2 (NDB, 'M', ISTEP, NVARNS, MAX(NVARNS,1), NUMNPS,
     &        MAX(NUMNPS,1), ISNSOK, VARNS, NNNPS,
     &        IDNPS, IVAR, IELB, *210)
      END IF

C     --Read sideset variables
      IF ((OPTION .EQ. '*') .OR. (INDEX (OPTION, 'S') .GT. 0)) THEN
         CALL DBIST2 (NDB, 'S', ISTEP, NVARSS, MAX(NVARSS,1), NUMESS,
     &        MAX(NUMESS,1), ISSSOK, VARSS, NEESS,
     &        IDESS, IVAR, IELB, *220)
      END IF

      RETURN

 180  CONTINUE
      WRITE (ERRMSG, '(A, I8)', IOSTAT=IDUM)
     &     'GLOBAL VARIABLES for TIME STEP', ISTEP
      GOTO 240

 190  CONTINUE
      WRITE (ERRMSG, '(A, I8, A, I8)', IOSTAT=IDUM)
     &     'NODAL VARIABLE', IVAR, ' for TIME STEP', ISTEP
      GOTO 240

 200  CONTINUE
      WRITE (ERRMSG, '(A, I8, A, I8, A, I8)', IOSTAT=IDUM)
     &     'ELEMENT VARIABLE', IVAR, ' of BLOCK', IELB,
     &     ' for TIME STEP', ISTEP
      GOTO 240

 210  CONTINUE
      WRITE (ERRMSG, '(A, I8, A, I8, A, I8)', IOSTAT=IDUM)
     &     'NODESET VARIABLE', IVAR, ' of NODESET', IELB,
     &     ' for TIME STEP', ISTEP
      GOTO 240

 220  CONTINUE
      WRITE (ERRMSG, '(A, I8, A, I8, A, I8)', IOSTAT=IDUM)
     &     'SIDESET VARIABLE', IVAR, ' of SIDESET', IELB,
     &     ' for TIME STEP', ISTEP
      GOTO 240

 240  CONTINUE
      CALL DBERR (IERR, ERRMSG)
      RETURN 1
      END

max_stars_repo_name: brian-kelley/seacas

max_issues_repo_name: brian-kelley/seacas

max_forks_repo_name: brian-kelley/seacas


text:
*DECK BESKS
      SUBROUTINE BESKS (XNU, X, NIN, BK)
C***BEGIN PROLOGUE  BESKS
C***PURPOSE  Compute a sequence of modified Bessel functions of the
C            third kind of fractional order.
C***LIBRARY   SLATEC (FNLIB)
C***CATEGORY  C10B3
C***TYPE      SINGLE PRECISION (BESKS-S, DBESKS-D)
C***KEYWORDS  FNLIB, FRACTIONAL ORDER, MODIFIED BESSEL FUNCTION,
C             SEQUENCE OF BESSEL FUNCTIONS, SPECIAL FUNCTIONS,
C             THIRD KIND
C***AUTHOR  Fullerton, W., (LANL)
C***DESCRIPTION
C
C BESKS computes a sequence of modified Bessel functions of the third
C kind of order XNU + I at X, where X .GT. 0, XNU lies in (-1,1),
C and I = 0, 1, ... , NIN - 1, if NIN is positive and I = 0, 1, ... ,
C NIN + 1, if NIN is negative.  On return, the vector BK(.) Contains
C the results at X for order starting at XNU.
C
C***REFERENCES  (NONE)
C***ROUTINES CALLED  BESKES, R1MACH, XERMSG
C***REVISION HISTORY  (YYMMDD)
C   770601  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   890531  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
C   900326  Removed duplicate information from DESCRIPTION section.
C           (WRB)
C***END PROLOGUE  BESKS
      DIMENSION BK(*)
      SAVE XMAX
      DATA XMAX / 0.0 /
C***FIRST EXECUTABLE STATEMENT  BESKS
      IF (XMAX.EQ.0.0) XMAX = -LOG (R1MACH(1))
C
      IF (X .GT. XMAX) CALL XERMSG ('SLATEC', 'BESKS',
     +   'X SO BIG BESSEL K UNDERFLOWS', 1, 2)
C
      CALL BESKES (XNU, X, NIN, BK)
C
      EXPXI = EXP (-X)
      N = ABS (NIN)
      DO 20 I=1,N
        BK(I) = EXPXI * BK(I)
 20   CONTINUE
C
      RETURN
      END


text:
submodule (Prng_Class) sm_xorshift128plus
!!# xorshift128+
!!This module contains routines to generate random numbers using the xorshift128+ method by Vigna's extensions to
!!Marsaglia, G., 2003. Xorshift RNGs. Journal of Statistical Software 8, 1-6.
!!
!!This module is a modified version of the public domain code written by Shun Sakuraba in 2015-2016.
!!The original code can be found [here](https://bitbucket.org/shun.sakuraba/xorshiftf90).
!!The modified code in coretran is distributed under the coretran license, see the repository for that information.
!!
!!Vigna's xorshift128+ pseudorandom generator.
!!Sebastiano Vigna. 2014. Further scramblings of Marsaglia's xorshift generators. CoRR, abs/1402.6246.
!! xorshift128+ is known as a fast pseudorandom generator with reasonably good resilience to randomness tests.
!! Since its primary imporance is its speed, do not forget to add inlining directives depending your compiler.
!!
!!
!!## Why xorshift128+ ?
!!The seed of the xorshift128+ can be jumped by k cycles, where each cycle jumps ahead by $2^{64}$ random numbers, but quickly.
!!This ability to jump is an important aspect for using random number generators in general but especially in parallel
!!in either OpenMP or MPI paradigms.
!!The time taken to generate a 64 bit integer was 1.06 ns on an IntelR CoreTM i7-4770 CPU @3.40GHz (Haswell) as shown in Vigna (2014).
!!xorshift128+ is the better for less massively parallel applications.
!!
!!If the period of the random number generator is too small for your application, consider using the xorshift1024*
!!generator which has a period of $2^{512}$ numbers. The time to generate a 64 bit integer is slightly slower at 1.36 ns
!!
!!### Important note: This next aspect needs to be tested.
!!@todo: test this
!!Because the Fortran standard does not specify an integer model, int64 overflow has an undefined behaviour.
!!Thus, there might be warnings issued while compiling this program, or worse, it may not work.
!!A workaround uses the int(int(a, kind=16) "op" b, kind=8) trick, which works well with gfortran,
!!but commercial compilers do not accept 128-bit arithmetics (incl. Intel, Cray, and Fujitsu).

implicit none

contains

  !====================================================================!
  module procedure rngInteger_128plus!(this, val)
  !====================================================================!
  !class(prng), intent(inout) :: this
  !integer(i64), intent(out) :: val

  integer(i64) :: s0, s1

  ! swap state variables
  s0 = this%seed(2)
  s1 = this%seed(1)

  this%seed(1) = s0

  ! xorshift
  s1 = ieor(ishft(s1, 23), s1)

  ! scramble
  s1 = ieor(ieor(s1, s0), ieor(ishft(s1, -18), ishft(s0, -5)))
  this%seed(2) = s1

  val = s1 + s0
  end procedure
  !====================================================================!
  !====================================================================!
  module procedure jumpState_128plus!(this)
  !====================================================================!
  !! Jumps the current prng state by $2^{512}$ pulls from the generator without drawing $2^{512}$ numbers.
  !class(Prng), intent(inout) :: this
  integer(i64) :: jump_table(2)
  integer(i64) :: t(0:1)
  integer(i64) :: tmp
  integer :: i, j, b
  data jump_table / -8476663413540573697_i64,& !0x8a5cd789635d2dff
                     1305993406145048470_i64/  !0x121fd2155c472f96

  t(:) = 0

  do i = 1, 2
      do b = 0, 63
        if(btest(jump_table(i), b)) then
            do j = 0, 1
              t(j) = ieor(t(j), this%seed(j))
            end do
        end if
        call rngInteger_128plus(this, tmp)
      end do
  end do
  this%seed(0:1) = t
  end procedure
  !====================================================================!
end submodule

max_stars_repo_name: leonfoks/coretran

max_issues_repo_name: leonfoks/coretran

max_forks_repo_name: leonfoks/coretran